Selection on grain shattering genes and rates of rice domestication

Plant material

We sampled 31 accessions of cultivated rice (O. sativa L.), including 17 accessions of O. sativa ssp. indica and 14 accessions of O. sativa ssp. japonica. This sample represents all five recognized subgroups within O. sativa (Garris et al., 2005). A total of 44 accessions of the wild progenitors were sampled, including 28 accessions of O. rufipogon and 16 accessions of O. nivara (Table 1). These collections covered a wide distributional range of the wild species from southeastern Asia to India, including all potential areas for rice domestication (see the Supporting Information, Fig. S1). Of the samples, 30 cultivars and 28 wild accessions were previously studied by Zhu et al. (2007). Oryza barthii was used as an outgroup for phylogenetic analyses. Seeds of the sampled accessions were germinated and DNA was isolated from seedlings.

DNA sequencing

A c. 2700-bp region of the sh4 gene was amplified using PCR primers InF 5′-ACCAGCTCAACTCAACAACG and BR 5′-GTGTCAAAGCCTTAATGTGC. This region includes the entire first exon and c. 1900 bp of 5′ upstream region of the gene. For qSH1, we amplified a c. 1600-bp region surrounding the previously reported functional single-nucleotide polymorphism (SNP) using primers q_rF1 5′-GGCTACGAGCGATAAATCTG and q_rR2 5′-GCTACAGGCTCCTACTATAC. The PCR and internal sequencing primers were designed based on rice genome sequences (GenBank accession AP008210 and AL606619 for sh4 and AB071330 and AB071331 for qSH1). The sh4 region was amplified with Takara LA Taq DNA polymerase with GC buffer I (Takara, Shiga, Japan). The qSH1 region was supplemented with ExTaq DNA polymerase (Takara). The PCR products were purified with either a Pharmacia purification kit (Amersham Pharmacia Biotech) or a Dinggou purification kit (Dingguo, Beijing, China). For cultivated rice, purified PCR products were sequenced directly on both strands. For all individuals of wild species O. rufipogon and O. nivara, PCR products were cloned into pGEM T-easy vectors (Promega). At least six clones were sequenced initially to obtain both alleles from the heterozygous individuals at a locus. Singletons were confirmed by reamplification and resequencing. All sequences were deposited into GenBank, with accession numbers EU999788–EU999948.

Sequence analyses and neutrality tests

Sequence alignments were performed using clustalx (Thompson et al., 1997) and were refined manually. Genetic variation of each gene region was estimated with average pairwise differences per base pair between sequences (π) (Nei & Li, 1979) and Watterson's estimate (θ_w) based on the number of segregating sites (Watterson, 1975) using DnaSP version 4.10 (Rozas et al., 2003). Under the standard neutral model for an autosomal gene from a random-mating population of a constant size, both π and θ_w are expected to be equal to 4N_eµ, where N_e represents the effective population size and µ represents the mutation rate per generation per site (Watterson, 1975).

Selection on the shattering genes was tested using multiple methods. Departure from neutrality was tested at segregating nucleotide sites with Tajima's D and Fu and Li's D* and F* using the program dnasp (Rozas et al., 2003). Tajima's D measures the difference between the mean pairwise differences (π) and Watterson's estimator (θ_w) (Tajima, 1989). Fu and Li's D* and F* test the discrepancy between the number of polymorphic sites in external branches, or polymorphisms unique to a sequence, and number of polymorphic sites in internal phylogenetic branches, or polymorphisms shared by the sequences (Fu & Li, 1993).

We also performed maximum likelihood Hudson–Kreitman–Aguade (MLHKA) test (Wright & Charlesworth, 2004). The maximum likelihood ratio test assesses departure from neutrality at a focal locus compared with neutral standards. The degree of increase or decrease of polymorphism caused by selection is measured as k. Here seven neutrally evolving rice genes (Adh1, GBSSII, Ks1, Lhs1, Os0053, SSII1, TFIIAγ-1; Zhu et al., 2007) were used as reference sequences and O. barthii was used as the outgroup for the MLHKA tests. The test was performed using the program MLHKA provided by S. I. Wright (http://www.yorku.ca/stephenw/). With different random number seeds, Markov chain lengths of 100 000 were run in the MLHKA tests to specify each gene in each accession, and at least three independent runs per model were performed to assess the maximum likelihood values.

To further characterize selection on sh4 and qSH1 in the context of the demographic dynamics of rice domestication, we conducted a coalescent simulation (CS) as detailed in Tenaillon et al. (2004). This evaluates whether the domestication bottleneck alone could account for the reduction of nucleotide diversity of cultivated rice in comparison with its wild progenitors. In CS analyses, the statistic is a likelihood ratio of the best-fitting demographic model for the focal genes (sh4 and qSH1) against an estimate of the demographic history of neutrally evolving genes, in this case Adh1, GBSSII, Ks1, Lhs1, Os0053, SSII1 and TFIIAγ-1. We used parameter estimates for the seven neutral genes from Zhu et al. (2007). To estimate demographic parameters for sh4 and qSH1 we assumed that rice cultivation began 10 000 yr ago, and that the domestication bottleneck spanned 3000 generations, consistent with archaeological evidence (Khush, 1997). Given this model, we explored the bottleneck severity K (Wright et al., 2005) over a grid of 25 values (0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1). The number of segregating sites (S) in cultivated rice was used as the goodness-of-fit statistic for all simulations. Results for sh4 and qSH1 were based on 10 000 simulations of the best-fitting coalescent model.

Phylogenetic analyses

The phylogeny of sequences was inferred by Neighbor-joining (NJ), maximum likelihood (ML), maximum parsimony (MP) and Bayesian inference (BI) methods. Neighbor-joining, ML and MP were conducted with paup 4.0b10 (Swofford, 2002). Bootstrap analyses were performed with 1000 replicates. Bayesian analyses were conducted using mrbayes version 3.1 (Huelsenbeck & Ronquist, 2001; Ronquist & Huelsenbeck, 2003) with Markov chain Monte Carlo (MCMC) estimation of posterior probability distributions. Four chains of the MCMC were run each for 1 000 000 generations and were sampled. An initial ‘burn-in’ phase of 100 000 generations was discarded. We then saved one tree every 100 generations and constructed a majority-rule consensus tree of the 9000 trees sampled in each run. The appropriate nucleotide substitution model for each data set was determined by modeltest 3.7 (Posada & Crandall, 1998). Optimization of the model and parameter search was made by using of Akaike information criterion (AIC, Sakamoto et al., 1986).

Estimation of the time since the fixation of the sh4 allele

We used two different approaches to estimate the time required to fix the nonshattering sh4 allele in cultivated rice. Based on the initial frequency of nonshattering sh4 allele (p = 1/2N, where N is the assumed population size) and selection coefficient (s), we first estimated the time required for sh4 fixation (T_f) using equations from Kimura & Ohta (1969):

where S = Ns, and .

To calculate T_f, s was estimated using the relationship d_phys = 0.01 s/c (Kaplan et al., 1989), where the physical distance of a selective sweep (d_phys) was obtained from the aligned genome sequences of rice chromosome 4 between an indica (GLA4) and a japonica (Nipponbare) cultivar and the local recombination rate c was calculated as described in Zhao et al. (2002).

Second, we used an approximate Bayesian method to investigate a simple model of a selective sweep in a bottlenecked, derived population (Thornton & Jensen, 2007). We first simulated a deterministic trajectory of the selected site forward in time, then a structured coalescent backward in time, conditional on the simulated allelic trajectory. For both indica and japonica, we simulated an ancestral population of size N_A =N₁ + fN₁ which splits at time t in the past into a modern wild population of size N₁ and a bottlenecked domesticated population of size fN₁ (see the Supporting Information, Fig. S2). At time t-d the bottlenecked population recovered to size N₁ (Fig. S2). We simulated a selective sweep in the bottlenecked populations starting shortly after the bottleneck began and ending at time t-d. We assumed t = 10 000 yr and N₁ = 250 000, in close agreement with archaeological data (Zhao & Piperno, 2000; Zong et al., 2007) and diversity at neutral loci (Zhu et al., 2007). Model parameters for each simulation were drawn from previous probability distributions, with f ∼ U (0.001, 0.5) and s ∼ U (0.0001, 1). The bottleneck duration d was set to be slightly longer than the expected duration of the selective sweep given f and s, using the equation d = f × ln(2fN₁) × (2[fN₁ + 0.5]s) + 10⁻⁵ (modified from Thornton & Jensen (2007) and Stephan et al. (1992)). The population mutation rate θ in each simulation was taken from previous estimates of nucleotide diversity (Zhu et al., 2007), and the population recombination rate ρ was set equal to θ. We assumed that the population recombination rate, calculated as ρ = 4N₁c, where c is the recombination rate per bp per generation, was equal to θ. We simulated a locus of 1800 silent sites and a number of alleles for each population identical to the observed sample sizes at sh4.

For each simulated dataset we calculated the number of segregating sites (S) and nucleotide diversity (π) at silent sites in both domesticated and wild population. Datasets with values of all four summary statistics within 20% of the observed values at locus sh4 were retained, and the selection coefficient (s), duration of the bottleneck (d), and severity of the bottleneck (f) were recorded. Recorded values of each parameter were used to estimate acceptance rates and form the posterior distribution.

Molecular diversity of sh4 and qSH1

In total we obtained sh4 sequences from 30 accessions of cultivated rice and 41 accessions of the wild progenitors, and qSH1 sequences from 31 accessions of rice cultivars and 32 accessions of the wild species (Tables 1, 2). These sequences were aligned with the outgroup, O. barthii, generating alignments of 2740 bp and 1575 bp for sh4 and qSH1, respectively.

Table 2. Estimates of nucleotide variation

Locus	Taxon	N a		H c	S d
sh4	Oryza sativa ssp. indica	34	2450.5	4	3	0.0003	0.0003	0.0004	0.0004	0
	O. sativa ssp. japonica	26	2450.3	3	3	0.0004	0.0003	0.0006	0.0004	0
	Oryza rufipogon	50	2454.1	32	64	0.0042	0.0060	0.0053	0.0080	8
	Oryza nivara	32	2449.6	20	45	0.0038	0.0046	0.0049	0.0060	8
qSH1	O. sativa ssp. indica	34	1542.0	12	21	0.0032	0.0033	0.0032	0.0033	3
	O. sativa ssp. japonica	28	1540.2	2	2	0.0005	0.0003	0.0005	0.0003	0
	O. rufipogon	38	1536.3	22	45	0.0054	0.0076	0.0054	0.0076	4
	O. nivara	26	1545.2	12	33	0.0057	0.0056	0.0057	0.0056	1

• a Total number of sequences.
• b Average length (bp) of the sequences from each species.
• c Total number of haplotypes at an individual locus in each taxon.
• d Total number of polymorphic sites.
• e Average number of pairwise nucleotide difference per site calculated based on the total number of polymorphic sites.
• ^f Watterson's estimator of θ per base pair calculated based on the total number of polymorphic sites.
• g Average number of pairwise nucleotide differences per site calculated based on silent sites.
• ^h Watterson's estimator of θ per base pair calculated based on silent sites.
• ⁱ R _M, the minimum number of recombination events.

It has been shown that the G to T substitution of sh4, which results in an amino acid substitution from lysine to asparagine, was responsible for the reduction of grain shattering during the domestication of rice (Li et al., 2006b). At the functional SNP site, all of the alleles from cultivars have a thymine (T) residue, whereas almost all of those from the wild accessions, including the outgroup, have a guanine (G) residue (Table 1). Two accessions of O. nivara, one from Sri Lanka (niv-LKA1) and the other from Cambodia (niv-KHM), are T/G heterozygotes, which is most likely the result of introgression of the nonshattering allele from cultivars into the wild populations (Li et al., 2006b).

Phylogenetic analysis of the sh4 dataset showed that all of the cultivars with T at the functional SNP site were grouped together with 65% bootstrap support (Fig. 1a). The alleles with T at the functional SNP site from the two heterozygous individuals of O. nivara were nested within cultivars. Thus, this monophyletic group of sh4 sequences with T at the functional SNP site contains all nonshattering alleles derived from domestication. Within this group, indica and japonica cultivars are intermixed. Immediately outside of this clade, two wild accessions, an O. nivara from India and an O. rufipogon from Laos, form a sister group to the cultivated alleles with 74% bootstrap support. Both of the alleles in the sister group contain a G at the functional SNP site.

At the qSH1 locus, three temperate japonica individuals were found homozygous T/T at the SNP site that confers the reduced shattering phenotype. The remaining accessions of cultivars, including two temperate japonica individuals, and all wild accessions, were homozygous G/G, conferring the shattering phenotype (Table 1; Fig. 1b). The three temperate japonica accessions with reduced shattering genotype form a monophyletic group with 88% bootstrap support. This group is nested in a clade that contains the rest of the japonica cultivars together with two indica, two O. nivara and one O. rufipogon accession. The remaining indica cultivars are intermixed with O. nivara and O. rufipogon in a well-supported clade.

With regard to genetic variation at the shattering loci, cultivars have lower levels of nucleotide diversity at both loci than the wild progenitors. However, the degree of reduction in nucleotide diversity differs between the two loci and between the cultivars. For sh4, the cultivars have the same level of nucleotide diversity (θ_sil = 0.0004 for both indica and japonica), which is 15- to 20-times lower than that of O. nivara or O. rufipogon (Table 2). At the qSH1 locus, however, the similar level of reduction in nucleotide diversity occurs only in japonica. The indica accessions contain about half of the nucleotide diversity of the wild species (Table 2).

Neutrality tests

We conducted several neutrality tests to determine whether the reduction in nucleotide diversity in rice cultivars could be caused by artificial selection. The MLHKA test, performed in reference to seven neutral genes (Zhu et al., 2007), indicated that a selective sweep could be the primary cause of the severe reduction in nucleotide polymorphism in both indica and japonica cultivars at the sh4 locus (P < 0.0001; Table 3). Because this could reflect selection on sh4 in wild populations (Yamasaki et al., 2005), we also applied MLHKA to the O. nivara and O. rufipogon samples (Table 3); the test of neutrality is borderline significant (P = 0.042) for O. nivara but not significant for O. rufipogon (P = 0.193). The results thus indicate that there was strong artificial selection on sh4 during rice domestication. At the qSH1 locus, the MLHKA test was marginally significant in japonica (P = 0.046) but not significant for either indica or the wild species (Table 3).

Unlike the HKA test, Tajima's D and Fu and Li's D* and F* rejected neutrality in only two instances of 24 statistics obtained (Tajima, 1989; Fu & Li, 1993; Table 3). These were Fu and Li's D* for qSH1 of indica and O. nivara (P < 0.01). However, those nonsignificant statistics are likely to be an artifact of population demographics (Tajima, 1989; Nielsen, 2001), and do not necessarily support neutrality. After populations undergo subdivision and bottlenecks (Garris et al., 2005; Zhu et al., 2007) there is an excess of alleles of intermediate frequency that leads to higher Tajima's D values (Das et al., 2004). In addition, theoretical work (Innan & Kim, 2004) has shown that the multilocus HKA test is more powerful in detecting artificial selection than the single-locus Tajima's D when polymorphisms is much reduced.

To determine whether the low nucleotide diversity at sh4 and qSH1 could be explained by a genetic bottleneck alone, we conducted the CS test in reference to seven neutrally evolving genes. We began by inferring a bottleneck model with the seven neutral genes, and then compared polymorphism in sh4 and qSH1 with the expectation of the model (Table 4). The results indicated that variation detected at sh4 could not be explained by a domestication bottleneck alone for indica (P = 0.016), japonica (P = 0.047) or a combined sample from both cultivars (P = 0.002). By contrast, variation at qSH1 was not significantly different from the expectation under a bottleneck scenario for indica (P = 0.939), japonica (P = 0.395) or the combined sample (P = 0.930). Together, the results of the MLHKA and CS tests suggest that sh4 was fixed in cultivated rice under artificial selection. The extent to which qSH1 was selected in japonica rice remains unclear.

Timing of sh4 fixation under artificial selection

Taking advantage of the wealth of genomic information for rice chromosome 4, we estimated the time required to drive the nonshattering allele to fixation by artificial selection. We first calculated the selection coefficient (s) using the equation d_phys = 0.01 s/c (Kaplan et al., 1989), where c is the recombination rate and d_phys is the physical distance between a selected and hitchhiking neutral site. Based on aligned genome sequences of rice chromosome 4 between an indica (GLA4) and a japonica (Nipponbare) cultivar, we observed hypervariable regions 21 kb upstream and 43 kb downstream of sh4. In this c. 66 kb region containing sh4, there were 34 SNPs between the indica and japonica cultivars. The SNP density of 0.52 SNPs per kb around the sh4 locus thus is much lower than 3.5 SNPs per kb for the entire chromosome 4 (Feng et al., 2002; B. Han, pers. comm.). This suggests that the selective sweep spanned a 66 kb region, thus signifying a 33 kb mean distance between the selected site and the farthest sweep site.

Given the local recombination rate around the sh4 locus of approx. 1.122 × 10⁻⁷/bp (B. Han, Pers. Comm.), s is calculated as 0.37. Using the previously reported nucleotide diversity at silent sites (θ_sil = 0.0024–0.00292) for O. sativa (Caicedo et al., 2007; Zhu et al., 2007) and a substitution rate of 5.9 × 10⁻⁹ substitutions per synonymous site per year for grasses (Gaut et al., 1996; White & Doebley, 1999), we estimated that the effective population size ranged from 10 200–12 300. To be cautious, we used a much wider range of the effective population size of O. sativa, between 1000 and 100 000, for the calculation of fixation time. With estimates of s and N_e, we calculated the fixation time (T_f) of sh4 according to Kimura & Ohta (1969). Assuming the initial allele frequency of 1/2N, T_f for sh4 was estimated to be between 76 and 112 yr.

Although our estimate of s was based on detailed genome sequence data, we still could not rule out potential bias associated with the calculation. Thus, a simulation method was also employed for estimating fixation time. Simulations incorporating a demographic model like that of Thornton & Jensen (2007) generate similar insights about the timing of the selective sweep. The joint posterior probability distributions of the bottleneck duration d and the selective coefficient s are shown in Fig. 2. The marginal posterior distributions for both parameters in both subspecies are leptokurtic, with the highest probability assigned to very low values but with long, relatively fat tails. The medians of the posterior distribution of f (japonica 0.034 and indica 0.043) and s (japonica 0.187 and indica 0.235) point to a dramatic reduction in population size and relatively strong selection. Under this model (see the Materials and Methods section), the duration of the bottleneck (d) is closely related to the length of time to fixation of the selective alleles. The median of the posterior distribution of d (japonica 118 yr, indica 97 yr) suggests that the fixation of sh4 could have occurred in an extremely short period of time, which is concordant with the earlier calculation based on estimated s. Estimates of the duration based on the joint posterior distributions of s and f are similarly rapid, whether based on median values (japonica 86 yr and indica 72 yr) or the peak (japonica 308 yr and indica 280 yr) of the distribution. This suggests that the fixation of sh4 could have occurred in a short period of time ranging from < 100 yr to only a few hundred years.