# data[].attributes.headline
How to read a scientific paper, part I: structure of a paper

# data[].attributes.subject
How to read a scientific paper, part I: structure of a paper

# data[].attributes.headline
How to read a scientific paper, part II: case study

# data[].attributes.subject
How to read a scientific paper, part II: case study

# data[].attributes.headline
There are no more spots for working on this article. 

# data[].attributes.subject
There are no more spots for working on this article. 

# data[].attributes.body
--To talk about the paper/presentation my email is: kiera.thurman@ucd.connect.ie
fionn.cahalane@ucdconnect.ie

# data[].attributes.body
A haplotype map of allohexaploid wheat reveals
distinct patterns of selection on homoeologous
genomes

Please put a special emphasize on exome sequencing

Presenters:Timothy Reed, Harshini M, Sarah Aldhafiri

# data[].attributes.attachment_link.readable_text
Abstract

Background

Bread wheat is an allopolyploid species with a large, highly repetitive genome. To investigate the impact of selection on variants distributed among homoeologous wheat genomes and to build a foundation for understanding genotype-phenotype relationships, we performed population-scale re-sequencing of a diverse panel of wheat lines.

Results

A sample of 62 diverse lines was re-sequenced using the whole exome capture and genotyping-by-sequencing approaches. We describe the allele frequency, functional significance, and chromosomal distribution of 1.57 million single nucleotide polymorphisms and 161,719 small indels. Our results suggest that duplicated homoeologous genes are under purifying selection. We find contrasting patterns of variation and inter-variant associations among wheat genomes; this, in addition to demographic factors, could be explained by differences in the effect of directional selection on duplicated homoeologs. Only a small fraction of the homoeologous regions harboring selected variants overlapped among the wheat genomes in any given wheat line. These selected regions are enriched for loci associated with agronomic traits detected in genome-wide association studies.

Conclusions

Evidence suggests that directional selection in allopolyploids rarely acted on multiple parallel advantageous mutations across homoeologous regions, likely indicating that a fitness benefit could be obtained by a mutation at any one of the homoeologs. Additional advantageous variants in other homoelogs probably either contributed little benefit, or were unavailable in populations subjected to directional selection. We hypothesize that allopolyploidy may have increased the likelihood of beneficial allele recovery by broadening the set of possible selection targets.

Background

Wheat genomic variation is shaped by the interplay of multiple factors including two recent polyploidization events [-] (Figure a), domestication [], spread from the sites of origin to new geographic regions, gene flow from the populations of wild and domesticated ancestors [], and post-domestication selection aimed at developing high-yielding locally adapted varieties. The eco-geographic habitats to which wheat is adapted span diverse environments ranging from low humidity regions in Nigeria, and the northern regions of Russia and Norway to the high-humidity regions of South America and Bangladesh []. It has been suggested that this broad adaptability likely results from the genetic diversity captured from the natural populations of its tetraploid ancestors [,] combined with a high rate of evolutionary changes in the wheat genome (particularly insertions and deletions), which are tolerated by its polyploid nature [,].

A detailed description of DNA sequence variation across the genome is a prerequisite for the systematic analysis of variants underlying trait variation in wheat and critical for understanding the role of various evolutionary factors in shaping genome diversity. Recently, low to medium density genotyping arrays were used to characterize SNP variation and linkage disequilibrium (LD) in wheat populations [,] and identify variants associated with phenotypic traits [,]. However, despite being a useful genotyping tool, these arrays are incapable of capturing the entire spectrum of DNA sequence variation and allele frequencies in wheat populations, and providing unbiased information that may help directly identify causal variants affecting phenotypes. Achieving this goal requires obtaining sequence data on a genome-scale from a diverse population of lines. To date, this has only been performed on limited samples of wheat lines used to discover SNPs in the parental lines of mapping populations [], or for SNP-based array design [].

Genome sequencing of populations of individuals has been undertaken in a number of species including humans, Arabidopsis, and several crops [-] and helped to detect alleles contributing to phenotypic variation and adaptation. Despite recent advances in next-generation sequencing (NGS), performing similar analyses of genomic variation in wheat is substantially complicated by allopolyploidy and large genome size (approximately 17 Gb). However, sequencing of DNA samples subjected to complexity reduction by exome capture [,] and genotyping by sequencing [,] was shown to be an effective strategy to analyze the complex genomes. In addition, the recent release of the chromosome-specific wheat genome assemblies [,] can now help to alleviate the problems associated with variant calling in the allopolyploid genome, and allow us to describe the chromosomal distribution of variants and their potential effect on gene function.

Here we used the newly developed genome assembly of the cultivar Chinese Spring [] based on flow-sorted chromosome survey sequence (CSS) contigs to create a diversity map of allopolyploid bread wheat. The data were generated by re-sequencing 62 diverse wheat lines using whole exome capture (WEC) and genotyping by sequencing (GBS) approaches. The panel of wheat lines was selected to capture the genetic diversity of the major global wheat growing regions and included landraces and cultivars (Figure b; Table S1 and Figure S1 in Additional file ). We used the obtained data to describe the effect of genetic variation on gene function, gain insights into the effect of selection on duplicated genes, and explore the LD landscape in each of the three wheat sub-genomes to better understand the role of selection in shaping the genetic diversity of wheat.

Results and discussion

Re-sequencing the allopolyploid wheat genome

The WEC assay probes were designed with 107 Mb of non-redundant low-copy genic regions [] targeting nearly 321 Mb of sequence in all three wheat genomes (Figure S2 in Additional file ). The capture probes covered 78% of the 124,201 high-confidence protein-coding genes (at the 95% similarity threshold) in the CSS contigs []. The GBS approach generated sequence data primarily outside the genic regions. We produced roughly 4.7 billion paired-end reads (4.5 billion WEC reads and 0.2 billion GBS reads), and 62% of WEC reads and 51% of GBS reads uniquely mapped to the CSS contigs of the individual chromosomes (Figure S3 and Tables S2, S3 in Additional file ), using alignment parameters optimized to separate reads from the different wheat genomes (Figures S4 and S5 in Additional file ). In the WEC dataset, similar relative read coverage across homoeologous targets indicated that non-redundant capture probes are capable of recovering sequences from the different genomes with equal efficiency (Figure c, Table S4 in Additional file ).

The overall genotype error rates for SNPs and indels assessed by comparing genotype calls generated for cultivar Chinese Spring with the CSS contigs of the same cultivar were 1.1% and 1.5%, respectively (for details see Materials and Methods). The error rate for rare SNPs and indels covered by >10 reads was 4.6% and 3.4%, respectively (Figure S7 in Additional file ). The majority (77%) of variants in the GBS dataset were found in intergenic regions (Table S5 in Additional file ), and only 4.3% of the variants overlapped with the WEC dataset (Figure d). This finding is supported by the in silico PstI digest of the CSS contigs, which showed that the regions targeted by the GBS cover only 6.8% of the regions targeted by the WEC.

Impact of purifying selection on genetic variation in the polyploid genome

In coding sequences we identified 83,622 non-synonymous and 76,361 synonymous SNPs (Table S5 in Additional file ). Based on high-confidence gene models in the CSS contigs, we determined that only 1,600 and 1,583 SNPs are predicted to produce premature termination codons (PTCs) and splice-site disruptions (SSDs), respectively, with a two- to three-fold lower incidence of functional mutations in the D genome than in the A and B genomes (Table S6 in Additional file ). Out of the 6,230 genes that have homoeologous copies in the wheat genome and harbor coding sequence-disrupting mutations including frame-shift indels and PTCs, 4,870 (78%) have at least one intact homoeologous copy suggesting that the deleterious effects of these variants, if any, could be compensated. The compensatory potential of the duplicated homoeologous genes is consistent with a higher density of chemically induced mutations (five- to eight-fold) that can be achieved in wheat compared to other diploid plant species []. However, in spite of the presence of intact functional homoeologous copies of genes, we found a reduced number of non-synonymous, PTC, and SDS variants with a high derived allele frequency in the population (Figure e; Figure S9 in Additional file ). This depletion of functional variants at higher allele frequencies is consistent with purifying selection acting against functional mutations, and suggests that the effect of purifying selection is not completely diminished by whole genome duplication. There are two plausible explanations for the retention of functional gene copies in young polyploids. First, selection acts on gene function that was partitioned among the homoeologs after the whole genome duplication. The functional partitioning is consistent with the studies of natural and artificial polyploids of wheat and other plants, often showing the tissue- and/or development-specific expression of homoeologs [-]. Alternatively, selection favors functional homoeologs to maintain the optimal stoichiometric ratios of gene products in macromolecular complexes, or in multistep regulatory cascades, which was proposed in the gene balance hypothesis [].

The fraction of non-synonymous changes varied significantly among different protein families (Table S8 in Additional file ) with the majority of PFAM domains involved in basic cellular functions showing a reduced proportion of non-synonymous changes compared to the genome-wide value, indicative of strong purifying selection. We detected a significant enrichment ( χ2 test P value <10⁻⁴, Table S8 in Additional file ) for non-synonymous changes in the LRR and NB-ARC domains of disease resistance genes. The enrichment for major effect SNPs in these genes appears to be common for plant genomes and was also found in Arabidopsis [] and peanut []. These observations are consistent with the hypothesis of an 'arms race' between the evolving populations of a pathogen and a plant defense system that results in fast evolution of genes with new disease-resistance specificities [].

Global patterns of genetic variation

The global patterns of genomic variation and distribution of inter-variant associations are impacted by historic selection and demographic events, and by variation in recombination rate []. We found a non-random variant distribution along the chromosomes with reduced variation near the centromeres and elevated variation at the telomeres (Figure a and b; Figures S10-15 in Additional file ), which is consistent with previous studies [,]. This pattern is similar to what was reported for maize and humans [,], but differs from Arabidopsis [], where regions of high polymorphism were located near the centromeres. Our data also showed reduced diversity and an excess of rare alleles in the D genome when compared to the A and B genomes (Figure a; Table S9 in Additional file ) []. These trends are consistent with the hypothesis that the limited number of ancestral genotypes of the D genome contributed to the origin of hexaploid wheat []. An elevated level of diversity in the A and B genomes, which otherwise would be expected to show the same levels of diversity as the D genome, could be attributed to the influx of allelic variation from the sympatric populations of wild tetraploid relatives [,].

Differentiation between landraces and cultivars ( F ST) varied along chromosomes with lower values found near the telomeres (Figure b; Figures S10-15 in Additional file ). Long stretches of elevated F ST were found along chromosome 4A and the short arm of chromosome 7B, two of the most structurally re-arranged chromosomes in the wheat genome [] (Figure S16 in Additional file ). Since these structural re-arrangements are fixed in wheat and, therefore, unlikely to affect gene flow between populations resulting in high F ST, the overlap of differentiated genomic regions with those showing the signal of positive selection (Table S10 in Additional file ) suggests that the detected differentiation could be associated with improvement selection.

Identification of genomic regions showing positive Tajima's D (excess of common alleles) and elevated diversity in the D genome (Figure b; Figures S10-15 in Additional file ) is indicative of the presence of highly diverged haplotypes in our sample most likely resulting from introgressions. To test this possibility, we compared the D genome haplotypes of 26 landraces not affected by modern breeding with the sequence of Ae. tauschii accession AL8/78 [], which is considered the most closely related to the wheat D genome []. The high proportion of AL8/78 alleles (74%) and their distribution along the wheat chromosomes confirms the ancestry of the D genome (Figure c) []. However, the fine-scale haplotypic structure also reveals regions carrying highly divergent haplotypes (Figure c) suggestive of significant levels of introgression from either the diverged Ae. tauschii lines or independently originated hexaploid wheat lineages founded by diverged Ae. tauschii genotypes. The preferential localization of introgressions in the high-recombining regions of the chromosomes indicates that gene flow between the different D genome lineages was uneven along the chromosomes.

Using re-sequencing data, we now have the possibility to assess historic recombination rate (parameter ρ) [], which is the product of meiotic recombination rate variation and effective population size. The median estimate of ρ in the D genome (ρ = 1.7 × 10⁻⁴/kb) was lower than in the A (ρ = 2.9 × 10⁻⁴/kb) and B (ρ = 3.2 × 10⁻⁴/kb) genomes (Figure a; Table S9 in Additional file ). Consistent with the observations made in Arabidopsis, maize, and humans [,,], we detected a positive correlation in the A (Spearman r ² = 0.23, P <10⁻⁹) and B ( r ² = 0.21, P <10⁻¹²) genomes between the meiotic (R = genetic distance in cM / physical distance in Mb) and historic recombination rates suggesting stability of chromosomal recombination rates over time. However, no significant correlation was found between R and ρ in the D genome. This fact is most likely explained by the polyploidy-associated population bottleneck, which can impact the estimates of historic recombination [] and also can result in reduced ρ in the D genome compared to that in the A and B genomes.

The historic recombination and diversity in wheat showed a positive correlation with relative distance from the centromere ( r  = 0.15, P <10⁻⁴), a trend previously reported in maize, and humans [,]. These relationships likely explain most of the inter-genomic diversity correlation along the homoeologous chromosomes (Table S11 in Additional file ). Except for a few cases, we found low inter-genomic correlation of the window-based Tajima's D and F ST estimates along the duplicated genomic regions of homoeologous chromosomes. This outcome is most likely a consequence of the sensitivity of these summary statistics to historic demographic and selection events that likely had different impacts on each wheat genome (Table S11 in Additional file ).

Genotype imputation and GWAS

In genome-wide association studies (GWAS), marker density affects the probability of finding variants in linkage disequilibrium (LD) with a causal variant. A SNP-hiding test [] showed that the probability of identifying high-LD SNPs ( r ²  > 0.8) in our population within a 2-kb window was 70% to 71% for all three genomes. Consistent with the observed LD levels in the wheat genomes (Figure a), for SNPs located from 2 kb to 4 kb apart, the probability of finding high-LD SNPs in the D genome (56%) was higher than in the A and B genomes (48% to 49%).

First, we sequentially selected each cultivar from our panel of 62 lines and, after 'hiding' all SNP sites besides those overlapping with the public 90 K SNP array [], we used the WEC SNP data in the remaining 61 lines to predict the 'hidden' variants. Depending on the selected wheat line, using a genotype calling probability cutoff of 0.6, the accuracy of genotype predictions assessed by comparing with the observed data was in the range of 93% to 97% (Figure a; Table S12 in Additional file ). This genotype probability cutoff value resulted in the removal of 5% to 15% of the data (Figure a), and allowed imputation of up to 549,918 SNPs. The accuracy of SNP imputation varied among the wheat genomes reflecting the inter-genomic differences in the extent of LD (Figure a). For example, the highest imputation accuracy was achieved for the D genome, which also showed the highest levels of inter-variant LD.

Second, we used our reference panel of 62 accessions to impute DNA polymorphisms in a GWAS. A panel of 678 diverse wheat landraces phenotyped for resistance to three rust diseases (Tables S13-S16 in Additional files , , , and ) [] was tested for marker-trait associations using genotyping data generated with the wheat 90 K SNP array. In this panel we were able to impute 344,544 SNPs, of which 210,017 SNPs with a MAF above 3% and the proportion of missing data less than 80% were used for GWAS. Three selected marker-trait associations were validated by mapping in the populations of recombinant inbred lines; two of these associations correspond to disease resistance loci Lr67 and Yr51, and one represents an uncharacterized locus (Figure b) associated with stem rust resistance (Tables S13-S16 in Additional files , , , and ). Two additional regions from our GWAS were shown to overlap with the positions of the previously mapped Lr37/Yr17/Sr38 and Lr68 disease resistance loci (Figure b). The markers associated with these loci showed highest similarity to CSS contigs 2AS-5264433 and 7BL-6748067, respectively [,]. Overall, comparison of marker-trait associations at non-imputed and imputed sites shows that imputed SNPs not only increase marker density but in most cases perform similar to or better than the SNPs directly genotyped using the 90 K assay (Figure b; Figure S17 in Additional file ). These results demonstrate the value of having a more complete ascertainment of DNA polymorphisms for GWAS that is achieved utilizing the high-density SNP variation data developed from 62 lines and the public 90 K SNP genotyping array [].

Signatures of selection in the polyploid genome

During the development of adapted lines, selection imposed by humans favored alleles controlling traits valuable for agriculture. When selection increases the frequency of beneficial alleles in a population, it impacts the standing variation of surrounding genomic regions resulting in reduced diversity, extended linkage disequilibrium, or strong inter-population allele frequency differentiation []. To detect these local patterns of variation, also referred to as 'selective sweeps', we investigated the patterns of genetic variation along the chromosomes and used two complementary approaches based on cross-population composite likelihood ratio (XP-CLR) [] (Figure S18 in Additional file ) and pair-wise haplotype sharing (PHS) [] tests (Tables S17-S19 in Additional files , , and ).

The reduced diversity observed near the wheat domestication genes Q and Tg was consistent with selection at the early stages of domestication (Figure b; Figure S11 in Additional file ) [,]. For example, the CSS contig 1750512_5AL harboring the Q gene harbors only one intronic SNP at position 1249. Regions harboring genes known to be associated with local adaptation ( Ppd, Vrn, Rht) also showed selection scan test statistic scores approaching the extremes (Table S20 in Additional file ). To further validate that our selection scans detect genomic regions associated with candidate loci controlling agronomic traits targeted by humans during the development of improved varieties, we compared selective sweeps with the marker-trait associations identified in mapping studies. We found 474 previously published associations that fell within the target regions (Table S21 in Additional file ). These markers showed association with major domestication and agronomic traits including spike length, rachis fragility, and compactness [,], heading date and flowering time [,], grain shape and yield characteristics [,], and nitrogen use efficiency []. The selective sweep regions also included markers associated with resistance to stripe rust [], bacterial leaf streak, and spot blotch [,]. The regions detected in the PHS scan overlapped with 459 marker-trait associations. Far fewer previously associated markers were located within the outliers of the F ST, XP-CLR, and diversity scans.

Non-synonymous, that is, likely functional, variants were significantly enriched ( χ2 test, P value <5.4 × 10⁻¹¹) in the extreme tails of the selection scans compared to synonymous variants. Regions identified by multiple selection scans were not common (Table S22 in Additional file ), and the regions where all three scans overlapped contained no genes with annotations. Among the overlapping regions detected using two different methods, one of the most common classes of genes were disease resistance genes (Table S10 in Additional file ). Consistently, the NB-ARC and LRR encoding domains of genes involved in disease resistance pathways [] were significantly over-represented in the extreme tail of the PHS scan, ( χ2 test; FDR <10⁻⁵ and <10⁻², respectively) (Table S23 in Additional file ) suggesting that some targets of selection can be associated with selection for disease resistance. To confirm this hypothesis, we performed GWAS of resistance to leaf, stem, and stripe rust (Figure b, Materials and methods) and tested for enrichment of marker-trait associations in the extreme tail of the selection scan. In the upper 2.5% tail of the PHS scan, we found three-, four-, and five-fold enrichment of SNPs ( χ2 test, P <6.5 × 10⁻⁶) associated with resistance to stem, leaf, and stripe rust, respectively, without the concomitant enrichment in the XP-CLR scan. Since the PHS test preferentially detects on-going selection events that have not reached fixation in a population [], our results suggest that multiple disease resistance genes undergo selection across wheat populations, which is likely associated with the spatial and temporal variation in pathogen populations, and consistent with the 'arms race' hypothesis [].

Among other candidates of selection is the WRKY transcription factor (Ta1dsLoc014113), identified by both selection scans and located on the short arm of chromosome 1D; its expression was shown to be associated with resistance to a fungal pathogen causing powdery mildew in wheat []. Glutathione-S-transferase encoding genes (Ta7dsLoc015003, Ta7bsLoc009692), that play an important role in drought response by reducing the toxicity of reactive oxygen species in wheat and other plants [,], were identified by the PHS and F ST scans on chromosomes 7B and 7D (Table S10 in Additional file ).

The development of the chromosome-specific wheat genome assemblies and population-scale whole exome re-sequencing data now provides the unique opportunity to investigate the impact of selection on duplicated copies of homoeologous genes. We have used the ordered sets of homoeologous genes to establish the syntenic relationships between the selection targets on different genomes. Inter-genomic comparison of selection targets associated with transition from landraces to cultivars revealed that only a single syntenic region shared the signature of selection between the B and D genomes (Table S24 in Additional file ). This region on the long arm of chromosome 1 shared two annotated genes encoding cellulose synthase (Ta1dlLoc001027, Ta1blLoc007155, Ta1blLoc025394), which plays a central role in cellulose biosynthesis, and trehalase (TH) (Ta1dlLoc015131, Ta1blLoc007983). The latter gene contains a PFAM domain that was significantly enriched in the tail of the XP-CLR scan (Table S23 in Additional file ), and along with trehalose phosphatase (TP), detected by both the XP-CLR and PHS scans (Table S10 in Additional file ), is involved in trehalose metabolism. The overexpression of TH was shown to increase drought stress tolerance in Arabidopsis [] and trehalose accumulation in transgenic rice expressing TP was associated with an increased tolerance to drought, salt, and cold stresses [].

The regions subjected to recent selection detected by the PHS scan showed a more substantial overlap in pairwise comparisons between the genomes than the regions detected using the XP-CLR and F ST scans (Table S24 in Additional file ) suggesting that during wheat improvement selection most likely operated on standing variation. However, the presumed adaptive variants with high-PHS in the overlapping homoeologous regions under selection are rarely found in the same line (Figure S19 in Additional file ). In the vast majority (75%) of overlapping homoeologous regions, no wheat lines in the panel possessed the high-PHS variants simultaneously in two genomes. Although there are known examples where allelic variation at the three wheat homoeologs affect the corresponding traits [-], our findings indicate that multiple parallel changes across homoeologous regions have rarely been favored by selection. Rather, it is likely that any favored variant at any one of the homoeologous regions may be sufficient to provide a fitness benefit, thereby expanding the target size of advantageous mutations. The rarity of individual genotypes with positively selected alleles in different homoeologous regions possibly indicate that additional advantageous mutations in other homoeologs either: (1) do not provide fitness benefit; or (2) are absent in a population subjected to improvement selection.

Gene expression studies have demonstrated that, in wheat, the homoeolog-specific transcriptional dominance affects up to 19% of genes [] with different homoeologs being preponderant in different groups of functionally related genes and showing the tissue- or development-specific patterns of expression []. These data are consistent with partitioning of gene function among the duplicated homoelogous genes [,], which could also affect the distribution of selective sweeps among the wheat genomes. One possibility is that selection acts on the preferentially expressed homoeolog. For example, most of the natural variation impacting the vernalization requirement in wheat is located in the VRN-A1 homoeolog [], which is the homoeolog expressed at the highest level [].

Conclusions

The sequence-based diversity map reported here is an important step towards the detailed characterization of DNA sequence polymorphism in the complex allopolyploid genome. The recently developed wheat genomic reference [] allowed us to catalogue common SNPs and small-scale indel polymorphisms from the low-copy fraction of the genome, describe their patterns of chromosomal distribution and inter-variant association, and identify variants that may have an impact on gene function based on the available annotation. A developed haplotype map will be a valuable tool for imputing genotypes and transferring sequence-level variation data across multiple gene mapping projects, thereby increasing the power and precision of trait mapping in GWAS and helping to understand better the basis of complex phenotypic traits.

Our data helped us gain insights into historic selective events and identify candidate selection targets associated with regions harboring genes controlling important agronomic traits or involved in response to biotic and abiotic stress stimuli. Our results suggest that directional selection in allopolyploids rarely acted on multiple parallel advantageous mutations across homoeologous regions. A favored variant at any one of the homoeologous regions appears to provide sufficient fitness benefit. By broadening the set of targets for selection, allopolyploidy may have played a critical role in the evolution of adaptation in wheat and contributed to wheat's success as a globally grown crop. Duplicated homoeologs may increase the likelihood of recovering beneficial alleles by expanding the advantageous mutation target size, and/or capturing allelic diversity present in different homoeologous genomes.

Materials and methods

Selection of wheat accessions

Capture assay design

A wheat exome capture (WEC) assay targeted the 107 Mb of non-redundant low-copy regions in the wheat genome. The capture probes were designed using the low-copy number genome assembly (LGC) of the wheat cultivar Chinese Spring []. The LGC had chloroplast, mitochondria, and transposon sequences removed and also contained homoeologous copies of genes collapsed into a single contiguous sequence []. The LGC was 3.8 Gb in size. We adopted two strategies to reduce the size over which we could design probes and target exonic regions. First, we used the BLASTN program (e-value <1e⁻¹⁰) to identify LGC contiguous sequences that were similar to Brachypodium exon sequences. Second, we used the same LGC sequence library (BLASTN e-value <1e⁻²⁰) to identify LGC contiguous sequences that matched a set of non-redundant wheat cDNA and EST sequences [] and transcriptome assemblies generated by 454 sequencing of nine diverse wheat cultivars []. Finally, to remove sequence duplications from the contiguous sequences set, we compared the set against itself using the BLASTN program. Similar sequences were identified (95% identity over 100 bp) and the longest contiguous sequence of a matching pair was retained. This process resulted in a design space of 110 Mb that was used to design a final probe set covering 107 Mb currently available from Roche NimbleGen []. The WEC assay was designed by the wheat-barley exome capture consortium that also designed an assay for the enrichment of the barley exome [].

DNA extraction and sequence capture

Each accession before DNA isolation was self-pollinated for two generations in the greenhouse. DNA was extracted using the DNeasy Plant Maxi Kit (Qiagen, USA) from a single 3-week-old seedling. One μg of DNA was fragmented with the Covaris S220 to obtain an average fragment length of 300 bp. The NEBNext DNA Library Prep Kit (NEB) for Illumina and Illumina TruSeq (TS) indexed (barcoded) adapters were then used for sample library preparation according to NEB protocol with the following exceptions. The PCR Enrichment step from the NEB was replaced with the Ligation Mediated PCR (LM-PCR). The TS-PCR Oligo1 (5′-AATGATACGGCGACCACCGAGA-3′) and TS-PCR Oligo2 (5′- CAAGCAGAAGACGGCATACGAG-3′) were used in the LM-PCR. The LM-PCR products were purified with the QIAquick PCR Purification Kit (Qiagen) followed by the size selection using Agencourt AMPure XP Beads (Beckman Coulter). The libraries were tested on 2100 Bioanalyzer (Agilent Technologies) and the NanoDrop Spectrophotometer (Thermo Scientific). Only samples with average fragment lengths of 200 to 400 bp, A 260/280 ratio 1.7 to 2.0, and LM-PCR yield >500 ng were pooled and used for sequence capture. Several levels of DNA sample pooling have been used in our study (Table S1 in Additional file ). One μg of non-pooled or pooled DNA (for example, 1 μg, 500 ng, 333 ng, 250 ng, or 125 ng of each component DNA sample library were used for 0, 2×, 3×, 4×, or 8× pools, respectively) was used in each of the sequence capture hybridizations. The sequence capture was performed as previously described [].

Genotyping by sequencing (GBS)

Complexity-reduced sequencing was performed according to the previously described protocol [] modified from Poland et al. []. A pooled library was sequenced on one lane of Illumina HiSeq 2000 (2 × 100 bp paired-end). All subsequent analyses were carried out using the same approaches for exome capture and GBS datasets.

Selection of alignment parameters for polyploid genome

Alignment to the CSS assemblies

Raw paired-end Illumina reads were quality filtered using the default setting of NGS QC toolkit v2.3 [], retaining reads if ≥70% of the bases had a quality score ≥20. Only paired-end reads were mapped to the CSS assemblies using Bowtie1 v.0.12.7 [] and Bowtie2 v.2.0.0 []. We applied the three-step iterative mapping strategy using Bowtie to perform ungapped read alignment (Figures S4 and S5 in Additional file ). Reads that do not align using more stringent criteria were reused for subsequent rounds of alignment with lower stringency. To find insertion/deletion (indel) polymorphisms, we performed gapped alignment using Bowtie2 v. 2.0.0 [] with the following parameters: -N 1, -L 75, -D 20, -R 3, --no-mixed, --end-to-end.

Efficiency of homoeologous target capture

The WEC assay included probes covering 107 Mb of non-redundant target space in the wheat genome. To assess the ability of the capture assay to enrich for targets from the three homoeologous wheat genomes, we compared the WEC against the CSS assemblies and retained only those that had three best BLASTN hits (e-value <1 × 10⁻¹⁰); one hit per genome. There were 47,739 homoeologous gene sets that fit these criteria. The log ₂ ratio of average coverage depth for each of these gene sets in pair-wise genome comparisons (A vs. B, B vs. D, and A vs. D) was distributed normally with the mean centered at 0, suggesting that homoeologous targets are captured with equal efficiency (Figure c).

To assess in silico the total size of the regions targeted by the WEC assay in the wheat genome we used the BLAT program [] to align the WEC design space against the CSS contigs (alignment length >100 bp, similarity >80-99%) (Figure S2 in Additional file ).

Capture efficiency

Analysis of alignments generated using the bowtie [] and BWA [] aligners showed that 95% of the 107 Mb WEC design space was covered at 30× depth and 99% of the design space was covered by at least one read. Approximately 78% of the annotated high confidence exons in the CSS [] were covered by at least one read. The average depth of read coverage for annotated exons in each accession was 8.1, 8.6, and 8.4 for the A, B, and D genomes, respectively (Table S4 in Additional file ), further suggesting no bias in the efficiency of homoeologous genome capture.

Variant calling and filtering

Each accession's BAM file was sorted and indexed using Samtools version 0.1.18 [] for variant calling. Next, GATK version 2.2-8 [] was used to realign reads around indels. The program Picard v. 1.62 [] was used to remove duplicate reads in the realigned BAM files. Finally, base quality recalibration was performed using the GATK program []. We identified 53.8 million raw variants using the Unified Genotyper from GATK following the GATK instructions for exome capture datasets. Subsequent recalibration of variant quality scores was performed using the default parameters of the VQSR tool in GATK []. As 'true variant' calls VQSR utilizes genotypes obtained with the wheat 90 K SNP assay []. All variants were filtered further to remove sites that had <46 (<75%) accessions genotyped or >2 alleles. We retained sites that had no more than one accession with a heterozygous call, or sites where the heterozygous call was due to a single read of the secondary allele (possibly sequencing error).

Due to the low-coverage depth obtained in the GBS dataset (on average of 1.04 read), genotype calling was performed only for alleles that were present in at least two different accessions in the population. The high level of genotype calling concordance between the WEC and GBS datasets was indicative of the high accuracy of the applied GBS genotype calling approach (see details in section ' Error rate estimation ').

Error rate estimation

We assessed the accuracy of genotype calling using four datasets. First, we compared variant calls generated for the cultivar Chinese Spring in our sequence capture experiment with the CSS assemblies that have also been generated using the same wheat cultivar. Of the filtered sites we found that 1,505,400 SNPs had genotype calls for cultivar Chinese Spring and 16,736 of these genotype calls were different from the reference sequence for an error rate of 1.1%. Similarly we estimated an error rate of 1.5% for indels where we found 1,576 indels that disagreed with the 106,438 indels that had a genotype call for Chinese Spring. The rare variant (MAF ≤1.6%) calling accuracy assessed in a set of 12,426 singletons unique to cultivar Chinese Spring depended on the read coverage depth (Figure S7 in Additional file ). By applying a ≥10 read coverage cutoff we obtained a singleton error rate of 4.6% for SNPs and 3.4% for indels.

Second, to exclude the possibility that the assessed error rates are not impacted by the quality of CSS reference assembly, we compared our variant calls generated for Chinese Spring with the published Chinese Spring BAC sequences generated using the Sanger approach []. We found 260 discrepancies out of 85,973 SNPs for an error rate of 0.3%.

Third, we estimated the concordance of our genotype calls with the genotype calls generated using the 90 K iSelect genotyping assay []. We selected only those SNP assays from the 90 K assay that were polymorphic in our population and whose flanking sequences could be unambiguously mapped to a single location in the CSS assemblies. Out of 27,147 homozygous genotypes called for these SNPs in our alignments, 520 were different from the SNP assay dataset, for an overall concordance rate of 98.1%.

Fourth, we have assessed the concordance of genotype calls by comparing the WEC and GBS datasets including the shared 10,028 SNP and 838 indel sites (Figure d). The level of genotype calling concordance between these two datasets was 97.2% for SNPs and 95.4% for indels.

SNP annotation

The impact of SNPs on coding sequences was assessed with the SNPEffect program [] using the 124,201 high confidence gene models predicted in the CSS contigs. We found that the proportion of SNPs identified in the intergenic regions is higher than that assigned to the coding regions. This pattern of SNP distribution is likely defined by several factors. First, WEC is designed using the low-copy fraction of wheat genome that overlapped with the sequences of wheat transcripts. These sequences do not always overlap with the high-confidence gene annotations of the CSS contigs used for annotating SNPs. It is possible that with better annotation of the wheat genome, new genes corresponding to these low-copy genomic regions could be found in WEC design. Second, the WEC co-captures regions outside the target region including introns and non-coding DNA. These regions are overall more genetically diverse and can also contribute to SNPs in the intergenic space.

Functional annotation and enrichment were determined using BLAST2GO software []. The ancestral state at SNP sites was inferred by comparing the SNP-flanking sequences among the A, B, and D genomes and with sequences of diploid wheat ancestors T. urartu [] and Ae. tauschii []. We were able to infer the ancestral allelic states for 754,080 of the 1.57 million SNPs.

Genetic diversity analyses

Due to variation in the depth of read coverage across the genome genotype calling rate for rare variants present only once in the population can vary from region to region. To reduce the effect of false negative rare variant calls on the local estimates of diversity, all analyses were performed using the variants with MAF >1.6%.

Recently, it was demonstrated that the estimates of commonly used diversity statistics ( F ST, Tajima's D, and π) obtained using restriction enzyme-based sequencing can deviate from true values []. Ascertainment bias in GBS data can result from the mutations at the restriction enzyme cut sites, sensitivity of the PstI enzyme to DNA methylation, or a high proportion of missing genotypes in low-coverage re-sequencing datasets. Allele frequency estimates obtained for the same sites in the GBS and WEC datasets showed good correlation (Figure S20 in Additional file ). However, we observe a slight decrease in correlation values with the increase of the proportion of missing data in the GBS dataset. To reduce the potential effect of missing data on our analyses this dataset was excluded from the estimates of diversity statistics.

The order of the CSS contigs along the wheat chromosomes was established using the combination of 'genome zipper' [] and a high-density genetic map developed by the population sequencing (PopSeq) of 90 recombinant inbred lines []. The approximate physical positions of ordered CSS contigs on the chromosomes were inferred using a framework created by cross-linking high-density SNP wheat genetic map [], wheat deletion bin map [], and barley genome sequence []. Comparison of the inferred physical positions with the positions of CSS contigs on the published sequence of the chromosome 3B pseudomolecule [] showed a high level of correlation (Pearson correlation r ² = 0.97). The ordered CSS contigs were used to perform sliding window analyses. π and Tajima's D estimates for 2 Mb sliding windows (1 Mb step) with >10 kb of CSS contigs covered by reads were calculated using LDHat (v. 2.2) []. Outliers of π and Tajima's D were defined by the 2.5 percentiles of the test statistic distribution. F ST was calculated by SNP using the R package adegenet (cran.us.r-project.org) by contrasting cultivars and landraces. A sliding window analysis (2 Mb with a 1 Mb step) was used to generate average estimates of F ST. Outliers were calculated as the 97.5 percentile of all window values within the genome. Adjacent outlier windows were merged. The thresholds for the cultivar-landrace comparison were 0.158 for the A genome, 0.108 for the B genome, and 0.0925 for the D genome.

For estimating recombination rate we selected 10,517 CSS contigs that contained at least 20 non-singleton SNPs. The historical recombination parameter (ρ = 4Ne) was estimated using a composite likelihood approach implemented in the Maxhap program []. Regions of the wheat genome showing elevated levels of historic recombination were defined as those falling above 97.5 percentile of the genome-wide ρ distribution. The critical values of ρ for each genome were 7.5 × 10⁻³, 8.5 × 10⁻³/kb, and 9.7 × 10⁻³/kb for the A, B, and D genomes, respectively.

Pair-wise estimates of LD were obtained for SNPs with a MAF ≥0.05 by measuring r ² as previously described []. The average length of pair-wise shared haplotypes in the population was calculated around each SNP according to the previously described procedure []. The genetic map distances were taken from the wheat SNP map developed using the 90 K SNP assay [].

Distribution of PFAM domains in the CSS contigs

The nucleotide sequences of the CSS gene models for T. aestivum were mapped against the wheat survey sequence using GMAP []. Output was filtered for hits with ≥98% identity and ≥80% coverage. PFAM domains were identified in the amino acid sequences of the CSS gene models using PFAM hidden Markov models (release 26) from the Sanger Institute [] and HMMER (version 3.0b3) from the Howard Hughes Medical Institute Janelia Farm Research Campus []. The PFAM searches were performed using the default parameters and the results were filtered for hits with ≥90% coverage.

Phylogenetic tree construction

Distribution of Ae. tauschii (genotype AL8/78) alleles across the D genome chromosomes

To determine which of the D genome alleles corresponds to the AL8/78 genotype, 100 bp flanking sequences of each D genome SNP were extracted and compared against the genomic sequence of Ae. tauschii AL8/78 genotype []. In total, we could determine AL8/78 alleles for 105,294 SNPs. The distribution of AL8/78 alleles along the wheat chromosomes was estimated in the population of the 26 landraces. The average frequency of the AL8/78 allele in a sliding window of 3 Mbp was plotted along the wheat chromosomes.

Detection of selective sweeps

The PHS statistic was calculated as described by Toomajian et al. []. Utilizing this statistic we have higher likelihood of detecting genomic regions harboring alleles present at intermediate frequencies in the population [] than those alleles that have reached high frequency. Thresholds of the PHS statistic were determined by taking the 97.5 percentile of the distribution of PHS values for each SNP within different allele frequency classes (in a window of 0.05). For detecting selective sweeps the wheat genome was split into 50-kb windows; each window was assigned a maximum PHS value for SNPs in the window. Neighboring windows located within 1 Mb that contained outlier SNPs were merged. Annotated genes harboring the outlier SNPs were used for PFAM domain enrichment analyses. Enrichment analysis of GWAS SNPs included sites with a minor allele frequency >3% and a significance level P <10⁻³.

The selection scan using the XP-CLR approach is robust to assumptions regarding recombination rates and demography [] and compares the allele frequency differentiation and the extent of linked variation between two populations (cultivars vs. landraces) to detect regions where change in frequency occurred too quickly to be caused by random drift. The XP-CLR scan was run with the grid size of 50 kb, the window size of 1 cM, and the maximum number of SNPs fixed at 500. The critical values for putative selection targets were estimated based on the 97.5 percentile of the test statistic distribution for each wheat genome. Since our population includes both spring and winter wheat lines, one of the concerns was that population differentiation between these growth habit groups can be mistaken for the signals of selection. However, we believe that at the genome-wide level, inclusion of spring and winter wheat lines should not have a large effect on the detected signals of selection because the level of genetic differentiation between the spring and winter wheat in our population was shown to be very low (mean genome-wide F ST = 0.03). This low F ST was observed previously [] and attributed to the common practice to use lines from both growth habit groups in the same breeding programs, as well as to the heterogeneity of the genetic basis of flowering time regulation in wheat where the same phenotypic outcome can be obtained by mutations at several independent genetic loci. In addition, we have identified regions of the wheat genome showing extreme genetic differentiation ( F ST) or XP-CLR test statistics between spring and winter wheat. Out of 372 XP-CLR outliers in the landrace-cultivar comparison only nine (2.4%) partially overlapped with the outliers in the spring-winter wheat comparison (shown in Table S18 in Additional file ). Out of 168 outliers in the F ST scan of cultivars and landraces only six (3.6%) overlapped with the F ST outliers in the spring-winter wheat comparison. No overlap of genetically differentiated genomic regions between spring and winter wheat with the outliers of the genetic diversity scan was found. None of the genes located in the genomic regions that showed the signature of selection in at least two scans reported in Table S10 (Additional file ) fall within the genomic regions differentiated between spring and winter wheat. Overall, these analyses suggest that the genetic differentiation between spring and winter wheat in our sample should not have significant impact on the results of our selection scans.

To test the significance of the observed overlap between the selective sweeps located on the homoeologous chromosomes, we have randomly permuted 10,000 times a genome-wide set of 50-kb windows. The proportions of windows that overlapped among each pairwise comparison of the wheat genomes was ranked and compared to observed data to calculate empirical P value.

Comparison of selective sweeps with previously characterized marker-trait associations

To check if the putative selective sweeps identified in our scans harbored loci associated with agronomic traits, we analyzed published marker-trait associations detected using the DArT markers [], and 9 K and 90 K iSelect SNP assays [,]. The DArT markers were mapped to the CSS contigs using the BLAT program (best hit with minimum alignment length >150 bp). Out of 57 DArT markers that could be mapped to the ordered CSS contigs, 37 markers fell into the regions detected in the PHS scan, and one fell into a region detected by both the PHS and F ST scans (Table S21 in Additional file ). Out of 555 SNP markers mapped to the ordered CSS contigs, 422 mapped to the regions identified in the PHS scan. We found five SNPs in the XP-CLR regions, three of which overlapped with the PHS regions. One SNP marker was located within the region showing reduced level of diversity. There were 30 SNP-trait associations that fell into high F ST regions, of which 18 overlapped with the high-PHS regions.

Imputation

A similar imputation strategy was used for imputing genotypes in the panel of 678 wheat lines used for GWAS of disease resistance.

Genome-wide association study (GWAS)

Plant materials

A total of 838 bread wheat accessions from the Arthur Watkins Collection were obtained from the Australian Winter Cereals Collection, Tamworth. This collection includes a large number of phenotypically diverse wheat landraces collected from 32 countries in the 1920s to 1930s []. The 838 accessions were grown under field conditions and single plant selections made for 678 accessions on the basis of plant type, rust resistance, and maturity. Despite some maturity differences, all 678 land race accessions flowered by the end of October. Seed from the purified 678 accessions can be obtained from the Plant Breeding Institute, Cobbitty, upon request (contact urmil.bansal@sydney.edu.au).

Phenotyping

Following seed bulk up for the single plant selections, the purified 678 accessions were grown under field conditions during 2006 and 2007 at the Lansdowne site of the University of Sydney. Each accession was grown as a single 1 m row. Field trials were artificially inoculated with Puccinia striiformis f. sp. tritici (Pst) pathotype 134 E16A+; P. triticina (Pt) pathotypes 104-1,(2),3,(6),(7),11 + Lr37; 104-1,(2),3,(6),(7),11,13; 104-1,(2),3,(6),(7),9,11; 76-3,5,9,10 + Lr37; 10-1,3,(7),9,10,11,12 and P. graminis f. sp. tritici (Pgt) pathotypes 98-1,2,3,5,6 and 34-1,2,7 + Sr38. The pathotype designations are provided according to McIntosh et al. []. These pathotypes carry partial virulence for the genes given in the parenthesis. Stripe rust, leaf rust, and stem rust responses were recorded using a 1-9 scale, where 1 was highly resistant and 9 was highly susceptible []. For stripe rust disease, two records were taken within each year. The consistency of phenotypic evaluations across years was tested by calculating the Pearson correlation coefficient. For leaf, stripe, and stem rust phenotyping datasets the correlation coefficients were 0.75, 0.71-0.87, and 0.57, respectively.

Genotyping

The 678 landrace accessions were genotyped using the Infinium iSelect 90 K SNP assay, the content of which is reported to have minimal ascertainment bias for the analysis of diverse wheat landraces []. Genotyping was performed on the iScan instrument according to the manufacturer's protocols (Illumina). SNP genotype calling was performed using GenomeStudio v2011.1 software (Illumina) and the genotype calling algorithm reported in Wang et al. []. Monomorphic markers and SNPs with more than 10% missing data (due to the presence of null alleles or poor genotype call rates) were removed. The genetic map developed by genotyping multiple mapping populations with the 90 K array was anchored to the GenomeZipper and PopSeq maps by comparing the sequences of 90 K array SNPs with the sequences of wheat chromosome assemblies [].

Association mapping

Mixed model variance component analysis of mean phenotypic values was performed using the R package GAPIT []. The information about the relationship among accessions in the population was provided as a kinship matrix (random effect). The effect of population structure was controlled by using the first three principal components from the principal component analysis (fixed effect). The P value <1 × 10⁻³ used to filter markers was selected based on the previous studies of agronomic traits in wheat and barley demonstrating that this threshold provides adequate accuracy for detecting marker-trait associations, as was validated independently in bi-parental mapping populations [,,]. In a follow-up analysis we successfully validated five GWAS regions by comparing with previously published studies, or by mapping a trait in a bi-parental mapping population in our study (shown on Figure b and in Table S14 in Additional file ).

Validation of GWAS signals in mapping populations

Plant materials

Mapping in the populations of recombinant inbred lines (RIL) was used to validate several marker-trait associations identified in the GWAS including two characterized and one previously uncharacterized disease resistance loci. The population segregating for leaf rust resistance gene Lr67 included 124 F ₃/F ₄ lines derived from a cross between Thatcher and RL6077 []. The population segregating for stripe rust resistance gene Yr51 was comprised of 89 F ₆ RILs derived from a single heterozygous F3 plant #5515 from a cross between Watkins' line PBI769 and Westonia []. The population segregating for a new stem rust resistance gene on chromosome 7A was comprised of 96 F ₆ RILs derived from a cross between Watkins' line PBI562 and Yitpi (henceforth, PBI562/Yitpi).

Identification of SNPs linked to rust genes

SNPs associated with rust resistance genes segregating in each mapping population were identified using bulk segregant analysis (BSA) [] and selective genotyping (SG).

For BSA, resistant and susceptible bulks were prepared by pooling equal amounts of genomic DNA from at least 20 plants for each phenotypic class. An artificial F1 sample was prepared by combining an equal amount of DNA from each of the two bulks. The bulked DNA samples, artificial F1, and parental lines were genotyped using a custom Infinium iSelect bead chip assay on the iScan instrument following the manufacturer's instructions (Illumina Ltd.). The Thatcher/RL6077 and PBI No. 769/Westonia were genotyped using the 9 K and 90 K iSelect genotyping arrays, respectively. The SNPs were assessed for putative linkage by comparing the normalized theta values for each sample as described in Hyten et al. []. Polymorphism was considered to be linked to a rust resistance gene when the normalized theta values for the resistant bulk and resistant parent, and susceptible bulk and susceptible parent were similar, and when the normalized theta value for the artificial F ₁ samples was about halfway between that of the other samples.

For SG, at least 15 resistant and 15 susceptible plants from PBI No. 562/Yitpi cross and its parents were genotyped using the 90 K iSelect bead chip assay. Polymorphism was considered to be linked to a rust resistant gene when the majority (>90%) of individuals within each phenotypic class were fixed for the expected parental allele.

Validation of GWAS-linked SNPs using bi-parental mapping crosses

The genomic location of SNPs associated with rust disease resistance in the GWAS analysis were compared with those linked to mapped rust resistance genes in each of the bi-parental mapping crosses. GWAS-associated SNPs were considered to co-locate with the mapped resistance genes when the linked SNPs in each study were located at (or very near) the same position in the 90 K consensus SNP genetic map []. In instances where the linked markers in either study were not present in the 90 K consensus SNP map, the GWAS-associated SNPs were considered to co-locate with the mapped resistance genes when the CSS contigs tagged by the SNPs [] co-located in the PopSeq map []. The co-location of SNPs associated with rust disease resistance in the GWAS analysis with mapped rust resistance genes in the bi-parental mapping populations is shown in Figure b and Table S14 in Additional file .

Data availability

Notes

Katherine W Jordan and Shichen Wang contributed equally to this work.

Declarations

Acknowledgments

This project was supported by the National Research Initiative Competitive Grants 2011-68002-30029 (Triticeae-CAP) (JD) and 2012-67013-19401 from the USDA National Institute of Food and Agriculture (EA), by the Bill and Melinda Gates Foundation grant (EA), by Genome Prairie, Genome Canada, Saskatchewan Ministry of Agriculture, Western Grains Research Foundation (CP and PH), by a BBSRC Career Development Fellowship BB/H022333/1 (AH) and Doctoral Training Grant (LG), by the Howard Hughes Medical Institute and the Gordon and Betty Moore Foundation (JD), and by Kansas Agricultural Experiment Station contribution no. 15-338-J. We thank Cyrille Saintenac for help at the initial phases of the project and Daniel Andresen for assistance with the computing resources of KSU Beocat cluster funded by NSF grant ACI-144054.

Additional files

Competing interests

JAJ recognizes a competing interest as an employee of Roche NimbleGen Inc.

Authors' contributions

Authors' information

International Wheat Genome Sequencing Consortium.

KJ and SW performed most of the data analyses with support from RM, CP, CT and TC. AH, NH, LG, CP, EA, and JJ developed and tested the wheat exome capture assay. LT and JD performed selection of wheat lines and participated in manuscript writing. MH, LG, and AH assessed the accuracy of genotype calling and participated in manuscript writing. AA, YL, CS, KW, and AS performed sequence capture and DNA sequencing and contributed to preparation of manuscript. DW, KLF, and MH performed genotyping of GWAS panel. UKB, HSB, and MJH performed GWAS. EA proposed the idea, coordinated data analyses, and wrote the first draft of the manuscript with the help from KJ. The gene annotations, CSS contig assemblies, and CSS contig order along the chromosomes are provided by the IWGSC (). All authors read and approved the final manuscript.

Authors' Affiliations

Department Plant Pathology, Kansas State University, Manhattan, USA

Integrated Genomics Facility, Kansas State University, Manhattan, USA

Institute of Integrative Biology, University of Liverpool, Liverpool, UK

Department Plant Sciences, University of Saskatchewan, Saskatoon, Canada

Department Environment and Primary Industries, Bundoora, Australia

National Research Council Canada, Saskatoon, Canada

Department Botany & Plant Sciences, University of California, Riverside, USA

Department Plant Sciences, University of California, Davis, USA

Howard Hughes Medical Institute, Chevy Chase, USA

Department Plant Sciences & Plant Pathology, Montana State University, Bozeman, USA

Plant Breeding Institute-Cobbitty, The University of Sydney, PMB4011, Narellan, Australia

Download PDF

Papers, Zotero, Reference Manager, RefWorks (.RIS)

See updates

# data[].attributes.attachment_link.title
A haplotype map of allohexaploid wheat reveals distinct patterns of selection on homoeologous genomes

# data[].attributes.attachment_link.description
Bread wheat is an allopolyploid species with a large, highly repetitive genome. To investigate the impact of selection on variants distributed among homoeologous wheat genomes and to build a foundation for understanding genotype-phenotype relationships, we performed population-scale re-sequencing of a diverse panel of wheat lines.

# data[].attributes.body
A map of rice genome variation reveals
the origin of cultivated rice

Please put a special emphasize on selective sweeps

Presenters: Kiera Thurman, Fionn Cahalane, Edward Gahan

# data[].attributes.attachment_link.readable_text
Abstract

Crop domestications are long-term selection experiments that have greatly advanced human civilization. The domestication of cultivated rice ( Oryza sativa L.) ranks as one of the most important developments in history. However, its origins and domestication processes are controversial and have long been debated. Here we generate genome sequences from 446 geographically diverse accessions of the wild rice species Oryza rufipogon, the immediate ancestral progenitor of cultivated rice, and from 1,083 cultivated indica and japonica varieties to construct a comprehensive map of rice genome variation. In the search for signatures of selection, we identify 55 selective sweeps that have occurred during domestication. In-depth analyses of the domestication sweeps and genome-wide patterns reveal that Oryza sativa japonica rice was first domesticated from a specific population of O. rufipogon around the middle area of the Pearl River in southern China, and that Oryza sativa indica rice was subsequently developed from crosses between japonica rice and local wild rice as the initial cultivars spread into South East and South Asia. The domestication-associated traits are analysed through high-resolution genetic mapping. This study provides an important resource for rice breeding and an effective genomics approach for crop domestication research.

Main

Cultivated rice ( Oryza sativa L.), which is grown worldwide and is one of the most important cereals for human nutrition, is considered to have been domesticated from wild rice ( Oryza rufipogon) thousands of years ago1,2,3,4. The differences between O. sativa and O. rufipogon are reflected in a wide range of morphological and physiological traits5,6,7,8,9. Despite the fact that rice is a major cereal and a model system for plant biology, the evolutionary origins and domestication processes of cultivated rice have long been debated. The puzzles about rice domestication include: (1) where the geographic origin of cultivated rice was, (2) which types of O. rufipogon served as its direct wild progenitor, and (3) whether the two subspecies of cultivated rice, indica and japonica, are derived from a single or multiple domestications.

A wide range of genetic and archaeological studies have been carried out to examine the phylogenetic relationships of rice, and investigate the demographic history of rice domestication10,11,12,13,14,15,16,17,18,19. Molecular phylogenetic analyses indicated that indica and japonica originated independently3,10,20. However, the well-characterized domestication genes in rice were found to be fixed in both subspecies with the same alleles, thus supporting a single domestication origin6,7,8,9,16. Recently, a demographic analysis of single-nucleotide polymorphisms (SNPs) detected from 630 gene fragments suggested a single domestication origin of rice17. Meanwhile, population genetics analyses of genome-wide data of cultivated and wild rice have tended to suggest that indica and japonica genomes generally appear to be of independent origin18,19, but many genomic segments bearing domestication alleles may have originated only once18. Despite these advances, wider sampling with population-scale whole-genome sequencing is needed to shed greater light on the evolutionary history of rice domestication. An in-depth investigation of the haplotype structure near the domestication sites will be critical for evaluating the direction of introgression. The specific ancestral population and the subsequent demographic event are yet to be identified.

Moreover, a comprehensive map of rice genome variation will facilitate genetic mapping of complex traits in rice. We recently collected diverse rice cultivars for sequencing, and carried out genome-wide association studies (GWAS) for many agronomic traits in cultivated rice21,22. Here we sequenced and analysed the genomes of 446 O. rufipogon accessions to investigate the phylogenetic relationships between cultivated and wild rice and identify the signatures of selection in rice domestication. This research also provides a robust foundation that will enable rice breeders to effectively exploit diverse genetic resources for rice improvement.

Analysis of wild rice populations

The strategy of this study is briefly described in Supplementary Fig. 1. The genus Oryza consists of 23 species23,24, and the wild rice O. rufipogon is believed to be the immediate progenitor of the cultivated rice O. sativa ( Supplementary Fig. 2 and Supplementary Table 1). From large collections of wild rice germplasm maintained in China and Japan, we selected 446 diverse O. rufipogon accessions including both perennial and annual (also called as Oryza nivara) forms from Asia and Oceania, spanning the native geographic range of the species ( Supplementary Tables 2 and 3). We sequenced these accessions with twofold genome coverage. After aligning the reads against the rice reference genome sequence, we identified a total of 5,037,497 non-singleton SNPs. Based on the SNP data, the sequence diversity ( π) of O. rufipogon was estimated at ∼ 0.003, which is higher than that in O. sativa-the sequence diversity is 0.0024 for O. sativa, and 0.0016 and 0.0006 for indica and japonica, respectively21. Approximately 82% of SNPs (minor allele frequency > 0.05) segregating in O. rufipogon also segregate in O. sativa ( Supplementary Table 4). This observation is consistent with previous suggestions that part of the genetic diversity in the progenitor would be lost because only a limited number of individuals were used during domestication11,25.

We investigated the population structure of the O. rufipogon accessions. On the basis of the neighbour-joining tree, as well as principal-component analysis (PCA), we classified the O. rufipogon species into three types, simply designated as Or-I, Or-II and Or-III in this study ( Fig. 1a and Supplementary Fig. 3). We found that the rice population structure strongly correlated with geographic distribution ( r ² = 0.2 between the first principal component and the longitude, and r ² = 0.3 between the second principal component and the latitude)26. Interestingly, the O. rufipogon accessions sampled from southern China mostly belong to the Or-III type ( Fig. 1b and Supplementary Fig. 4). The level of population differentiation, F ST, was estimated at 0.18 among the groups of O. rufipogon, which is much lower than that of O. sativa (∼ 0.55 on average21). The differentiation was not evenly distributed across the rice genome ( Supplementary Fig. 5). We scanned the whole genome for highly differentiated loci and found 68 loci with F ST > 0.3, which covered ∼ 3% of the complete rice genome ( Supplementary Table 5). Of these loci, we found several known genes or quantitative trait loci (QTLs) that included DPL2 (hybrid incompatibility27) ( Fig. 1c), OsSOC1 and Ghd7 (flowering time28,29,30) and qCTS12 (cold tolerance31), all of which have been reported to be closely related to indica- japonica differentiation in rice. Hence, the highly differentiated loci can provide important clues for searching the genes involved in local adaptation.

Figure 1: Genetic structure and association analysis in the wild rice population.

a, Neighbour-joining tree of 446 O. rufipogon accessions, which was calculated from ∼ 5 million SNPs, identifies the three groups of Or-I (red), Or-II (grey) and Or-III (blue). b, Geographic origins of wild rice accessions. c, The level of genetic differentiation ( F ST) in O. rufipogon population around the DPL2 gene that underlies indica- japonica hybrid incompatibility in rice. d, Regional Manhattan plots of GWAS for tiller angle in O. rufipogon population identify a known gene, PROG1, using a compressed mixed linear model. The genome-wide significance threshold (1 × 10 −6) and the position of the peak SNP are indicated by a horizontal dash-dot line and a vertical red line, respectively.

Because O. rufipogon is an out-crossing species, it was expected to have a relatively high decay rate of linkage disequilibrium. We found that the decay rate of linkage disequilibrium in O. rufipogon, expressed as r ², dropped to half of its maximum value at ∼ 20 kilobases (kb) on average ( Supplementary Fig. 6), which is much more rapid than that in O. sativa (∼ 123 kb and ∼ 167 kb in indica and japonica, respectively)21. To perform GWAS in O. rufipogon, we used the k-nearest neighbour algorithm for data imputation ( Supplementary Table 6), and phenotyped the O. rufipogon population for two traits, leaf sheath colour and tiller angle. The strongest associations for sheath colour and tiller angle were found to be just around the known loci OsC1, for colouration32, and PROG1, for prostrate growth7,8 ( Fig. 1d and Supplementary Figs 7-9). Through computational simulations we predicted that the mapping resolution of GWAS in O. rufipogon was approximately three times greater than that in O. sativa on average ( Supplementary Figs 10-12). Hence, the wild rice population and accompanying comprehensive sequence resource should be of great utility for directly dissecting agronomic traits in rice.

Phylogenetic relationships of rice

We used whole-genome sequencing data for a large panel of accessions containing 446 O. rufipogon accessions and 1,083 O. sativa varieties to explore their phylogenetic relationships. The 1,083 diverse O. sativa varieties were collected throughout the world and sequenced with one-fold genome coverage ( Supplementary Table 7). Of these, 950 genomes had previously been reported21,22, and the remaining 133 genomes, including many representative varieties (for example, aromatic rice), were sequenced and first reported in this study. A total of 7,970,359 non-singleton SNPs were identified from the 1,529 genome sequences. To determine the ancestral states of the SNPs, the close relatives of O. rufipogon and O. sativa were also sequenced with a total of approximately 50-fold genome coverage ( Supplementary Table 8). Using genome sequences of these outgroups, we were able to identify the ancestral states of 6,119,311 SNPs out of the approximately 8 million SNPs (77%).

We used the genotype data set of the ∼ 8 million SNP sites from 1,529 genomes to infer genome-wide relationships. Both the phylogenetic tree and the PCA plots indicate that O. sativa indica and japonica are descended from Or-I and Or-III, respectively ( Fig. 2a and Supplementary Fig. 13). We then investigated the detailed relationship between indica and Or-I and that between japonica and Or-III separately ( Supplementary Figs 14-16). The level of genetic differentiation between indica and Or-I was modest ( F ST = 0.17), and indica contains approximately 75% of the genetic diversity in Or-I ( Fig. 2b). A small number of indica and Or-I accessions seemed to be intermediate between cultivated and wild rice ( Fig. 2a). In contrast, the japonica groups ( temperate japonica, tropical japonica and aromatic) were all clustered together, and had descended from wild rice in southern China (sub-clade Or-IIIa). There was an obvious genetic distinction between japonica and Or-IIIa ( Fig. 2a), which had a relatively high level of population differentiation ( F ST = 0.36). We found that only approximately 33% of the genetic diversity of Or-III persisted in japonica ( Fig. 2b). The strong genetic bottleneck indicates that a small effective population from Or-III was used for domesticating japonica cultivars.

Figure 2: Genome-wide relationship between cultivated rice and its wild progenitor.

a, Phylogenetic tree of the full population (446 O. rufipogon accessions and 1,083 O. sativa varieties) calculated from ∼ 8 million SNPs in O. rufipogon and O. sativa. The double-layer rings indicate O. rufipogon (outer ring: Or-I, Or-II and Or-III are coloured in red, grey and blue, respectively) and O. sativa (inner ring: indica and japonica subspecies are in pink and sky blue, respectively). b, Illustration of genetic diversity and population differentiation in O. rufipogon and O. sativa. The size of the circles represents the level of genetic diversity ( π) of the groups, and the F ST values between the groups are indicated. ind, indica; jap, japonica. c, The spectrum of allele frequencies at the causal polymorphisms of Ghd7, DPL2 and GS3.

On the basis of our SNP data, most sequence variants observed between indica and japonica, which were estimated to have diverged hundreds of thousands of years ago20, already existed in the progenitor gene pools. However, the differentiation was enhanced during domestication, with F ST expanding from 0.18 in O. rufipogon to 0.55 in O. sativa ( Fig. 2b). The high level of genetic differentiation between indica and japonica has resulted from a combination of the modest differentiation present within their ancestor population and the recent domestication bottleneck ( Fig. 2c). For example, at the causal SNP site for GS3 (a major QTL for grain shape33,34), the null allele mainly existed in the wild rice populations in the Guangxi and Guangdong provinces of southern China. This allele rapidly extended to become the major allele (∼ 98%) in the japonica population, thus generating a distinct phenotypic difference between indica and japonica. We screened all SNPs that were highly differentiated in frequency between indica and japonica, and found a total of 213,188 indica- japonica-differentiated SNPs. With regard to the differentiated SNPs, the differences present within the O. rufipogon, the domestication bottleneck in japonica and in indica were estimated to contribute ∼ 50%, ∼ 40% and ∼ 10% of total alterations in allele frequency, respectively ( Supplementary Fig. 17).

We examined the ancestral states of indica- japonica-fixed SNPs, and found that 55% of SNPs of the ancestral alleles are fixed in indica and the other 45% are fixed in japonica. The d N/d S ratio (d N, number of non-synonymous substitutions per non-synonymous site; d S, number of synonymous substitutions per synonymous site) of the SNPs was calculated to be 0.34, which was almost equal to the average level of the total SNPs. In contrast, for the 9,595 SNPs that were fixed between O. rufipogon and O. sativa, the ancestral alleles of 93.3% SNPs are identical to O. rufipogon, indicating that O. rufipogon has retained more ancestral states than O. sativa, which also further supported our conclusion that O. rufipogon is likely to be the ancestral progenitor of O. sativa. The d N/d S ratio of the fixed SNPs between O. rufipogon and O. sativa is calculated to be 1.04, indicating positive selection during domestication ( P  < 0.001, chi-square test).

Screens and annotation of domestication loci

Selective signatures from domestication include a reduction in nucleotide diversity and altered allele frequency in the domestication loci19,25. We measured the ratio of the genetic diversity in wild rice to that in cultivated rice ( π w/ π c) across the rice genome, and determined the cutoff on the basis of permutation tests ( Supplementary Information section 2). We performed whole-genome screening in indica, japonica and the full population using the diversity ratios ( Fig. 3 and Supplementary Figs 18-20). In total we identified 60 loci in indica, 62 in japonica and 55 in the full population ( Supplementary Tables 9-11). We noticed that many loci with strong signals of selection were nearly identical in both indica and japonica where F ST between indica and japonica was extremely low, indicating that introduction of traits during domestication has in many cases involved introgression events. We noted that most well-characterized domestication genes, including Bh4 (hull colour9), PROG1 (tiller angle7,8), sh4 (seed shattering5,6), qSW5 (grain width35) and OsC1 (leaf sheath colour and apiculus colour32), were among the 55 loci detected in the full population ( Fig. 3). Another three well-characterized domestication genes, qSH1 for seed shattering36, Waxy for grain quality37 and Rc for pericarp colour38,39, which showed strong selection signals in the japonica panel, were not fully shared in the indica population.

Figure 3: Whole-genome screening and functional annotations of domestication sweeps.

a, Whole-genome screening of domestication sweeps in the full population of O. rufipogon and O. sativa. The values of π w/ π c are plotted against the position on each chromosome. The horizontal dashed line indicates the genome-wide threshold of selection signals ( π w/ π c > 3). b- d, A large-scale high-resolution mapping for fifteen domestication-related traits was performed in an O. rufipogon ×  O. sativa population. The domestication sweeps overlapped with characterized domestication-related QTLs are shown in dark red, and the loci with known causal genes are shown in red. Among them, three strong selective sweeps were found to be associated with grain width ( b), grain weight ( c) and exserted stigma ( d), respectively. In b- d, the likelihood of odds (LOD) values from the composite interval mapping method are plotted against position on the rice chromosomes. Grey horizontal dashed line indicates the threshold (LOD > 3.5).

We investigated the genetic patterns around the well-characterized domestication loci, which were quite different from those at the whole-genome scale. According to phylogenetic trees calculated from SNPs on these loci, most cultivars were clustered together, and the Or-III population from southern China tended to be the closest ancestral progenitor of all the cultivars ( Supplementary Fig. 21). These patterns were much clearer when expanded to the total regions of the 55 domestication-related loci ( Fig. 4a). The indica cultivars were generally close to Or-I across the genome, but were closer to Or-IIIa than to Or-I at the 55 selected loci ( Supplementary Fig. 22). We further calculated the genetic distances between cultivated rice and the wild rice populations from each geographic sampling region through analysis of the genomic regions around the 55 domestication-related loci ( Supplementary Fig. 23). Our genetic approach showed that the middle area of the Pearl River in Guangxi province, southern China, was probably the place of the first development of cultivated rice ( Fig. 4b-d), although an archaeological finding had identified the Lower Yangtze region in eastern China as one of the centres of rice cultivation14,15. These results suggest a model in which japonica was first domesticated from Or-III in southern China, and was subsequently crossed to local wild rice in South East Asia and South Asia, thus generating indica after many cross-differentiation-selection cycles ( Fig. 4e). Furthermore, we performed computational simulations to generate in silico data sets under various demographic scenarios ( Supplementary Fig. 24), using the forward simulator SFS_CODE40 ( Supplementary Information section 5). We calculated the genetic distances of indica and japonica with different clades of O. rufipogon for both the real and simulated data ( Supplementary Fig. 25), which showed that there were significant differences between the real data and the simulated ones in all the scenarios except for the model proposed here ( P  < 0.01, rank sum test), providing further support for the model that we propose.

Figure 4: Genetic and geographic origins of rice domestication.

a, Phylogenetic tree of 446 O. rufipogon accessions and 1,083 O. sativa varieties calculated from SNPs in the overall regions of the 55 major domestication sweeps. b, Geographic locations of 62 O. rufipogon accessions, whose phylogenetic positions during domestication are indicated. Colour index represents the average of the genetic distance of O. rufipogon accessions to all cultivated rice accessions. Two major rivers in southern China are labelled in grey in the map. c, The average distance of O. rufipogon accessions from different countries to all cultivars. The distance was estimated by simple matching distance of SNPs around the Bh4 locus or all SNPs within the 55 domestication sweeps. d, The average distance of O. rufipogon accessions from different provinces in southern China to all cultivars. e, Schematics of the origin of cultivated rice. The aus and aromatic rice are minor groups of rice accessions with small geographic distributions.

To investigate why these chromosomal loci were selected during rice domestication, QTL mapping of domestication-related traits was performed using a population that was developed from a cross between O. sativa indica Guangluai-4 and O. rufipogon (Or-IIIa accession W1943). Using a sequencing-based genotyping method41, we sequenced 271 lines in the population with 0.5× coverage for each, and constructed an ultra-dense genotype map. A total of 15 domestication-related traits were phenotyped, with 58 QTLs detected using a composite interval mapping method ( Supplementary Fig. 26 and Supplementary Tables 12, 13). For ten QTLs for which causal genes have been reported, almost all the peak signals were within 200 kb from the known genes, indicating that a high mapping resolution was achieved ( Supplementary Fig. 27).

Of the 58 QTLs detected, 32 QTLs were located within domestication sweeps ( Fig. 3b-d and Supplementary Figs 28-30). The overlap of the QTLs for domestication traits with domestication sweeps is very significant ( P < 0.001, chi-square test). Among them, the QTLs for traits like exserted stigma (controlling mating system) and grain size (controlling output yield) show stronger signals than those for shattering ( sh4) and plant architecture ( PROG1). Mating systems have been proposed to have a fundamental part in crop domestication42,43. We identified three major loci and two minor loci responsible for exserted stigma, and found that all the five QTLs were located within domestication sweeps.

Identification of domestication-associated variants

Characterization of domestication-associated genes and their genetic variation relies on high-quality genome sequences of both the cultivated species and its immediate progenitor44,45. Hence, we further sequenced the O. rufipogon accession W1943 with 100-fold genome coverage and carried out de novo genome assembly. The total length of the assembly is 406 megabases (Mb) and the N50 length of the initial contigs was 16 kb. After aligning the assembly against the reference genome of cultivated rice, we identified 2,621,077 SNPs, 619,132 small indels and 140,075 structural variants of large size ( Supplementary Table 14). We examined the sequence variants for their potential effects on protein coding, and identified a total of 128,010 non-synonymous SNPs and 49,236 sequence variants with large effect. Moreover, we surveyed the allele information at all the polymorphic sites across all rice accessions, and constructed a comprehensive sequence-variant frequency spectrum map.

We investigated the patterns of sequence polymorphism to detect variants that affect gene coding and have differential frequency in different populations. For the domestication loci with well-characterized genes, this approach allowed us to narrow down the causative polymorphism to a limited number of sequence variants ( Supplementary Figs 31-34). For the domestication loci for which known genes have not yet been identified, we found a total of 273 novel functional variants within 204 genes that showed high differentiation in allele frequencies between cultivated and wild rice ( Supplementary Table 15). Moreover, we identified 305 sequence variants in the promoter regions with a differential frequency between cultivated and wild rice ( Supplementary Table 16), and 1,120 functional variants that were fixed in either japonica or indica panel ( Supplementary Tables 17 and 18).

Discussion

Our study has provided new insights into how and where rice was likely to be domesticated, and we have identified a set of domestication sweeps and putative causal genes. Such endeavours will be enhanced by continuing improvements in assembly and annotation of wild rice genome sequences, generation of functional genomics data sets of wild rice46, advances in mapping of rice domestication traits47, and biological follow-up of the putative causal genes. The understanding of past domestication, including the selections on critical traits and the recent rapid speciation, will further guide future breeding efforts48. Moreover, the great diversity in the wild rice populations, which have much more natural allelic variation than domesticated rice, will further facilitate breeding to modify crops in the post-domestication era.

Methods Summary

The cultivated and wild rice accessions were all from large collections of rice accessions preserved at the China National Rice Research Institute in Hangzhou, China, and the National Institute of Genetics in Mishima, Japan. The DNA samples were sequenced on the Illumina Genome Analyzer IIx or HiSeq2000.

Online Methods

Sampling and sequencing

The cultivated and wild rice accessions were all from large collections of rice accessions preserved at the China National Rice Research Institute in Hangzhou, China, and the National Institute of Genetics in Mishima, Japan. The accessions were selected on the basis of the germplasm database records of phenotypic data and sampling localities to maximize genetic and geographic diversity. The collection was maintained by selfing in the laboratories. For each accession, genomic DNA from a single plant was used for sequencing, and seeds derived from the same plant were used for following field trials. In total, the genomes of 1,083 O. sativa accessions, 446 O. rufipogon accessions and 15 accessions of outgroup species were sequenced on the Illumina Genome Analyzer IIx generating 73-bp (or 117-bp) paired-end reads, each to approximately onefold (for O. sativa accessions), twofold (for O. rufipogon accessions) or threefold (for outgroup species) coverage. The detailed information, including geographic origin and sequencing coverage of the rice accessions, was listed in Supplementary Tables 2, 7 and 8. Library construction and sequencing of these accessions were performed as described21. One representative accession of O. rufipogon, W1943, was sequenced on the Illumina HiSeq2000, generating 100-bp paired-end reads with 100-fold genome coverage. An amplification-free method of library preparation49 was used in deep sequencing of the rice accession, which reduced the incidence of duplicate sequences, thus facilitating genome assembly and variation analysis.

Read alignment and SNP calling

The paired-end reads of all the rice accessions were aligned against the rice reference genome (IRGSP 4.0) using the software Smalt (version 0.4) with the parameters of '-pair 50, 700' and '-mthresh 50'. SNPs were called using the Ssaha Pileup package (version 0.5) with detailed procedure described previously21. Genotypes of the rice accessions, including 1,083 O. sativa accessions, 446 O. rufipogon accessions and 15 accessions of outgroup species, were further called at the SNP sites from the Ssaha Pileup outputs. The genotype calls in 15 accessions of outgroup species were used to determine the ancestral states of SNPs in O. sativa and O. rufipogon. SNPs in coding regions, which were defined based on the gene models in the RAP-DB (release 2), were then annotated to be synonymous or non-synonymous for calculating the non-synonymous/synonymous ratio and d N/d S ratio. The genotype data set of the 1,529 rice accessions (1,083 O. sativa accessions and 446 O. rufipogon accessions) was generated on the basis of the calls in each rice accession. Seven sets of genome sequences, which included bacterial-artificial-chromosome-based Sanger sequences and high-coverage resequencing data, were used to assess the accuracy of the genotype data sets ( Supplementary Table 6). The wild rice accessions with sequencing coverage >9 were selected to investigate the heterozygosity based on the overlapped reads that were aligned onto the reference sequence. For each accession, the proportion of heterozygosity genotypes was calculated at the polymorphic sites ( Supplementary Table 3).

Population genetics analysis

The software Haploview was used to calculate linkage disequilibrium with default settings, using SNPs with information in 446 O. rufipogon accessions50. Pairwise r ² was calculated for all the SNPs and then averaged across the whole genome. The matrix of pairwise genetic distance derived from simple SNP-matching coefficients was used to construct phylogenetic trees using the software PHYLIP51 (version 3.66). The software TreeView and MEGA5 were used for visualizing the phylogenetic trees. Principal component analysis of the SNPs was performed using the software EIGENSOFT52. The sequence diversity statistics ( π) and the population-differentiation statistics ( F ST) were computed using a 100-kb window. The value of π was calculated for each group in O. rufipogon and O. sativa, respectively, and the ratio of π in the full population (or each clade) of O. rufipogon to that in the full population (or corresponding subspecies) of O. sativa was used to detect selective sweeps. The genomic regions where both O. rufipogon and O. sativa show a low level of genetic diversity were excluded for further analysis. To adopt appropriate thresholds to reduce the false-positive rate but also retain true selection signals, thresholds were chosen on the basis of both whole-genome permutation tests and signals at known loci. Permutation tests were performed to estimate the genome-wide type I error rate and determine the threshold to call selective sweeps (see Supplementary Information section 2 for details)53. The method cross-population extended haplotype homozygosity (XP-EHH) was also tested for detecting selective sweeps using the software xpehh54 ( http://hgdp.uchicago.edu/Software/) ( Supplementary Fig. 20). The genetic distance between two clades was computed based on the matrix of pairwise genetic distance, where the distance of all pairs of accessions from the two clades were retrieved and averaged. A custom Perl script was developed to plot all O. rufipogon accessions, using the public geographic information of world borders from the 'Thematic Mapping' data set (version 0.3). The computational simulations under different demographic scenarios were performed using the program SFS_CODE40.

Planting, crossing and phenotyping

For the O. rufipogon population, approximately five seeds for each accession from the collection of wild rice were germinated and planted in the experimental field (in Sanya, China at N 18.65°, E 109.80°) from March 2011. The leaf sheath colour was observed and scored directly and the tiller angle was measured for each plant. The mapping population of 210 backcross inbred lines (BILs) and 61 chromosome segment substitution lines (CSSLs) was derived from a cross between O. sativa ssp. indica cv. Guangluai-4 and O. rufipogon accession W1943. The BILs were developed by one generation of backcross to Guangluai-4 followed with six generations of self-fertilization. The CSSLs were developed by five generations of backcross to Guangluai-4 followed with three generations of self-fertilization. Phenotyping was conducted in the experimental field (in Shanghai, China at N 31.13°, E 121.28°) from May to October, 2011. The fifteen traits that we phenotyped for this study include germination rate, tiller angle, heading date, stigma colour, the degree of stigma exsertion, plant height, panicle length, the degree of shattering, awn length, grain number per panicle, grain length, grain width, grain weight per 1,000 grains, hull colour and pericarp colour. The degree of stigma exsertion was scored based on the observation of ∼ 20 randomly sampled spikelets of each line, on a scale of 1-3 (no, incomplete or complete exertion). Seed germination rate was measured by using mature seeds which were placed in a plastic Petri dishes kept at 30 °C in the dark for 48 h5. Other traits were phenotyped and scored as described previously21,22,29.

Imputation and association analysis

For the genotype data set in O. rufipogon, genotypes of 446 O. rufipogon accessions were called specifically at the ∼ 5 million SNP sites that were polymorphic in the O. rufipogon population. In the panel for GWAS, only the SNPs that have a minor allele frequency (MAF) of more than 5% and contain genotype calls of more than 100 accessions were left for subsequent imputation. The k-nearest neighbour algorithm-based imputation method was used for inferring missing calls21. The specificity of the genotype data set before and after imputation was assessed using three sets of genome sequences ( Supplementary Table 6). Association analysis was conducted using the compressed mixed linear model55. The top five principle components were used as fixed effects and the matrix of genetic distance was used to model the variance-covariance matrix of the random effect. Permutation tests were used to define the threshold of association signals of the GWAS in the wild rice population. A total of 20 permutation analyses were performed (10 independent permutation tests for each of the two traits, sheath colour and tiller angle), which resulted in two 'association signals' with the thresholds we set53. Hence, there were an average of 0.1 false positives (that is, totally two false positives in 20 permutation tests) in a single whole-genome scanning analysis. Simulation tests were used to compare the performance of GWAS between the populations of cultivated and wild rice.

Genotyping and linkage analysis

Genomic DNA of each line in the mapping population was sequenced on the Illumina Genome Analyzer IIx, each to approximately 0.5× coverage. Both parents of the population, Guangluai-4 and W1943, were sequenced with at least 20× genome coverage, in a previous work21 and in this study, respectively. SNP identification between parents was conducted as described previously41. Genotype calling, recombination breakpoint determination and bin map construction was performed using the software SEG-Map ( http://www.ncgr.ac.cn/software/SEG/). QTL analysis of the fifteen traits was conducted with the composite interval mapping (CIM) method implemented in the software Windows QTL Cartographer56 (version 2.5) with a window size of 10 cM and a step size of 2 cM. QTL with LOD value higher than 3.5 were called, of which the location was described according to its LOD peak location. The phenotypic effect ( r ²) of each QTL was computed using Windows QTL Cartographer. QTLs located within selective sweep regions were further used to associate the selected regions with their biological functions. It needs to be noted that we adopted a stringent threshold in the QTL calling (LOD > 3.5), and the genomic regions with LOD ranging from 2.5 to 3.5 may include many minor QTLs (the threshold was set to 2.5 in most studies).

Genome assembly and contig anchoring

The genome of W1943 was assembled by using a custom pipeline integrating Phusion2 (clustering the raw reads into different groups)57 and Phrap (then assembling all the reads in each group to generate contigs)58. The N50 length of the entire assembly was calculated for the initial contigs with small contigs of <200 bp excluded. All the full-length complementary DNA sequences46 of W1943 were aligned with the final assembly of W1943 genome sequence using the software GMAP59 (version 6) with the parameters '-K 15000' and '-k 0.97'. The resulting contigs from whole-genome de novo assembly were anchored to the rice reference genome sequence (IRGSP4.0) using the software MUMmer60 (version 3).

Genome annotation and variant detection

Gene models of the genome of the wild rice W1943 were predicted using the software Fgenesh that was set for a monocot model61 (version 2.0). The resulting proteome of W1943 was compared with protein sequences in Rice Genome Annotation Project (version 7.0) using BLASTP with a cutoff of a minimum of 95% identity. Sequence variants, including SNPs, indels and imbalanced substitutions, were called using the diffseq program in the EMBOSS package62. Indels of large size were called from the alignment results of MUMmer. Effects of the sequence variants were predicted according to the gene models of Nipponbare in the RAP-DB (release 2) across the rice genome. For indels in genic regions and SNPs with large effect around the domestication loci, the effects were mainly based on the reference gene models.

Population-scale sequence comparison

The sequence reads of 1,083 O. sativa accessions, 446 O. rufipogon accessions and 15 outgroup accessions were then aligned against assembled genome sequences of W1943 using the same parameters with those against the reference Nipponbare genome sequences. Genotypes of each accession were called at all sequence variant sites (including SNPs, indels and imbalanced substitutions that were detected from assembled sequences), based on the alignment outputs against the two genome sequences. The allele frequencies at the sequence variant sites were calculated for each clade of O. sativa and O. rufipogon. In each clade, variant sites with information of less than 10 accessions (less than 2 for the outgroups) were then excluded for computing allele frequencies, namely no data available.

Accessions

Primary accessions

Data deposits

DNA sequencing data are deposited in the European Nucleotide Archive ( http://www.ebi.ac.uk/ena/) under accession numbers ERP001143, ERP000729 and ERP000106. De novo assembly and genome annotation of wild rice W1943, the genotype dataset of 1,529 rice accessions and the imputed dataset of 446 O. rufipogon accessions for GWAS are available at the Rice Haplotype Map Project database ( http://www.ncgr.ac.cn/RiceHap3).

References

1. 1.
Oka, H. I. Origin of cultivated rice. (Japan Scientific Societies Press, 1988)

2. 2.
Khush, G. S. Origin, dispersal, cultivation and variation of rice. Plant Mol. Biol. 35, 25-34 (1997)

3. 5.
Li, C. , Zhou, A. & Sang, T. Genetic analysis of rice domestication syndrome with the wild annual species, Oryza nivara. New Phytol. 170, 185-194 (2006)

4. 10.
Londo, J. P. , Chiang, Y. C. , Hung, K. H. , Chiang, T. Y. & Schaal, B. A. Phylogeography of Asian wild rice, Oryza rufipogon, reveals multiple independent domestications of cultivated rice, Oryza sativa. Proc. Natl Acad. Sci. USA 103, 9578-9583 (2006)

5. 23.
Ge, S. , Sang, T. , Lu, B. R. & Hong, D. Y. Phylogeny of rice genomes with emphasis on origins of allotetraploid species. Proc. Natl Acad. Sci. USA 96, 14400-14405 (1999)

6. 24.
Vaughan, D. A. , Morishima, H. & Kadowaki, K. Diversity in the Oryza genus. Curr. Opin. Plant Biol. 6, 139-146 (2003)

7. 28.
Tadege, M. et al. Reciprocal control of flowering time by OsSOC1 in transgenic Arabidopsis and by FLC in transgenic rice. Plant Biotechnol. J. 1, 361-369 (2003)

8. 31.
Andaya, V. C. & Tai, T. H. Fine mapping of the qCTS12 locus, a major QTL for seedling cold tolerance in rice. Theor. Appl. Genet. 113, 467-475 (2006)

9. 32.
Saitoh, K. , Onishi, K. , Mikami, I. , Thidar, K. & Sano, Y. Allelic diversification at the C (OsC1) locus of wild and cultivated rice: nucleotide changes associated with phenotypes. Genetics 168, 997-1007 (2004)

10. 38.
Sweeney, M. T. , Thomson, M. J. , Pfeil, B. E. & McCouch, S. Caught red-handed: Rc encodes a basic helix-loop-helix protein conditioning red pericarp in rice. Plant Cell 18, 283-294 (2006)

11. 40.
Hernandez, R. D. A flexible forward simulator for populations subjects to selection and demography. Bioinformatics 24, 2786-2787 (2008)

12. 42.
Chen, K. Y. , Cong, B. , Wing, R. , Vrebalov, J. & Tanksley, S. D. Changes in regulation of a transcription factor lead to autogamy in cultivated tomatoes. Science 318, 643-645 (2007)

13. 50.
Barrett, J. C. , Fry, B. , Maller, J. & Daly, M. J. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics 21, 263-265 (2005)

14. 51.
Felsenstein, J. PHYLIP: phylogeny inference package (version 3.2). Cladistics 5, 164-166 (1989)

15. 56.
Wang, S. , Basten, C. J. & Zeng, Z. B. Windows QTL Cartographer 2.5. (Department of Statistics, North Carolina State Univ., 2007)

16. 61.
Salamov, A. A. & Solovyev, V. V. Ab initio gene finding in Drosophila genomic DNA. Genome Res. 10, 516-522 (2000)

Download references

Acknowledgements

We thank the China National Rice Research Institute for providing all cultivated rice germplasm and Chinese wild rice accessions. The rest of the wild rice accessions were distributed from the Rice Collection of National Institute of Genetics jointly supported by the National Bioresource Project, MEXT, and Systems Functional Genetics Project of the Transdisciplinary Research Integration Center, ROIS, Japan. We thank Z. Ning, Y. Minobe and A. Osbourn for their advice and assistance in this work. This work was supported by the Ministry of Science and Technology of China (2011CB100205, 2012AA10A302 and 2012AA10A304), the Ministry of Agriculture of China (2011ZX08009-002 and 2011ZX08001-004), the National Natural Science Foundation of China (31121063) and the Chinese Academy of Sciences (to B.H.).

Author information

Xuehui Huang

, Nori Kurata

, Xinghua Wei

& Zi-Xuan Wang

These authors contributed equally to this work.

Affiliations

1. National Center for Gene Research, Institute of Plant Physiology and Ecology, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai 200233, China

Xuehui Huang

, Zi-Xuan Wang

, Ahong Wang

, Qiang Zhao

, Yan Zhao

, Kunyan Liu

, Hengyun Lu

, Wenjun Li

, Yunli Guo

, Yiqi Lu

, Congcong Zhou

, Danlin Fan

, Qijun Weng

, Chuanrang Zhu

, Tao Huang

, Lei Zhang

, Yongchun Wang

, Lei Feng

, Qilin Zhan

, Canyang Li

, Tingting Lu

, Qi Feng

& Bin Han

2. Plant Genetics Laboratory and Comparative Genomics Laboratory, National Institute of Genetics, Mishima, Shizuoka 411-8540, Japan

Nori Kurata

, Zi-Xuan Wang

, Hiroyasu Furuumi

, Takahiko Kubo

, Toshie Miyabayashi

, Asao Fujiyama

& Atsushi Toyoda

3. State Key Laboratory of Rice Biology, China National Rice Research Institute, Chinese Academy of Agricultural Sciences, Hangzhou 310006, China

Xinghua Wei

, Xiaoping Yuan

, Qun Xu

, Guojun Dong

& Qian Qian

4. National Center for Plant Gene Research, State Key Laboratory of Plant Genomics, Institute of Genetics and Developmental Biology, Chinese Academy of Sciences, Beijing 100101, China

5. Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing 100029, China

Contributions

B.H. conceived the project and its components. X.H. and B.H. designed studies and contributed to the original concept of the project. X.W., X.Y. and Q.X. contributed the collection of rice cultivars and Chinese wild rice accessions. N.K., H.F., T.K., T.M., A.F. and A.T. contributed the collection of other wild rice accessions and analysed geographical distributions of wild rice. Z.-X.W., A.W., Y.W., L.F., Qilin Z., C.L., G.D. and Q.Q. contributed in phenotyping of rice. Z.-X.W., A.W. and L.F. contributed in phenotyping of the backcross inbred line population and genetic mapping of domestication-related traits. W.L., Y.G., Y.L., C.Z., D.F., Q.W. and Q.F. performed the genome sequencing. X.H., Y.Z., K.L., C.Z., T.H., L.Z. and T.L. performed genome data analysis. Qiang Z. and H.L. performed de novo genome assembly. X.H. performed evolutionary study and genome annotation. X.H., Y.Z. and K.L performed GWAS, population genetics, and statistical analyses. J.L. contributed to functional analyses. X.H. and B.H. analysed whole data and wrote the paper.

Competing interests

The authors declare no competing financial interests.

Corresponding author

Correspondence to Bin Han.

Supplementary information

PDF files

1. 1.

This file contains Supplementary Text 1-8, Supplementary References and Supplementary Tables 1, 3-6, and 8114 (see separate zipped file for Supplementary Tables 2, 7, 15, 16, 17 and 18.

2. 2.

This file contains Supplementary Figures 1-34.

Zip files

1.

This file contains Supplementary Tables 2, 7 and 15-18.

About this article

Publication history

DOI
https://doi.org/10.1038/nature11532
Rights and permissions

This article is distributed under the terms of the Creative Commons Attribution-Non-Commercial-Share Alike licence ( http://creativecommons.org/licenses/by-nc-sa/3.0/), which permits distribution, and reproduction in any medium, provided the original author and source are credited. This licence does not permit commercial exploitation, and derivative works must be licensed under the same or similar licence.

To obtain permission to re-use content from this article visit RightsLink.

# data[].attributes.attachment_link.title
A map of rice genome variation reveals the origin of cultivated rice

# data[].attributes.attachment_link.description
Crop domestications are long-term selection experiments that have greatly advanced human civilization. The domestication of cultivated rice ( Oryza sativa L.) ranks as one of the most important developments in history. However, its origins and domestication processes are controversial and have long been debated.

# data[].attributes.body
Genetic control and comparative genomic analysis
of flowering time in Setaria (Poaceae)

Please put a special emphasize on QTL analyses

Presenters: Dominika Brodzisz, Charlie Clarke & John Savage

# data[].attributes.attachment_link.title
Genetic Control and Comparative Genomic Analysis of Flowering Time in Setaria (Poaceae)

# data[].attributes.attachment_link.description
We report the first study on the genetic control of flowering in Setaria, a panicoid grass closely related to switchgrass, and in the same subfamily as maize and sorghum.

# data[].attributes.body
Molecular characterization of the major wheat domestication gene Q

Please put a special emphasize on yeast two-hybrid assays

Presenters: Jorge Alejo, Alba Beltrán and Manuel Beltrán

# data[].attributes.attachment_link.readable_text
Plant materials and Southern analysis:

The genomic sequence of WAP2 was obtained from the genotypes listed in Table 1, where GenBank accession numbers for the WAP2 sequences are also presented. The sequence of the T. monococcum WAP2 gene was identified previously (). The CS fast neutron-induced speltoid deletion mutant fndel-143 (null for Q on chromosome 5A, but retains homeoalleles on 5B and 5D) () was used to ensure all PCR amplifications were specific to the chromosome 5A copy of WAP2 and to evaluate the specificity of the Taqman system (see Transcription analysis). The T. aestivum cv. "Bobwhite" ( QQ genotype) was used for transformation experiments. CS homeologous chromosome group 5 nullisomic-tetrasomic (NT) lines, where a pair of missing chromosomes is partially compensated by an extra pair of homeologous chromosomes, were used to assign restriction fragments to individual chromosomes by Southern analysis. DNA extraction, digestion, Southern transfer, and hybridization were performed as previously described ().

Generation of mutants:

Mutants were generated by treating CS seeds with 0.4% ethyl methanesulfonate (EMS) in 0.05 M phosphate buffer as described (). Selected M ₂ families were then planted and screened in a greenhouse for plants with speltoid spikes.

BAC sequencing and analysis:

A fragment of the T. monococcum WAP2 gene was used to screen the LDN BAC library (). Positive BACs were digested with HindIII and DraI, Southern blotted, and hybridized with the same probe to assign BACs to chromosome 5A or 5B on the basis of restriction fragment sizes. One BAC from chromosome 5A (376H15, 142 kb) and one from 5B (1004P5, 152 kb) were sequenced by Myriad Genetics (Salt Lake City). BAC sequences were subjected to BLASTn and BLASTx (Altshul et al. 1997) searches of the public databases to identify the WAP2 gene sequence from the BACs. The sequences were then submitted to RICEGAAS ( http://ricegaas.dna.affrc.go.jp/) to identify predicted open-reading frames using multiple gene prediction programs. Further details regarding the analysis of the chromosome 5B WAP2 homeoallele are included as supplementary material available at http://www.genetics.org/supplemental/.

WAP2 sequencing and analysis:

The T. monococcum WAP2 sequence was used to design primers to amplify WAP2 from the genotypes in Table 1. The genomic WAP2 gene sequence was PCR-amplified as three separate overlapping fragments using primers AP5P.11-3 (5′-GCCCTCGCAGCCCGCGGCCACCGCGCTCCCA) or AP2startF (5′-ATGGTGCTGGATCTCAATGTGGAGTCGCCGGCGGA) in combination with AP2.8R (5′-CGCGGCCAAATCGGGGCAAAGGAATTCAAACGA) for fragment 1, WAP2.2F (5′-CACTGGATAATTTCTTCAGGTGGTTTCGACACTGC) with AP2.15R (5′-ACATGGAACCTTAATTTCAGGAACGAACTTGTCG) for fragment 2, and AP2.16F (5′-CTGCTTGGTGCGCTGCTCCACCAGCTTACTGAAA) with AP45.1R (5′-CAGAAGGCCCAACGGTTAACGCAACAATGGC) for fragment 3. PCR conditions for fragment 1 were 250 ng Dra I-digested genomic DNA, 0.4 μm of each primer, 200 μm dNTP, 1× reaction buffer, 8% DMSO, and 2.5 units Herculase Hotstart DNA Polymerase (Stratagene, La Jolla, CA) in a 50-μl volume. Cycling conditions for fragment 1 amplification were the following: 98° for 3 min, 10 cycles of 98° for 40 sec, 68° for 30 sec, 72° for 2 min, 30 cycles of 98° for 40 sec, 68° for 30 sec, 72° for 2 min with extension time increasing 10 sec each cycle. One unit of Biolase DNA polymerase (Bioline USA, Randolph, MA) was added and allowed to finish cycling at 72° for 10 min. Fragments 2 and 3 were amplified using 1 μl Advantage cDNA polymerase mix (CLONTECH, Palo Alto, CA), 300 ng of DNA, 200 μm dNTP, 1× buffer, and 0.4 μm of forward and reverse primers in a 50-μl reaction. Cycling conditions for fragment 2 were the following: 94° for 4 min, 35 cycles of 94° for 40 sec, 68° for 30 sec, 72° for 3 min, and a 10-min final extension step. Cycling conditions for fragment 3 were identical to those of fragment 2 except that the annealing temperature was 72° and the extension time was reduced to 2 min.

Products from five independent PCR reactions for each segment of the gene from each genotype were separated on 1% agarose gels and gel-purified using the QIAquick gel extraction Kit (QIAGEN, Chatsworth, CA). Fragments were cloned using the TA cloning kit (Invitrogen, Carlsbad, CA) and sequenced. The sequences were reassembled using CLUSTAL W to obtain the full-length genomic sequences.

The coding sequence of WAP2 was determined by identifying database ESTs ( http://wheat.pw.usda.gov/NSF/; http://www.tigr.org/tigr-scripts/tgi/T_index.cgi?species=wheat) and 5′ RACE (gene-specific primer AP2race2: 5′-GGGCGGCGACGCGGGGAAGAGCTGCCTCGTG). The 5′ UTR was obtained using the BD SMART RACE cDNA amplification kit (CLONTECH) in combination with the BD Advantage-GC 2 PCR kit (CLONTECH). Reverse transcriptase (RT)-PCR for the three EMS-induced mutants, mq36, mq125, and mq194, was done with primers flanking the regions containing single base mutations (AP45.6F: 5′-ATGGGGCAGCAGGCCCCGGCGCCTGCGCCGATGGC and AP2.13R: 5′-CTCTTGGGATCGTGCGCGGTGGGTTGCGACATC), using the TITANIUM one-step RT-PCR kit (CLONETECH).

The 5′ region extending from position −140 to −1121 from the transcription start site was PCR amplified with primers AP5P.16F: 5′-GGATCACGTGGGTGGTTCTTTGTCCATGCC and AP5P.12R: 5′-GTCGGGGAGGCCAAGGGCATCAGAGG. Cycling conditions were the following: 94° for 4 min, 35 cycles of 94° for 40 sec, 66° for 30 sec, 72° for 2 min, and a 10-min final extension step.

Transformation:

Tissue culture, particle bombardment of immature embryos, and initial screening for the bar gene were done following previously described methods (; ; ). A 5.8-kb subclone containing the WAP2 gene and its native promoter from the LDN 5A BAC was ligated into the pSMART vector (Lucigen, Middleton, WI). The T. aestivum cv. Bobwhite, which is presumed to have two endogenous copies of the Q allele, was chosen for its ease of transformation. Calli were cobombarded with the subclone and the pAHC20 vector (contains the bar gene as a selectable marker driven by the maize ubiquitin promoter) using microprojectile bombardment. Positive individuals were grown to maturity and the spike morphology recorded. Selected T ₁ plants derived from speltoid and compactoid T ₀ plants were analyzed for transgene integration and transcription.

Relative quantitative-PCR:

Tissue for RNA extractions was collected, immediately frozen in liquid nitrogen, and stored at −80°. RNA was isolated using the QIAGEN RNeasy plant mini kit (QIAGEN) and on-column DNA digestion. cDNA was prepared using TaqMan RT reagents with random hexamer primers and multiScribe reverse transcriptase (Applied Biosystems, Foster City, CA) as suggested by the manufacturer. Relative quantitative (RQ)-PCR was performed on a 7500 Real Time PCR System (Applied Biosystems) using a TaqMan system for Q and either the 18S rRNA or actin genes. 18S TaqMan system (): left primer, CGGAGAATTAGGGTTCGA; right primer, CCGTGTCAGGATTGGGTA; probe, VIC-CTACCACATCCAAGGAAGGCAG-TAM. Q TaqMan system: left primer, CCCTGAATCGTCAACCACAATG; right primer, CCGTGCCATGTTGATGCA; probe, FAM-CTTCGTCCCAGTGGCCTG-NFQ. Actin Taqman system: left primer, ATGGAAGCTGCTGGAATCCAT; right primer, CCTTGCTCATACGGTCAGCAATAC; probe, VIC-CCTTCCTGATATCCACATCACACTTCATGATAGAGT-TAM. At least four replications of each reaction were performed using 1X TaqMan Universal PCR master mix (Applied Biosystems). Each reaction consisted of 5 μl of cDNA diluted 1:50, 100 mm of probe, and 200 mm of each primer. Templates to determine the amplification efficiency of 18S, actin, and Q consisted of five twofold dilutions. Raw C T values were averaged for each dilution and each sample. Amplification efficiencies of 18S, Q, and actin were tested for similarity by plotting RNA concentration vs. C T value. The difference between the slopes was <0.1 and thus considered to have similar amplification efficiencies. Sample averages were linearized using the 2 −ΔΔCT method ( Livak and Schmittgen 2001).

Yeast two-hybrid analysis for detecting dimerization of Q gene products:

The yeast two-hybrid MATCHMAKER system (BD Biosciences CLONTECH, Palo Alto, CA) was used for analysis of protein-protein interactions according to the manufacturer's protocols. First, the full-length coding sequences of Q from CS and q from CS-DIC 5A were amplified from cDNA using primers Q_Nde1_F: 5′-CATATGGTGCTGGATCTCAATGTGGAGTC and Q_BH1344_R: 5′-GGATCCTCAGTTGTCCGGCGGGCGGGGGAA. Serial deletions of the coding sequences were then constructed and ligated into vectors pGBKT7 and pGADT7 to be bait and prey, respectively. Yeast Y187 cells were transformed with the bait plasmids, and AH109 cells were used for the prey plasmids. The two yeast strains containing bait and prey were mating on YPD plates for overnight, and the yeast cells were cultured in SD medium without leucine (L) and tryptophan (W) for 2 days. To test the protein interactions, the SD/-L-W cultured yeast cells were inoculated on SD plates containing X-Gal but lacking leucine, tryptophan, histidine, and adenine. The empty vectors of pGBKT7 and pGADT7 were used as a negative control, and the interaction between p53 and LTA served as a positive control.

Phylogenetic analysis:

Phylogenetic trees were constructed from CLUSTALW alignments of the complete genomic DNA sequences of Q/q using the neighbor-joining method and multiple distance-based methods available in the MacVector v7.2 software. Confidence values for nodes were calculated using 1000 bootstraps.

Structure and validation of Q:

We obtained the genomic sequence of the wild-type WAP2 gene from CS using PCR, and determined the full-length coding region using a combination of RT-PCR, 5′-RACE, and searches against the wheat EST databases ( Figure 1A). The WAP2 gene consists of 10 exons and 9 introns, extends 3229 bp from start to stop codon, and has a GC content of 54%. The coding sequence is 1344 bp long plus a 128-bp 5′ UTR (supplemental Figure 1 at http://www.genetics.org/supplemental/) and a 255-bp 3′ UTR (supplemental Figure 2 at http://www.genetics.org/supplemental/). It encodes a protein of 447 amino acids.

The M ₂ generation of a population of CS EMS mutants was screened for the speltoid phenotype to identify putative knockouts. Three speltoid mutants harboring point mutations within WAP2 were identified ( Figure 1B). Mutant mq194 had a single base substitution in an AP2 DNA binding domain (exon 5) that resulted in the change of a cysteine to a tyrosine. The mutants mq125 and mq36 had point mutations in the donor site of intron 2 and the acceptor site of intron 7, respectively. RT-PCR experiments showed that the point mutations in mq36 and mq125 resulted in larger PCR products compared to the wild-type WAP2 gene, indicating alternate splicing in these mutants (data not shown). The altered transcription of the mutated WAP2 gene in mq36 produced a predicted protein longer than the wild type due to translation of intron 7, and mq125 likely produced a truncated protein due to a stop codon encountered in the translation of intron 2. The sequence analysis of the mutants validated that WAP2 is Q and therefore, the gene will be referred to as Q hereafter.

Structural comparison of Q and q alleles:

Previous research ( Kuckuck 1959) suggested that Q evolved as a duplication of q due to unequal crossing over. We performed Southern analysis by hybridizing a fragment of the Q gene with multiple nonfree-threshing and free-threshing hexaploid and tetraploid genotypes ( Figure 2). Restriction patterns and fragment hybridization intensities suggested that a single copy of the Q gene is present on each of the homeologous group 5 chromosomes. In addition, sequence analysis of the 142-kb BAC harboring the Q allele from LDN indicated a single ORF corresponding in size and structure to the q ORF from T. monococcum (data not shown).

To further investigate structural differences between Q and q alleles, we aligned the gene sequences from the 12 Triticum genotypes listed in Table 1. Comparison of the genomic sequences revealed six conserved differences between Q- and q-containing genotypes ( Figure 3). Four of these differences, including a variable microsatellite, were present in introns and one was in the 3′ UTR. One conserved nucleotide difference changed a predicted amino acid where, at position 329, all Q-containing genotypes possessed an isoleucine while all q-containing genotypes possessed a valine.

To identify regulatory elements in the promoter region of Q, we PCR amplified and sequenced a fragment extending from position −140 to −1121 from the transcription start site in the Q genotypes CS, T. turgidum ssp. carthlicum, and the Iranian spelt, and from q genotypes CS-DIC 5A, T. turgidum ssp. dicoccum, and both European spelts. We were unable to PCR amplify the 140-bp fragment immediately upstream from the transcription start site using whole genomic DNA as template. However, the sequence of this segment was available from the T. monococcum and LDN BACs, and we confirmed the sequence from these two genotypes by PCR amplification directly from the BACs followed by sequencing. We generated consensus sequences of the Q and q promoter regions on the basis of ClustalW alignments to identify conserved differences. Five conserved differences were identified (supplemental Figure 3 at http://www.genetics.org/supplemental/), but searches of plant promoter databases PlantCARE ( http://intra.psb.ugent.be:8080/PlantCARE/index.html), PlantProm ( http://www.softberry.com), and PLACE ( http://www.dna.affrc.go.jp/PLACE/) indicated that only two of these differences are likely to be associated with putative regulatory elements.

Transcription analysis:

We compared the level of transcription between the hexaploids CS ( Q), CS-DIC 5A ( q), and fndel-143 (null for Q) in various stages of spike development. Spikes were collected and assigned fraction-length values on the basis of the length of corresponding fully matured spikes. Transcription in both CS and CS-DIC 5A peaked in the early stages of spike development and gradually declined as the spikes matured ( Figure 4). While transcription patterns for both Q and q were quite similar, the transcription level of Q was consistently higher than that of q.

To investigate transcription levels in other tissues, we compared CS ( Q), CS-DIC 5A ( q), and fndel-143 (null for Q) using root tips, the youngest fully expanded leaf from 1-month-old plants (Feekes stage 4-5), and flag leaves from adult plants (Feekes stage 10). Both CS and CS-DIC 5A exhibited transcription in all tissues tested ( Figure 5), but levels of transcription were less compared to immature spikes. Higher levels of transcription were observed in CS relative to CS-DIC 5A in all tissues.

Yeast two-hybrid analysis:

The yeast two-hybrid system was used to test for potential homodimer formation and to investigate the biological significance of the point mutation responsible for changing the valine in q to an isoleucine in Q at amino acid position 329 (I329V). The results indicated that the full-length Q protein (Q1344) has the ability to form a homodimer ( Figure 6). An N-terminal deletion construct (Q800) containing the C-terminal 277 amino acids and the putative coil-forming region, which links the two AP2 domains, also had the ability to form a homodimer. Further deletions of the remaining C-terminal 247 amino acids lacking the coil region significantly decreased the level of homodimerization. The point mutation (I329V) of q greatly reduced homodimer formation of the full-length q protein (q1344), and it completely disrupted the formation of a q800 homodimer. In addition, while Q800 was able to form a homodimer, it was not able to form a dimer with q800 (data not shown).

Genetic transformation and investigation of Q dosage effects:

We transformed Bobwhite with Q from the durum wheat variety LDN and obtained 17 putative transformants. We screened the T ₀ generation for individuals with either speltoid or compactoid spikes. One speltoid T ₀ plant was identified, and five T ₁ progeny derived from this plant were grown for further analysis. In addition to having speltoid spikes ( Figure 7A), the T ₀ and all five T ₁ plants were taller than Bobwhite ( Figure 7E), had keeled tenacious glumes that adhered strongly to the seed, and fragile rachises, in which the disarticulation pattern was such that the rachis broke above the junction of the rachis and rachilla leaving a portion of the rachis at the base of the spikelet ( Figure 7F). This tendency mimics the disarticulation pattern of the shatter-prone wild wheats that lack the Q allele. Southern analysis showed each T ₁ individual harbored an abundance of transgene copies in addition to the endogenous gene (data not shown).

Two partially sterile T ₀ plants with compactoid spikes were identified. Three T ₁ plants derived from one compactoid T ₀ and two T ₁'s from a second compactoid T ₀ showed a range of phenotypes. One plant was of normal height and had square spikes ( Figure 7B), two plants were slightly shorter than Bobwhite ( Figure 7E) and had subcompactoid spikes ( Figure 7C), and two were very short ( Figure 7E) with compactoid spikes ( Figure 7D). All plants had tough rachises and round, soft glumes that loosely held the seed, and were free-threshing ( Figure 7F). Southern analysis indicated the T ₁ plants segregated for copy number and either were lacking the transgene (square spike) or possessed multiple copies of the transgene in addition to the endogenous gene (subcompactoid and compactoid spikes) (data not shown).

RQ-PCR indicated the level of Q transcription in the T ₁ speltoid plants was less than that observed in untransformed Bobwhite ( Figure 8). T ₁ plants with subcompactoid and compactoid spikes had greater levels of Q transcription compared to Bobwhite.

Phylogenetic analysis:

Phylogenetic analysis of the 12 Triticum genotypes indicated that all Q genotypes were highly similar, while q genotypes were more divergent ( Figure 9). The genomic sequence of the AP2-like homeoallele from T. turgidum ssp. durum cv. LDN chromosome 5B was used to root the tree.

View Abstract

# data[].attributes.attachment_link.title
Molecular Characterization of the Major Wheat Domestication Gene Q

# data[].attributes.attachment_link.description
The Q gene is largely responsible for the widespread cultivation of wheat because it confers the free-threshing character. It also pleiotropically influences many other domestication-related traits such as glume shape and tenacity, rachis fragility, spike length, plant height, and spike emergence time. We isolated the Q gene and verified its identity by analysis of knockout mutants and transformation.

# data[].attributes.body
Genome-wide association mapping reveals
a rich genetic architecture of complex traits
in Oryza sativa

Please put a special emphasize on GWAS

Presenters: Andrew Tighe, Daniel James, Roisin Glancy.

# data[].attributes.attachment_link.readable_text
Diversity panel and genotyping array

A rice diversity panel consisting of 413 inbred accessions of O. sativa collected from 82 countries ( Fig. 1; Supplementary Data 1) was genotyped using an Affymetrix single nucleotide polymorphism (SNP) array containing 44,100 SNPs (hereafter referred to as the 44 K chip). With a genome size of 380 Mb (ref. 13), this custom-designed genotyping chip provides high quality data (less than 4.5% missing data), with 1 SNP per 10 kb across the 12 chromosomes of rice (Methods; Supplementary Data 2). The diversity panel was evaluated for 34 traits related to plant morphology, grain quality, plant development and agronomic performance using field-grown plants with replications within and between years ( Supplementary Table S1; Supplementary Data 3).

Figure 1: Population structure in O. sativa.

( a) The large pie chart summarizes the distribution of subpopulations in the 413 O. sativa samples in our diversity panel, and the smaller pie charts on the world map correspond to the country-specific distribution of subpopulations sampled (note: large countries such as China, India and the US were divided into several major rice growing regions). The size of the pie chart is proportional to the sample size and colours within each pie chart are reflective of the percentage of samples in each subpopulation. Seeds representing each subpopulation are displayed with and without hull in the centre, with 1 cm scale bar. ( b) Principal component analysis was used to provide a statistical summary of the genetic data, and the top four principle components are illustrated in the bottom panels.

Population structure and linkage disequilibrium estimation in rice

Using principle component analysis (PCA) to summarize global genetic variation in the diversity panel, we observed clear, deep subpopulation structure in this collection of germplasm ( Fig. 1a). The top four principal components (PCs) explained almost half of the genetic variation ( Fig. 1b). The five subpopulations indica, aus, temperate japonica, tropical japonica and aromatic formed clear clusters based on the top four PCs, and were well differentiated from each other, with pairwise Fst (F-statistic) values ranging from 0.23-0.53. This is in agreement with previous findings where global germplasm collections have been used in combination with much smaller numbers of SNP or simple sequence repeat (SSR) genotypes. Because the array was designed to assay variation in all O. sativa groups, most SNPs are shared or polymorphic across subpopulations ( Table 1).

Table 1: Polymorphism summary of Affymetrix 44 K SNPs in each subpopulation.

We examined allele sharing across the panel by calculating 'identity by state' coefficients among all pairs of accessions ( Fig. 2a). We find that whereas allele sharing clearly tracks subpopulation ancestry as identified by the PCA analysis, there is also a substantial number of admixed accessions, highlighting the complex history of rice varieties grown throughout the world. Excluding the small sample of aromatic accessions, the mean observed identical by state (IBS) sharing is greatest between the closely related tropical japonica and temperate japonica accessions (0.80), followed by indica and aus (0.64), with relatively little IBS sharing between the two major subspecies, Indica and Japonica (0.47) ( Fig. 2a). The fact that most of the admixture occurs within (rather than between) subspecies underscores the existence of genetic and cultural barriers to genetic exchange between these two major groups of Asian rice, despite documented cases of targeted Japonica-Indica introgression mediated by artificial selection.

Figure 2: Identity by State and phenotypic variation among subpopulations.

( a) Individuals are ordered according to their genotypic distance (1-IBS, identified by state) clustering with the tree shown on the right. The upper diagonal shows the IBS-sharing between individuals (values rescaled from 0 to 1). The lower diagonal shows the individual correlation coefficients based on all phenotypes. Coloured bars along the bottom of the panel reflect the sample subpopulation assignment as labelled; dark colour within each subpopulation indicates admixed individuals. ( b) Summary of phenotypic distributions among all individuals, with phenotypes grouped by trait category and individuals grouped by subpopulation as in ( a).

The amount of genomic variation tagged by our SNP array was calculated by measuring the pairwise SNP linkage disequilibrium (LD) among the 44 K common SNPs. On average, LD drops to almost background levels around 500 kb-1 Mb, reaching half of its initial value at 100 kb in indica, 200 kb in aus and temperate japonica, and 300 kb in tropical japonica ( Supplementary Fig. S1). Given that our average inter-marker distance is 10 kb, we expect to have reasonable power to identify common variants of large effect associated with our traits of interest, even if we have not queried the causal variant for association in the domesticated varieties.

Phenotypic variation

The phenotypes we examined in our GWAS can be classified broadly into six categories: plant morphology-related traits; yield-related traits; seed and grain morphology-related traits; stress-related phenotypes; cooking, eating and nutritional-quality-related traits; and plant development, represented by flowering time, which we measured in three geographic locations that differed in day-length and ambient temperature. Canonical correlation analysis demonstrated that phenotypes within a category are often correlated, ranging from a low of −0.41 between brown rice seed width and brown rice seed length, to a high of 0.9 between hulled and dehulled seed morphology ( Fig. 2b; Supplementary Fig. S2).

For all the phenotypes evaluated in this study, we observed global similarities among members of the same subpopulation, consistent with the domestication and breeding history of these varieties. Correlation coefficients between accession pairs across all phenotypes were significantly higher for accession pairs from the same subpopulation than from different subpopulations ( P<2.2e−16, one-sided Mann-Whitney U-test) (lower triangle of Fig. 2a). Consistent with this observation, the top four PCs (based on the 44 K SNPs mentioned above) explained a large proportion of phenotypic variation, with values ranging from 20-40% ( Supplementary Table S1). In the case of rice grain, morphological and cooking-quality traits are key to varietal identity and have been under strong diversifying selection by humans in different parts of the world. Physical grain characteristics in rice are salient because they serve as indicators of local and regional eating preferences in a crop that, unlike wheat or maize, is consumed largely as whole kernel. Traits such as flowering time and disease resistance are also strongly correlated with region and environment, meaning that genotypic, phenotypic and environmental variation in O. sativa are all correlated to some degree, posing significant challenges for GWAS.

The strong confounding effect of population structure

The results of our genome-wide association scans are summarized in Supplementary Figures S3-S36 where we show SNP-trait associations discovered in the diversity panel as a whole, as well as in each subpopulation individually. As can be seen in the quantile-quantile plots ( Fig. 3a; Supplementary Figs. S3-S36), the distribution of observed −log10 P-values from the naïve analysis (no population structure adjustment) departed quite far from the expected distribution under a model of no association (that is, the P-values should lie on the diagonal line), with significant inflation of nominal P-values leading to a high level of false positive signals. Use of a modified mixed model strategy allowed us to consider different levels of population structure and relatedness in our diversity panel. This effectively eliminated the excess of low P-values for most traits, but it also likely eliminated true positives. This is a common problem seen in other systems as well; for example, geographic coordinates correlate closely with flowering time in plants. For this reason, we believe a combination of naïve and population structure-adjusted hits, coupled with subpopulation-specific analyses in rice, is the most thoughtful way to identify potential variants for follow up.

Figure 3: Phenotypic distribution and genome-wide association scan for plant height.

( a) Quantile-Quantile plots for both naïve and mixed model for plant height in all samples. ( b) Boxplot showing the differences in plant height among subpopulations. Box edges represent the upper and lower quantile with median value shown as bold line in the middle of the box. Whiskers represent 1.5 times the quantile of the data. Individuals falling outside the range of the whiskers shown as open dots. ( c) Histogram of plant height in all samples. Dashed black line represents the null distribution. ( d) Genome-wide P-values from the mixed model and naïve method. x axis shows the SNPs along each chromosome; y axis is the −log ₁₀ ( P-value) for the association. Coloured dots in ( a) and ( c) indicate SNPs with P-values <1×10 −4 in the mixed model and the top 50 SNPs in the naïve method; SNPs within 200 kb range of known genes are in red; other significant SNPs are in blue. Candidate gene locations shown as red vertical dashed lines with names on top.

Using the mixed model to analyse the associations between 34 phenotypes and 44 K SNP genotypes evaluated in our 413 O. sativa rice lines, we successfully identified both known associations (for example, enrichment in a priori candidate genes and previously reported QTLs from rice and other species) as well as new candidate loci in the rice genome. Detailed results for each of the 34 phenotypes can be found in the Supplementary Data 3 as well as online in the Gramene database (www.gramene.org) and on our project website (www.ricediversity.org/44kgwas).

Trade offs between the mixed model and naïve model

Plant height is an important developmental and yield-related trait. Dozens of genes regulating plant height in rice have been identified previously including dwarfing mutants, QTLs, orthologues from other plant species, and genomic targets of fine-mapping experiments related to harvest index and yield. Both the naïve and the mixed model consistently detected strong signal linked to the Green Revolution semi-dwarf gene, SD1, on chromosome 1 ( Fig. 3d). Interestingly, several SNPs near other height-controlling genes such as OsBAK1 on chromosome 8 (ref. 27), DGL1 on chromosome 1 (ref. 28) were only detected by the naïve approach ( Fig. 3d). This suggests that, in the case of rice, the mixed model may overcompensate for population structure and relatedness, leading to false negatives. Therefore, the many mapping resources derived from crosses between parents belonging to different subpopulations and Oryza species will be needed to complement GWAS, helping to reduce the rate of false positives and false negatives, yielding QTLs that cannot be identified by mapping within subpopulations.

Genetic heterogeneity across subpopulations

In our diversity panel, aromatic varieties had the longest mean panicles (30 cm), temperate japonica had the shortest (21 cm), aus and indica had intermediate panicle length, and the greatest range of panicle length was observed among tropical japonica varieties ( Fig. 4a).

Figure 4: Genetic heterogeneity of panicle length across subpopulations.

( a) Histogram showing distribution of panicle length across the diversity panel and boxplot showing differences in panicle length among subpopulations. In boxplot, the box edges represent the upper and lower quantile with median value shown as bold line in the middle of the box. Whiskers represent 1.5 times the quantile of the data. Individuals outside of the range of the whiskers shown as open dots. ( b) Genome-wide P-values from the mixed model for panicle length for all 413 accessions in top panel ( all), and for tropical japonica, temperate japonica, indica and aus subpopulations individually in subsequent panels. Note: the aromatic subpopulation was not included because of the small sample size. X-axis indicates the SNP location along the 12 chromosomes; y axis is the −log ₁₀ ( P value) from each method. Coloured dots indicate SNPs with P-values <1×10 −4 in the mixed model; SNPs within 200 kb range of known genes are in red; other significant SNPs are in blue. Candidate genes near peak SNP regions known to be previously associated with panicle, stem and internode elongation in rice are shown along the top.

To determine whether different networks of alleles were associated with trait variation in the different subpopulations, we performed GWAS on each subpopulation independently and in the panel as a whole, and compared results. As summarized in Figure 4a,b, the genetic architecture of panicle length differs significantly among subpopulations and different GWAS peaks are observed when the subpopulations are analysed individually or when the diversity panel is analysed as a whole. For example, in the indica population, we see clusters of highly significant SNPs near OsTB1 [ TEOSINTE BRANCHED1 (ref. 30)], SLR1 [ SLENDER RICE1 (ref. 31)] and OsBRI1 [syn. DWARF61, or BRASSINOSTEROID-INSENSITIVE1 (ref. 32)], in the aus subpopulation we observe significant SNPs near FZP [ FRIZZY PANICLE] and SSD1 [ SWORD SHAPE DWARF1 (ref. 34)], and in the tropical japonica population, we see SNPs near OsLIC [LEAF AND TILLER ANGLE INCREASED CONTROLLER] and MOC1 ( MONOCLUM 1 (ref. 36)).

From these results, we conclude that different networks of genes regulate panicle length in different subpopulations and propose that subpopulation-derived genetic heterogeneity is a general pattern in O. sativa. This suggests that the Indica and Japonica varietal groups should be properly treated as true sub-species for association analyses, and helps explain why crosses between members of divergent subpopulations, as well as between cultivated and wild species, often give rise to transgressive offspring. We also demonstrate that the subpopulations of O. sativa contain alleles with vastly different effect-size on many traits of interest (that is, allele effects that are in the opposite direction to mean subpopulation differences for those traits). This conforms to the general mechanism that explains the production of extreme, or transgressive, phenotypes at both the species level and below and suggests a blueprint for harnessing natural variation to liberate transgressive phenotypes in the context of plant improvement.

Genotype by environment effects

To investigate how environmental variation affected the performance of GWAS, we evaluated flowering time in three different environments and compared results. One experiment was conducted during 2007 in the field in Stuttgart, Arkansas, USA (34°4′) under long-day conditions (14-12 h during May-September); one was conducted in the field in Faridpur, Bangladesh (23°5′) under 12-13 h days (January-May); and the third was conducted in the greenhouse in Aberdeen, Scotland, UK (57° 9′) across a nine-month period during which the days became very long and then very short (a range of 18-6 h during the period spanning March-December). The GWAS peaks explained between 5-50% of the phenotypic variation for flowering time in each environment ( Supplementary Data 4). As seen in Figure 5a, 10 genomic regions were associated with candidate genes for flowering time under one or more daylengths while only the HEADING DATE 1 ( HD1) region on chromosome 6 was detected in more than one environment.

Figure 5: Genome-wide association scan for flowering time.

( a) Genome-wide P-values from the mixed model for flowering time in three geographic locations are shown in the three panels. Association analysis in each subpopulation is shown in each row of the matrix. X axis indicates the SNP location along the 12 chromosomes, with chromosomes separated by vertical grey lines; y axis is the −log ₁₀ ( P value) from each method. Candidate genes previously shown to determine flowering time near peak SNPs are shown along the top, rice genes are in red, Arabidopsis homologues are in black. SNPs with P value <1×10 −4 are indicated by coloured dots. SNPs within 200 kb range of known rice flowering time genes are in red; SNPs within 200 kb range of Arabidopsis flowering-time homologues are in magenta; other significant SNPs are in blue. ( b) GWAS regions associated with photoperiod sensitivity, calculated as the ratio of days-to-flowering across pairs of environments.

The most significant signal was observed under very long days in Aberdeen around HD1, the major photoperiod-sensitivity locus, (synonym: SE1, or OsCONSTANS, OsCO) on chromosome 6 (ref. 39). A well-defined peak in the same location was observed under long days in Stuttgart, AR when either the entire diversity panel or the temperate japonica subpopulation was analysed. The significant SNPs detected in Aberdeen covered an extensive region of 2.3 Mb around HD1, corresponding to a 'mountain range' as described by Atwell et al. The 'mountain range' distribution may be due to the presence of several linked genes that contribute to flowering time across the region, and/or to the presence of multiple alleles at the HD1 locus, along with multiple introgression events that have been documented within a 5.5 Mb region around the HD1 gene. In domesticated species like rice, loci that are critical to both local adaptation and yield performance are often the targets of both natural and artificial selection, leading to complex forms of allele sharing and admixture in diverse varieties.

Some varieties were highly sensitive to daylength and others, mostly temperate japonica accessions, were insensitive to photoperiod and flowered at similar times across the three environments. When photosensitivity (the ratio of days-to-flowering across pairs of environments) was used as a derived trait for GWAS, the most significant SNPs ( P<10 −6) were consistently found near HD1, followed by a region on chromosome 7 containing homologues of Arabidopsis genes known to regulate the circadian rhythm [that is, TIME FOR COFFEE, TIC] and light sensing [ PLASTIC MOVEMENT IMPAIRED 15, PMI15] ( Fig. 5b), which have not previously been shown to be associated with natural variation for flowering time in rice.

We also demonstrate the effect of genotype by environment (GxE) interaction for flowering time by comparing GWAS results over two years in the same location in Stuttgart, Arkansas ( Supplementary Fig. S37). In this case, we observe extensive year-to-year variation between 2006 and 2007, mostly due to the different weather patterns experienced during the two growing seasons. Several GWAS peaks associated with candidate genes are significant in only a single year. The HD1 locus was significant in 2007, but not in 2006, though it was significant using the average flowering time from two years.

Although the genetic complexity and low heritability of flowering time as well as several other traits evaluated here in field grown plants ( Supplementary Figs S37-S40) tend to confound the interpretation of GWAS results, this study provides an opportunity to look carefully at a range of ecologically and agronomically important traits evaluated under natural growing conditions and compare GWAS results with prior QTL and mutant studies to better understand plant growth and development.

Gene linkage or pleiotropy

A matrix summarizing the QTL regions associated with all traits, as well as the percent of the phenotypic variation explained by significant SNPs for each trait, can be found in Figure 6. For many traits, the maximum-effect locus falls within a 200 kb region containing a previously identified functional gene (as highlighted in Fig. 6b and summarized in Supplementary Data 4). When our results are compared with those of Huang et al., the same known genes showed clear signal for the same phenotypes (for example, GS3 and qSW5 for grain length and width, SSII-3 and Waxy for alkaline spreading value and amylose content, Rc for pericarp colour). The significant SNPs in our study explained up to 58% of the variance compared with values up to 68% reported by Huang et al. In addition, we evaluated traits not previously documented and identified known genes associated with those traits (for example, SD1 for plant height, OsMADS13 for flowering time and Pi-ta for blast resistance). This demonstrates that our 44 K SNP array is capable of capturing the major common variants responsible for critical agronomic traits.

Figure 6: Summary of trait associations across genomic regions and percentage of variance explained by significant locus.

( a) Each row represents a trait, and each column corresponds to a genomic region containing multiple SNPs that are significantly associated with a trait. Significance is colour-coded based on the P value of the association. ( b) The x axis represents the trait, the y axis shows the contribution (%) of significant loci. Candidate genes detected within 200 Kb region of significant loci are labelled on top of the maximum effect locus.

In several cases, the same SNPs were significantly associated with multiple traits. This could be the result of pleiotropy or closely linked genes (local LD). For example, we observed SNPs at 31 Mb chromosome 4 that were significantly associated with both rice blast disease resistance and flag leaf width, and SNPs associated with rice blast disease resistance, amylose content and flowering time at 4.2-4.6 Mb on chromosome 6. These associations were also supported by Canonical correlation analysis based on traits measured in Arkansas ( Supplementary Fig. S2, r=−0.3 for blast resistance and flag leaf width, r=−0.31 for blast resistance and flowering time, r=0.37 for amylose content and flowering time). Similar trait associations have been previously reported in these and other regions in rice. Linkage among favourable alleles is a strong determinant of phenotypic value under both natural and artificial selection, a fact long appreciated by plant breeders. Validation studies involving joint linkage and association mapping, coupled with fine-mapping to identify the exact genes and alleles underlying our GWAS hits, will be required to more clearly understand the relationship between these candidate genes and the phenotypes observed in our panel, as well as to provide breeders with the appropriate genomic tools needed to break deleterious linkages and liberate valuable alleles in this region.

# data[].attributes.attachment_link.title
Genome-wide association mapping reveals a rich genetic architecture of complex traits in Oryza sativa

# data[].attributes.attachment_link.description
Understanding the genetics and physiology of domesticated species is important for crop improvement. By studying natural variation and the phenotypic traits of 413 diverse accessions of rice, Zhao et al. identify many common genetic variants that influence quantitative traits such as seed size and flowering time.

# data[].attributes.body
Duplication of a well-conserved homeodomain-leucine zipper transcription factor gene in barley generates a copy with more specific functions.

Please put a special emphasize on phylogenetic analyses

Presenters: Anthony Rowan, Dianah Sithole & Dearbháile de Buitleár

# data[].attributes.body
Six-rowed barley originated from a mutation in a
homeodomain-leucine zipper I-class homeobox gene

Please put a special emphasize on in situ hybridization

Presenters: Sirui Pan, Katarzyna Sleczka & Abdulaziz Albogami

# data[].attributes.attachment_link.readable_text
Abstract

Increased seed production has been a common goal during the domestication of cereal crops, and early cultivators of barley ( Hordeum vulgare ssp. vulgare) selected a phenotype with a six-rowed spike that stably produced three times the usual grain number. This improved yield established barley as a founder crop for the Near Eastern Neolithic civilization. The barley spike has one central and two lateral spikelets at each rachis node. The wild-type progenitor ( H. vulgare ssp. spontaneum) has a two-rowed phenotype, with additional, strictly rudimentary, lateral rows; this natural adaptation is advantageous for seed dispersal after shattering. Until recently, the origin of the six-rowed phenotype remained unknown. In the present study, we isolated vrs1 ( six-rowed spike 1), the gene responsible for the six-rowed spike in barley, by means of positional cloning. The wild-type Vrs1 allele (for two-rowed barley) encodes a transcription factor that includes a homeodomain with a closely linked leucine zipper motif. Expression of Vrs1 was strictly localized in the lateral-spikelet primordia of immature spikes, suggesting that the VRS1 protein suppresses development of the lateral rows. Loss of function of Vrs1 resulted in complete conversion of the rudimentary lateral spikelets in two-rowed barley into fully developed fertile spikelets in the six-rowed phenotype. Phylogenetic analysis demonstrated that the six-rowed phenotype originated repeatedly, at different times and in different regions, through independent mutations of Vrs1.

Results

Map-Based Cloning of Vrs1.

We previously mapped the vrs1 locus to a 0.90-cM interval between cMWG699 and MWG865 ( 15, 16). In the present study, we narrowed the locus to a 0.07-cM interval between e40m36-1110S and BC12348 by using segregating progeny equivalent to 9,831 gametes ( Fig. 1 J). The candidate genomic region was covered completely with six bacterial artificial chromosome (BAC) clones by means of chromosome walking. The contig containing vrs1 is composed of 518,343 bp. Annotation showed three predicted genes [see supporting information (SI) Table 1]: HvHox1 and HvEP2 appeared to be intact genes, whereas HvEP1 is highly degenerated and interrupted by several insertions of transposable elements. HvHox1 is the only gene that lies between the flanking markers e40m36-1110S and 53F18-T7 and is thus a likely candidate for Vrs1 ( Fig. 1 J). The ORF of two-rowed barley encoded a polypeptide composed of 222 amino acid residues, including a homeodomain-leucine zipper motif (HD-ZIP) (see SI Fig. 5).

Expression of Vrs1 in Lateral Spikelet Primordia.

The expression pattern of the HD-ZIP I gene ( Vrs1) was highly tissue- and stage-specific in two-rowed barley. Transcription of Vrs1 was abundant during the early developmental stages of the immature spike (1-4 mm long; Fig. 2 A), when the central and lateral spikelet primordia became differentiated. (For a description of the stages of spike development, see ref. 17.) Transcript levels were lower, but still moderately abundant, at a length of 5 to 10 mm, but decreased greatly at later stages. In situ hybridization revealed that Vrs1 was expressed only in the lateral spikelet primordia ( Fig. 2 B-F). Vrs1 was not detectable at the double-ridge stage, when the central and lateral spikelets remain undifferentiated ( Fig. 2 D), but Vrs1 expression was clearly detected at the triple-mound stage ( Fig. 2 E) and at the glume primordium stage ( Fig. 2 F). During these Vrs1-expressing stages (when the immature spikelets are 1.5 to 2.0 mm long), the primordium of the central spikelet became larger than the primordia of the lateral spikelets. Only the lateral spikelet primordia showed expression of the gene ( Fig. 2 B); the central spikelet primordia showed no expression ( Fig. 2 C).

Mutant Vrs1 Gene.

To demonstrate the biological function of Vrs1, we analyzed mutant lines. Among 91 mutant lines reported ( 14), we included 57 mutant lines ( Fig. 3), which were derived from five two-rowed cultivars, of which both original and mutant seeds were available. The hexastichon ( hex-v) mutants have six-rowed spikes with fully fertile, well developed, and long-awned lateral spikelets ( 18) and thus resemble normal six-rowed barley. The Intermedium spike-d ( Int-d) mutants produced sterile or partially fertile lateral spikelets with variable awn length, which appeared intermediate to those of two-rowed and six-rowed barley ( 18). Because allelism of these mutations with vrs1 was confirmed in a previous study ( 19), morphological changes in the mutants can be attributed to changes in the Vrs1 gene. Lesions in Vrs1 were correlated with morphological changes in 48 mutant lines ( Fig. 3 A). Of the 48 mutants, 21 revealed a single amino acid substitution, 12 revealed truncation of the protein by a new stop codon, 3 showed a single nucleotide substitution in the conserved splicing sites of introns (change of splicing was confirmed by sequencing their transcripts), 5 had a frameshift mutation caused by a deletion, and 7 revealed complete deletion of the Vrs1 region (based on analysis using many flanking markers; see SI Table 2). These deletions (>182 kb), which were generated by means of irradiation, always resulted in hex-v-type six-rowed spikes under a range of growing conditions (see SI Table 3). The phenotypes observed in mutants that consistently exhibited six-rowed spikes support our hypothesis that complete deletion of Vrs1 occurred. Because the 7 deletion mutants did not show any developmental lesions, Vrs1 appears to be dispensable in barley. We found 9 mutants without any DNA changes throughout the coding region, and RT-PCR analysis showed that 7 of the 9 exhibited reduced or no expression of Vrs1 ( Fig. 3 B). This result suggests the occurrence of mutational events in the regulatory regions for Vrs1. The remaining 2 mutants showed Vrs1 expression at almost the same level as in the two-rowed cultivars, suggesting that posttranscriptional regulation is involved in these mutant phenotypes. The two mutation lines ( hex-v.08 and hex-v.46) were derived from cv. 'Bonus' by ethyleneimine exposure and neutron bombardment, respectively, and were allelic to vrs1, as confirmed by crosses with hex-v.3 and hex-v.4 ( 19).

The Origin of Six-Rowed Barley.

To infer the origin of six-rowed barley, we analyzed the sequence of the Vrs1 region from 15 two-rowed and 16 six-rowed cultivars, one variety ( deficiens), and three wild barley lines chosen from different geographical areas. Two alleles in two-rowed barley ( Vrs1.b2 in 1 cultivar and Vrs1.b3 in 14 cultivars) and three alleles in six-rowed barley ( vrs1.a1 in 10 cultivars, vrs1.a2 in 5 cultivars, and vrs1.a3 in 1 cultivar) were identified by means of haplotype analysis ( Fig. 4 A). The sequences of each allele were identical except that vrs1.a1 had three subhaplotypes. For vrs1.a2 and vrs1.a3, we could deduce a direct descent from Vrs1.b2 and Vrs1.b3, respectively, as a result of point mutation. Phylogenetic analysis ( Fig. 4 B) supported the hypothesis that the six-rowed alleles were derived from two-rowed alleles, rather than vice versa, because the wild barley lines (OUH602, OUH630, and OUH743) were outgroups. The vrs1.a2 allele has an insertion of one nucleotide in exon 2, which results in a frame shift of the deduced amino acid sequence ( Fig. 4 A). The vrs1.a3 allele has a substitution of one nucleotide in which the phenylalanine (F) at amino acid position 75 becomes leucine (L)( Fig. 4 A). The F is highly conserved in the DNA-binding domains of plants, animals, and yeasts, as shown by BLASTP analysis against the database of conserved domains ( 20). Three hex-v mutants ( hex-v.12, hex-v.13, and hex-v.14) also revealed the replacement of F by another nucleotide at position 75, supporting the importance of this amino acid in the DNA-binding domain. The vrs1.a3 allele occurred in one six-rowed cultivar ("Natsudaikon Mugi") in East Asia. The most common allele, vrs1.a1, which is distributed around the world (10 cultivars), has a deletion of one nucleotide in exon 3, which results in a different frame shift ( Fig. 4 A). The progenitor of the vrs1.a1 allele remains to be identified ( Fig. 4 B). It was noteworthy that the vrs1.a1 allele has been differentiated into three subhaplotypes by point mutations in the 5′ noncoding region ( Fig. 4 A). Two- and six-rowed barleys showed equal levels of transcription of Vrs1 ( Fig. 4 C), thus single-base-pair mutations in vrs1 appear to be responsible for the functional changes observed in all three natural six-rowed alleles.

Discussion

The evolution of barley was not a single event but has been an ongoing process of continuous specialization, as suggested by the concept of diffuse centers of differentiation ( 12). All of the present barley forms have resulted from specialization as a result of cumulative mutations, and tracing the origin of six-rowed barley would require us to find the least-specialized form of barley either conceptually (based on genetic analysis) or archaeologically. One of the main arguments in the debate over the evolution of barley relates to genetic changes between the brittle and nonbrittle rachis and between the two- and six-rowed spikes. We previously inferred the phylogeny of rachis brittleness ( 21) and of the number of rows of spikelets ( 22) based on molecular markers that are tightly linked to the btr1/btr2 and vrs1 genes. The genealogy based on each domestication gene (here, a domestication gene is one that motivates humans to domesticate a plant) allows adequate inferences about the origins of cultivated barley. In the present study, cloning of the Vrs1 gene substantiated our proposed phylogeny for the number of rows of spikelets. Wild and cultivated barley are interfertile ( 9, 23), and mutual introgression of genes appears to be substantial between the two forms ( 24). Therefore, random markers do not appear to represent the two- and six-rowed types of barley, and evolutionary studies based on DNA markers that are inherited independently of the domestication genes might provide different insights into the specialization of barley ( 25, 26). In the present paper, we focus exclusively on the origin of six-rowed barley. In a more strict sense, we focus on the Vrs1 gene rather than the whole genome of barley to elucidate the origin of six-rowed barley, even though our data might have some implications for evolution of the species.

An HD-ZIP I Homeobox Gene Determines the Number of Rows of Spikelets in Barley.

Our study revealed that Vrs1 encodes a member of the HD-ZIP class of transcription factors. Analysis of induced mutants clearly confirmed the identity of Vrs1 and demonstrated its biological function. Although complementation analysis by transformation is another option ( 27), our data, based on as many as 57 independent mutant lines, together with the specific gene expression patterns, confirm that Vrs1 encodes an HD-ZIP I protein. Although HD is universal, the HD-ZIP is unique to the plant kingdom. HD-ZIP protein forms dimers via ZIP and correctly binds by means of HD to dyad-symmetrical recognition sequences of DNA based on the strict spatial relationship between HD and ZIP ( 28, 29), and the DNA-binding activity of HD-ZIP was demonstrated ( 30). HD-ZIP proteins have been grouped into four families (I-IV) ( 31), and the barley VRS1 protein belongs to family I (data not shown). The HD-ZIP III and IV genes are functionally well characterized, being involved in ( i) development of the apical meristem, ( ii) vascular development, and ( iii) establishment of cell fates in epidermal ( 32). The biological function of the HD-ZIP I and II genes has not yet been clarified, although involvement in plant growth and development was suggested ( 32). Our study shows a conspicuous association between an HD-ZIP I gene and a plant developmental process. The spatial and temporal specificity of Vrs1 gene expression suggests that VRS1 is a transcription factor involved in the development (suppression) of lateral spikelets in two-rowed barley. Most domestication genes have been found to encode transcription factors ( 33), and our results agree with these observations.

During the domestication of barley and wheat, rudimentary spikelets or flowers of wild species have been restored to their full functionality to increase seed number ( 1, 8). In the present study, we demonstrated that loss of function of Vrs1 is sufficient to convert the rudimentary spikelets of two-rowed barley into the fertile spikelets of six-rowed barley. The mutations found in six-rowed cultivars were commonly represented by single-nucleotide polymorphisms, but even more drastic changes such as complete deletion of Vrs1 were also found in several induced mutants ( Figs. 3 and 4). This loss of function agrees with the facts that six-rowed spikes are genetically recessive and that two-rowed spikes are controlled by the wild-type allele, which occurs naturally in wild barley. The creation of six-rowed spikes is similar to the gigantism that occurs during domestication because both types of change can be caused by loss of function of genes ( 34). The dominant nature of Vrs1 and the potential DNA-binding activity of HD-ZIP I proteins suggest that VRS1 is a repressor protein that binds to the DNA of genes that regulate the development of lateral spikelets. Further investigation of the subcellular localization of VRS1 proteins will be necessary to test this hypothesis.

Multiple Origins for Six-Rowed Barley Implied by the Mutant Vrs1/HD-ZIP I Gene.

The number of rows of spikelets is a key characteristic in inferring the origin of cultivated barley, and it has been debated whether the progenitor of cultivated barley was six-rowed ( 35), two-rowed ( 12, 36), or both ( 37, 38). Recently, the two-rowed progenitor hypothesis was supported by ( i) archaeological specimens showing preagricultural utilization of wild barley and the existence of domesticated two-rowed remains that were older than six-rowed barley remains, ( ii) the dominance of the two-rowed phenotype, and ( iii) the rather weedy nature of brittle six-rowed barley in Israel and Tibet ( 9- 12).

Our phylogenetic tree ( Fig. 4 B) shows two lineages of domesticated two-rowed barley that diverged considerably from each other: Vrs1.b3 stood alone, and the other lineage further differentiated into Vrs1.b2 and Vrs1.t ( Fig. 4 B). Owing to its provenance in the western Mediterranean, the Vrs1.b2 allele is likely to be an immediate progenitor of the vrs1.a2 six-rowed allele, which predominates in the western Mediterranean. This is an indication that one of the six-rowed alleles originated in this region. Similarly, it is clear that the vrs1.a3 six-rowed allele derived directly from the Vrs1.b3 two-rowed allele by means of a point mutation ( Fig. 4 A). Indeed, the Korean vrs1.a3 cultivar has great similarities to old European two-rowed cultivars, such as a spring growth habit, nonbrittle rachis due to the recessive btr1 gene ( 38), lax spikes, a long stem, and the K haplotype of cMWG699, which is specific to two-rowed barley ( 22). It is interesting to speculate that the vrs1.a3 six-rowed cultivar might have originated from European two-rowed barley introduced into Eastern Asia. Our results suggest that, at least for the vrs1.a2 and vrs1.a3 alleles, six-rowed barley was derived from cultivated two-rowed barley. Because the number of lines of wild barley was limited in the present study, further collection of wild barley will be necessary to find a potential wild ancestor for each allele.

We were unable to resolve the origin of the most common vrs1.a1 six-rowed allele in two-rowed barley ( Fig. 4 B). A notable feature of the vrs1.a1 allele is that it has three subhaplotypes ( Fig. 4 A). If we assume the minimum number of changes, the vrs1.a1 allele must have originated from a hypothetical Vrs1.b1 allele in a single, unknown, two-rowed ancestor, and the first vrs1.a1 allele (which today exists in cv. 'Soren Oomugi 19329'; Fig. 4 A) has differentiated into two additional haplotypes by means of background point mutations in the 5′ noncoding region. The vrs1.a1 cultivars correspond to the major Type I cultivars classified by use of a marker tightly linked to vrs1, and the dominant Type I cultivars are distributed throughout the world ( 22). Because of its predominance and higher haplotype diversity, we hypothesize that the vrs1.a1 allele may represent the most ancient six-rowed allele (dating to 7,000-6,000 yBP), which was widely present in Neolithic agriculture in the Near East.

Our study demonstrates three independent processes for the origins of six-rowed barley ( Fig. 4 B). (This result does not necessarily mean that barley has a polyphyletic origin as a species, a distinction that is often misunderstood.) We previously hypothesized two independent origins of six-rowed barley, one of which gave rise to the six-rowed barley that has spread around the world and the other of which gave rise to the six-rowed barley endemic to the western Mediterranean ( 22). Our present study supports the concept of multiple origins of six-rowed barley. We believe these origins are straightforward because we used the DNA sequence encoding the Vrs1 gene itself. Creation of the three alleles ( vrs1.a1, vrs1.a2, and vrs1.a3) cannot be deduced through intragenic recombination by mutual hybridization of six-rowed barleys possessing different alleles ( Fig. 4 A), which suggests that their origins were independent. Divergence of these alleles cannot be explained by simple mutation. Thus, six-rowed barley must have originated repeatedly from cultivated two-rowed or wild barley by means of a loss-of-function mutation of an HD-ZIP I-class homeobox gene. Six-rowed barley originated independently at different times and in different regions ( Fig. 4 B), probably as a result of conscious selection by delighted farmers who greatly appreciated the improved seed yield.

Implications for the Evolution of Barley.

The DNA sequence of Vrs1 may trace the origin and migration of two-rowed barley. Both two- and six-rowed barleys emerged in Greece between 8,000 and 6,000 yBP, but six-rowed barley dominated in the Balkans and Central Europe (5,000-2,000 yBP) and in Southern Europe and North Africa (7,000-4,000 yBP) ( 1). After its disappearance from ancient Mesopotamia and Egypt, two-rowed barley does not reappear in the archaeological record in these regions until 1,100 yBP ( 11). Cultivation of two-rowed barley is virtually unknown in Central and Northern Europe until 1,000 yBP, and it has been assumed that two-rowed barley was introduced into Europe at 900-800 yBP by Crusaders from the Near East, as suggested by Fischbeck ( 39). In the present study, 14 two-rowed cultivars carrying the Vrs1.b3 allele did not show any polymorphism in the DNA sequence ( Fig. 4 A). This indicates that the origin of the Vrs1.b3 allele is relatively recent and that the allele has expanded rapidly throughout the world and become dominant, which partially supports Fischbeck's hypothesis ( 39). The interpretation that Vrs1.b3 is a recent allele agrees with the fact that only one vrs1.a3 six-rowed cultivar (derived from Vrs1.b3) has been detected thus far ( Fig. 4 A). This hypothesis does not necessarily contradict the hypothesis of Helbaek ( 11) that two-rowed barley was cultivated earlier than six-rowed barley, because the DNA sequences of archaeological specimens remain unknown. Moreover, the type of lateral spikelet in the ancient domesticated two-rowed cultivars ( deficiens, Vrs1.t; distichon, Vrs1.b) is not clear. We favor the hypothesis that the common progenitor of Vrs1.t and Vrs1.b2 occurred earlier than that of Vrs1.b3, based on the genetic distances seen in the allele coalescence and polymorphism results for the Ethiopian and western Mediterranean cultivars ( Fig. 4 A). Detailed genetic analysis of archaeological barley specimens would help to answer these questions.

The inflorescence architecture in the Poaceae is a continuous story of reduction from a more original "panicle" of spikelets (as seen in rice and oats) to a "spike" of spikelets ( 40), resulting in three sessile spikelets per node in barley and a single sessile spikelet per node in wheat and rye. In wild Hordeum species, the three spikelets and their slender awns form a light dispersal unit that permits both anemochory and zoochory ( 9). In two-rowed barley, strict temporal and spatial regulation of Vrs1 expression leads to reduction and sterility of the lateral spikelets. We speculate that either strong alleles or differential expression or regulation of Vrs1 orthologs could lead to complete repression of lateral spikelet formation at inflorescence nodes of other species in the Poaceae, as is found in wheat and rye. A Poaceae-wide assessment of the variability and regulation of Vrs1 orthologs would be an exciting and productive way to improve our understanding of plant development and the evolution of grass species.

Materials and Methods

Plant Materials.

Barley cultivars representing geographical distribution and three wild lines used in this study are maintained in the laboratory of T.K. (see SI Table 4). Mutants consisting of 39 hex-v and 16 Int-d lines were obtained from the Nordic GenBank (Alnarp, Sweden). Two more hex-v-type mutants derived from Japanese two-rowed cultivars New Golden ( 41) and Misato Golden ( 42) by gamma irradiation were obtained from T. Makino, National Institute of Crop Science (Tsukuba, Japan).

Genetic Mapping.

Three fine-mapping populations of 6,269 gametes were analyzed previously ( 15, 16). An additional 1,781 F ₂ plants of vrs1.a1 × Vrs1.b3 were screened with cMWG699 and BC12063, and recombinants were analyzed further (see SI Table 5). ESTs contigs BC12348 (AJ468022 and CB881790) and BC12063 (BU973565 and BU987455) were converted to sequence-tagged site markers. Other markers were generated from BAC DNA sequences.

BAC Contig and Sequencing.

A BAC library of cv. Morex ( 43) was obtained from Clemson University Genomic Institute (Clemson, SC) and has been screened by PCR ( 27). BAC contig was made by chromosome walking (see SI Methods). Shotgun sequencing of the BACs was carried out by the standard method ( 44). Assembled sequence contigs were correctly ordered, oriented with information of the cloning sites of the vector and the end-sequences of the bridge subclones, and annotated (see SI Methods). Contig gaps were filled by full-sequencing the bridge subclones either by primer walking or by GPS-1 transposon sequencing (New England Biolabs, Ipswich, MA).

RNA Analysis.

Immature inflorescences were developmentally staged by measuring under a microscope. RNAs were extracted from leaves or immature inflorescences excised from 2-month-old plants by using TRIzol (Invitrogen, Carlsbad, CA). First-strand cDNA was generated by using SuperScript II (Invitrogen). RT-PCR and full-length cDNA sequencing were carried out with gene-specific primers (see SI Methods).

In Situ Hybridization Analysis.

A 315-bp RT-PCR fragment including only the 3′-UTR of Vrs1 was amplified by using primers M111L15F84380 (5′-CAT ACT TAA CGC ACG CCT AGA GATC-3′) and M111L15R84670 (5′-TAG CTG CTG CCG CCG CCA AAT CCTC-3′) and cloned into pCR4-TOPO vector (Invitrogen); after digestion by EcoRI, the insert was subcloned into pBluescript II KS (+) vector (Stratagene, La Jolla, CA). Two clones with different insert orientations were linearized by using NotI and were used as templates to generate antisense and sense probes by using T3 RNA polymerase. Hybridization with a digoxigenin-labeled RNA probe, and immunological detection, were conducted. In situ hybridization was conducted according to the methods of Kouchi and Hata ( 45), with thorough modifications for improved reactions (see SI Methods).

Haplotype Analysis.

Genomic DNAs (2.1-kb) including Vrs1 were amplified from wild barley, cultivars, and mutants by PCR using the gene-specific primer pair M111L15F84204 (5′-GAA AGA TGA TTG CCA ACT ACC-3′) and M111L15R86329 (5′-GTC ATA ACT CGG CAA ACA TAG-3′), and their DNA sequences were determined by using internal primers designed by taking advantage of the BAC sequence (see SI Methods).

Phylogenetic Analysis.

A phylogenetic tree was constructed by the neighbor-joining method ( 46) using PAUP 4.0b10 ( 47). The DNA sequences of single haplotypes or subhaplotypes represented the cultivars sharing identical sequences. Phylogenetically informative insertions/deletions (indels) were included for analysis. All substitutions and indels were weighted equally. Bootstrap analysis with 1,000 replicates was performed.

Acknowledgments

We thank S. Okamura, M. Sameri, T. Ueda, T. Akihiro, N. Hara, S. Toki, and Y. Nagamura for help and advice and H. Kitano, M. Ashikari, D. Saisho, N. Minaka, D. Vaughan, and H. Knüpffer for discussions. This work was funded by research grants MP1113b (Rice Genome Project) and GD3006 (Green Techno Project) from the Ministry of Agriculture, Forestry, and Fisheries of Japan.

Footnotes

†To whom correspondence should be addressed at:
National Institute of Agrobiological Sciences, Kannondai, Tsukuba, Ibaraki 305-8602, Japan.
E-mail: takao{at}affrc.go.jp

Author contributions: T.K. and M.P. contributed equally to this work; T.K., M.P., A.G., M.M., and M.Y. designed research; T.K., M.P., C.H., P.A., H.K., D.P., N.S., A.T., U.L., T.F., and T.M. performed research; T.K., M.P., H.K., T.W., M.M., T.M., and M.Y. analyzed data; and T.K., M.P., and N.S. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS direct submission.

Data deposition: The sequences reported in this paper have been deposited in the EMBL/GenBank/DDBJ database (accession nos. AB259782, AB259783, and EF067844).

This article contains supporting information online at www.pnas.org/cgi/content/full/0608580104/DC1.

Abbreviations:

HD-ZIP,

homeodomain-leucine zipper;

yBP,

years before present.

Freely available online through the PNAS open access option.

View Abstract

# data[].attributes.attachment_link.title
Six-rowed barley originated from a mutation in a homeodomain-leucine zipper I-class homeobox gene

# data[].attributes.attachment_link.description
Increased seed production has been a common goal during the domestication of cereal crops, and early cultivators of barley ( Hordeum vulgare ssp. vulgare ) selected a phenotype with a six-rowed spike that stably produced three times the usual grain number. This improved yield established barley as a founder crop for the Near Eastern Neolithic civilization.

# data[].attributes.body
Global dissemination of a single mutation conferring white pericarp in rice

Please put a special emphasize on haplotype analysis

Presenters:Chris Byatt, and Salouf Almadani, Deirdre Sumski

# data[].attributes.attachment_link.readable_text
Megan T Sweeney,

Michael J Thomson,

Yong Gu Cho,

Yong Jin Park,

Scott H Williamson,

Carlos D Bustamante,

Susan R McCouch

Subject Areas

For more information about PLOS Subject Areas, click here.

We want your feedback. Do these Subject Areas make sense for this article? Click the target next to the incorrect Subject Area and let us know. Thanks for your help!

Abstract

Here we report that the change from the red seeds of wild rice to the white seeds of cultivated rice (Oryza sativa) resulted from the strong selective sweep of a single mutation, a frame-shift deletion within the Rc gene that is found in 97.9% of white rice varieties today. A second mutation, also within Rc, is present in less than 3% of white accessions surveyed. Haplotype analysis revealed that the predominant mutation originated in the japonica subspecies and crossed both geographic and sterility barriers to move into the indica subspecies. A little less than one Mb of japonica DNA hitchhiked with the rc allele into most indica varieties, suggesting that other linked domestication alleles may have been transferred from japonica to indica along with white pericarp color. Our finding provides evidence of active cultural exchange among ancient farmers over the course of rice domestication coupled with very strong, positive selection for a single white allele in both subspecies of O. sativa.

Author Summary

Understanding the history and origin of genetic mutations that have changed wild plants into crops can help us understand the history of the people who cultivated these plants. Rice is one of the oldest crops grown in Asia and it contains two different subspecies that are believed to have been domesticated in different locations by different people. Surprisingly, some of the genetic mutations responsible for domestication are common in all rice. We here show that a mutation in the Rc gene that changed the red seed of wild rice into the white seeds of modern rice is shared by a large majority of all rice varieties, regardless of subspecies. This transfer of genes requires contact among rice types and implies contact among the people who cultivated the different subspecies. We have traced the origin of the mutation in Rc to the japonica subspecies. As additional domestication genes are cloned and their evolutionary history described, we will see how many times and in how many directions such gene transfers have occurred.

Citation: Sweeney MT, Thomson MJ, Cho YG, Park YJ, Williamson SH, Bustamante CD, et al. (2007) Global Dissemination of a Single Mutation Conferring White Pericarp in Rice. PLoS Genet 3(8): e133. https://doi.org/10.1371/journal.pgen.0030133

Editor: Jonathan Pritchard, University of Chicago, United States of America

Received: May 4, 2007; Accepted: June 25, 2007; Published: August 10, 2007

Copyright: © 2007 Sweeney et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Funding: We gratefully acknowledge financial support from the Plant Genome Program of the National Science Foundation (Award Number 0606461 to SRM and CDB, 0110004 to SRM, 0319553 to SMc and CDB, and 0301919 to MJT), the Generation Challenge Program (GCP 4748/A001 to SRM), and Chungbuk National University/National Institute of Agricultural Biotechnology (NIAB-05-6-11-30-2).

Competing interests: The authors have declared that no competing interests exist.

Abbreviations: RID, rice insertion/deletion polymorphism; RM, rice microsatellite

Introduction

Human efforts to domesticate plant and animal species began thousands of years before recorded history, leaving us to guess at the methods that transformed wild species into agriculturally important crops and livestock. Rice, Oryza sativa, was domesticated in Asia and is now grown on every continent throughout the world, with the exception of Antarctica. The species consists of two major subspecies, indica and japonica, whose separate genetic identities were recognized in ancient times and are maintained by sterility barriers coupled with the inbreeding habit of O. sativa [ 1- 3]. The indica and japonica subspecies can be further subdivided into five genetically distinct subpopulations. Several estimates have been made for the time of divergence between indica and japonica based on intron sequence and retrotransposon insertions, and all of these calculations place the time to the recent common ancestor at more than 100,000 years ago [ 4- 6]. This divergence time is an order of magnitude larger than the oldest estimates for rice domestication. The data strongly support at least two independent domestications of O. sativa from predifferentiated pools of the O. rufipogon wild ancestor [ 3, 4, 7, 8].

With such evidence for independent domestications, we would expect to see different mutations conditioning domestication traits fixed within the different subspecies. We recently cloned the Rc gene, a bHLH protein, which is the only reported locus in rice that effects a change from red grain to white grain [ 9]. The change in pericarp color from red to white is an important hallmark of rice domestication. Most rice cultivars grown and consumed throughout the world today have a white or beige pericarp, the color of unpolished "brown rice" while all known accessions of the wild ancestor, O. rufipogon, invariably have grains with a red pericarp. Here, we investigate the number and origin of mutations that give rise to diverse white-grained rice cultivars from around the world, and the extent to which they were disseminated throughout the genetically and geographically diverse range of O. sativa. The evolutionary history of the Rc locus provides clear documentation that a single mutation moved rapidly from the japonica into the indica subspecies, while an independent mutation in the aus subpopulation was not widely disseminated during rice domestication.

Results

How Many Independent Mutations Are Responsible for White Pericarp in O. sativa?

To determine the frequency and distribution of the 14-bp deletion in the Rc gene that has been shown to cause white pericarp in rice [ 9, 10], a set of 440 geographically and genetically diverse rice cultivars were genotyped using the rice indel RID 12 primers reported in Sweeney et al. (2006). Our panel of varieties was selected to represent the range of diversity within Asian cultivated rice. It included both landrace and modern varieties from the five well-defined subpopulations of O. sativa [ 3] ( Figure 1) (204 indica, 33 aus, 87 temperate japonica, 99 tropical japonica, and 17 aromatic rices) collected in 24 different countries on three different continents ( Table S1). Molecular polymorphism data groups the indica and aus subpopulations within the indica subspecies and the tropical japonica, temperate japonica, and aromatic subpopulations within the japonica subspecies [ 2, 3].

For details about the geographical distribution of accessions used in this study, see Table S1.

One hundred percent of the varieties from the indica, tropical japonica, and temperate japonica populations ( n = 311) with white pericarp contained an identical 14-bp deletion in the bHLH gene; this deletion was found in none of those with red pericarp ( n = 103) ( Figure 2). This mutation was also found in 15 of the 17 white aromatic and four of the nine white aus varieties. Thus, a single 14-bp deletion in the Rc gene is present in both subspecies and in all five subpopulations and is responsible for white pericarp in 97.9% (330/337) of rice varieties surveyed.

The gene model for Rc, containing seven exons and comprising a region of approximately 8 kb, is shown horizontally along the top; arrows on the gene model and promoter region of Rc indicate positions of markers: RID-number markers amplify rice insertion/deletion polymorhisms, RM-number markers amplify rice microsatellite polymorphisms, and the RS-number marker detects a rice SNP. Physical positions along the rice Chromosome 7 pseudomolecule are indicated by the ruler below the gene model. In the tables below, allele sizes for RID marker loci are given in bp and allele sizes for RM markers show the number of repeats at each microsatellite locus. The number of varieties in each subpopulation that carry each haplotype is provided in boxed squares in the table to the right. Haplotypes found in varieties with red pericarp are shaded red.

To identify mutations that could explain the lack of pigment in the seven white rice varieties that did not contain the 14-bp deletion, we sequenced the entire coding region of Rc in two of the white aus varieties and one white aromatic variety that lacked the 14-bp deletion. As a control, nine red pericarp varieties from different subpopulations were also sequenced. A SNP (C→A) in exon 6 of the Rc gene distinguished these white- and red-grained rices. This SNP introduces a premature stop codon that truncates the protein before the bHLH domain. Further sequencing confirmed that this SNP was predictive of white pericarp in all seven white rices that did not contain the 14-bp deletion in Rc (DQ902346-DQ902352). Thus, the C→A SNP in exon 6 represents an independent mutation in the Rc gene that was not found in any of the indica, tropical, or temperate japonica varieties but exists at moderate frequency in white aus varieties (five of nine, 55.5%) and in two accessions out of 17 aromatic cultivars ( Figure 2).

Ancestral Haplotypes and the Origin of Rc Mutations

To determine the subpopulation origin of the Rc mutations, we examined ancestral haplotypes across the Rc coding sequence and promoter region in 103 genetically diverse, red-grained rices. Four rice insertion/deletion polymorphism (RID) and two rice microsatellite (RM) polymorphisms were used to construct haplotypes across the 6.5-kb region containing the Rc gene ( Figure 2). A single haplotype, H1, was detected among the 44 red-grained tropical and temperate japonica plants; we will refer to this as the ancestral japonica haplotype. No red rices belonging to the aromatic subpopulation were available for in this study. Five closely related haplotypes were identified among the 24 red aus and 35 red indica landraces, and all are clearly differentiated from the japonica haplotype. Haplotypes 6-7 and 9-11 ( Figure 2) are derivatives of each other and the similarity among them reflects the close evolutionary relationship between the indica and aus subpopulations [ 3]. When all red-grained indica, aus, and japonica accessions ( n = 103) are compared, it can be seen that they share the wild-type alleles at both functional markers (the non-deletion allele at RID12 and the C allele at SNP marker RS40, Figure 2). All flanking marker alleles were unique to one of the subspecies or subpopulation groups. Additionally, we collected ∼6,350 bp of sequence data for the Rc gene and ∼650 bp of downstream DNA from four red indicas and four red japonica varieties. Within this 7-kb region there are 21 SNPs whose alleles are subspecies specific. Thus, there are clearly differentiated haplotypes that distinguish the ancestral japonica gene pool (Haplotype Group A) from the ancestral indica and aus gene pools (Haplotype Group B) across this genic region. Based on these ancestral haplotypes, we were able to trace the origin of both functional mutations leading to white pericarp.

The most common haplotype of white pericarp rices, H2, contained the 14-bp deletion at RID12, and differed from the ancestral japonica haplotype, H1, only at this functional mutation. In contrast, the white H2 haplotype differed from the aus and indica red haplotypes (Haplotype Group B) at every marker tested ( Figure 2). This provides strong support for the conclusion that the H2 haplotype originated in a japonica ancestor. White haplotypes H3 and H4 are derived from H1, differing at a single marker locus in each case. Haplotype 5 provides evidence of a putative recombination event. It is found only in indica accessions and cultivars carrying this haplotype have ancestral indica alleles across the promoter and first intron of Rc and ancestral japonica alleles across the remaining sequence of the gene, including the 14-bp deletion. These data strongly support a single origin of the 14-bp deletion from a japonica background that then was introgressed into indica and aus.

Sequence data (7 kb across Rc, Table S2) from white accessions with the deletion (eight indica, two aus, eight tropical japonica, two temperate japonica, and one aromatic) and eight red accessions from both the indica and japonica subspecies were used to create a gene tree ( Figure 3). If the Rc gene tree and the tree created using genome-wide SSR data [ 3] had similar typologies, it would suggest an independent origin of white pericarp in indica and japonica. However, these trees are not consistent, suggesting introgression from one subspecies into the other. Accessions with red pericarp share a similar distribution pattern in both tress, but all whites with the deletion have japonica alleles at the 21 subspecies specific SNPs within the Rc gene and form a cluster with the red japonicas regardless of their subpopulation identity ( Figure 3). Thus, phylogenetic analysis confirms the japonica origin of the 14-bp deletion. Together, the lack of concordance between the Rc gene and rice genomic trees and the marker analysis provides strong support for the conclusion that the 14-bp deletion conferring white pericarp in rice arose once in the japonica gene pool and was widely introgressed into indica and aus landraces. The presence of this deletion within 97.9% of white-grained rice varieties found throughout the world today suggests either that the gene was dispersed during the early phases of domestication and is common by descent in modern varieties or a that very strong, positive selection for the allele lead to its introgression and maintenance in already established gene pools.

Sequence alignment adjusted manually in MacClade 4.05. A heuristic parsimony search was conducted on the data matrix, excluding gapped characters, with TBR branch swapping. This identified a single tree at 53 steps, with a consistency index of 0.98 (0.97 without autapomorphies) and a retention index of 0.99. Bootstrap support was estimated using 100 replications of the same search strategy. Pericarp color and subpopulation of each variety sequenced is labeled.

The second mutation Rc-s is only found in haplotype 8 and is restricted to six aus and two aromatic accessions. These white-grained varieties share no marker alleles with any of the tropical or temperate japonica accessions, suggesting that the origin of this allele is different from that of the 14-bp deletion. H8 differs from H7, a red aus haplotype, by only the functional C→A SNP. Since H7 is restricted to the aus subpopulation, we conclude that the C→A functional polymorphism originated in the aus subpopulation and that this allele was not widely disseminated during the domestication of O. sativa.

How Extensive Was Genetic Hitchhiking around the rc Allele?

How much DNA of japonica ancestry was introgressed into indica cultivars along with the 14-bp deletion conferring white pericarp is of particular interest because the Rc locus falls within the map positions of several quantitative trait loci associated with other domestication-related traits (dormancy, shattering, tillering, and panicle architecture) on Chromosome 7 [ 11- 17]. It is possible that an array of domestication alleles arose in the japonica background and were transferred as a block into indica cultivars. In order to define the extent of the introgression, we designed a series of indel and SNP markers that clearly distinguished the ancestral indica and japonica gene pools across Chromosome 7. To identify these markers, we evaluated polymorphism on a panel of 30 japonicas and 15 red-seeded indicas. Red indicas contain ancestral indica sequence in this region while both red and white japonicas (as well as white indicas) are expected to contain ancestral japonica sequence. F ST values provide a quantitative estimate of the degree of allelic differentiation between subpopulations. The indel and SNP markers used to distinguish the ancestral gene pools all had F ST values above 0.8 (see upper portion of Figure 4).

Arrows indicate the location of markers on the physical map of Chromosome 7, numbers are in Mb. F ST values for red indica versus japonica are plotted across the chromosome using diamonds, white indica versus japonica are plotted using squares. The expanded region shows a detailed view of the 1-Mb window around the Rc gene. All varieties genotyped are white indicas containing the 14-bp deletion. DNA segments containing alleles of japonica ancestry (in the indica background) are shown in white, segments containing alleles of indica ancestry are shown in black, and segments containing a breakpoint between alleles of japonica and indica ancestry are shown in grey. Locations of the Rc and Hsh genes indicated by brackets.

These markers were used to genotype 88 diverse white indica varieties to explore the extent of japonica DNA in the Rc region ( Table S3). In this study, nine extended haplotype patterns were identified and they showed varying sizes of japonica introgressions ( Figure 4, upper portion). Ninety-one percent of the indica varieties contained less than 1 Mb of japonica derived DNA in the region, with asymmetric distribution around rc. This is similar to haplotype patterns described for the Arabidopsis FRI gene [ 18] The haplotype with the smallest introgression, H2, which was observed in two varieties, was only 247-371 kb in size ( Figure 4, lower portion) and contained approximately 100 genes. Ten white indica varieties, H8 and H9, contained an introgression that included the recently mapped shattering locus, Hsh [ 14], although in most white indicas the japonica introgression does not extend that far. In an extreme case, H9, eight indica landraces collected from an isolated population in Kalimantan, Indonesia (on the island of Borneo) had japonica-derived alleles across the entire length of Chromosome 7 (29.7 Mb) ( Figure 4).

When F ST values were calculated using data from 88 white indicas and 30 japonicas, we observed a dramatic drop around Rc, from values of over 0.75, down to 0 and back above 0.75 within a 1-Mb region ( Figure 4, upper portion). F ST values around 0 indicate no differentiation between subpopulations and the drop around Rc illustrates the location and size of the japonica introgression ( Figure 4). The difference in F ST values outside the Rc region between the red or white indicas compared to japonica is due to the inclusion of the eight white indicas from Kalimantan that carry japonica alleles across the entire length of Chromosome 7. As more domestication genes are cloned from this region, the pattern and extent of introgression reported here can be used to determine if japonica alleles for other domestication traits hitchhiked along with the rc mutation when it was introduced into indica landrace varieties.

Genetic Diversity in Rc

To investigate the reduction in diversity around the 14-bp deletion in Rc, we used sequences from a portion of the Rc gene in 21 diverse varieties of white rice carrying the 14-bp deletion and eight red varieties representing indica, aus, tropical, and temperate japonica subpopulations ( Table S2). A summary of DNA polymorphism in the Rc gene is given in Table 1. Levels of polymorphism are reduced by 98% in the white pericarp rices, compared to the red landraces. Notably, the levels of DNA polymorphism in red landraces is comparable to randomly chosen loci in an unbiased sample of O. sativa landraces. For red landraces in the Rc region π (the level of nucleotide diversity) = 0.0025 per base pair, whereas in a survey of randomly chosen loci in O. sativa, π = 0.0032 (Ana Caicedo and Scott Williamson, personal communication, University of Massachusetts, Amherst, Massachussetts and Cornell University, respectively). Thus, other than the extreme reduction in diversity in white pericarp rices, the Rc region is not atypical of the rice genome. Furthermore, the site-frequency spectrum of the white rice sample is completely skewed towards polymorphisms with rare derived alleles. All of these observations are consistent with very strong, recent selection in favor of the 14-bp deletion in the Rc gene.

DNA Sequence Polymorphism in the Rc Gene

Discussion

Here, we show that the previously described 14-bp deletion in Rc leading to white pericarp in rice is found in over 97% of white pericarp rice varieties surveyed, including landrace representatives from all of the major rice subpopulations. This allele arose in a japonica background and was introgressed into the other subpopulations. We determined that the size of the japonica introgression surrounding rc in many indica varieties is less than 1 Mb.

A second mutation leading to white pericarp in 3% of rice varieties surveyed was identified and shown to have arisen in an aus background. This mutation is identical to the one described by [ 9] for the Rc-s allele conditioning light red or amber pericarp in the aus variety Surjamkuhi. As previously reported, the C→A transversion in exon 6 of the Rc gene results in a premature stop codon that truncates the protein before the bHLH domain. While it is not surprising that a truncated version of the bHLH protein would result in white pericarp, it is curious that the same truncated protein leads to light red pericarp in some genetic backgrounds. The explanation is expected to lie in the fact that bHLH proteins are found in complexes in maize, Arabidopsis, and petunia [ 19- 21], and different alleles of the interacting proteins are likely to determine the phenotype when Rc is truncated by the C→A substitution. A segregating population derived from a cross between white and light red varieties, both carrying the Rc-s allele, could pinpoint the genetic cause(s) of the difference in phenotype.

The lack of global dissemination of the Rc-s allele leading to either white or light-red pericarp may be partly a question of timing. If the Rc-s allele occurred after the rc allele (14-bp deletion) had become prevalent, there would have been no strong selective pressure driving its spread. The fact that in some genetic backgrounds the Rc-s allele produces a light red grain may be another reason why this allele was considered less desirable for early farmers, reducing the selective advantages and limiting its distribution.

How Was the rc Allele Disseminated across Asia?

In an outcrossing species, it is not surprising to see the rapid spread of a favorable mutation. However, O. sativa is 97%-98% inbreeding [ 1, 22], has a closed-flower morphology, and short-lived pollen grains. Gene flow is further constrained by the presence of a complex network of sterility barriers between the two independently domesticated indica and japonica subspecies of rice [ 1, 23]. Even among largely out-crossing wild Oryza relatives, viable pollen rarely travels more than 10 m [ 24], making it difficult for a favorable allele to be widely disseminated through natural pollen flow. Thus, it is unlikely that the rc mutation would have traveled far beyond its point of origin if not for the activities of humans who valued the white grains as a source of food and presumably as a commodity for trade. The importance of the rc mutation to early agriculturalists is evidenced by the fact that it moved around the Himalayan mountain range that is found between the proposed centers of indica and japonica domestication [ 4, 8], and, having traversed this substantial geographic barrier, was rapidly introgressed into all major subpopulations of rice despite an emerging fertility barrier.

Three key genetic features of white pericarp contributed to the rapid spread of this trait. First, rc is a single gene mutation that causes a qualitative change in phenotype so it is straightforward to visually distinguish red from white grains. However, the color of the seed coat is not visible on plants in the field because each grain of rice is completely covered by a hull, or glume. As dehulling reduces the germination rate, seeds that were to be planted the following season would be maintained with the hulls on. The pericarp color would only have been obvious to those who dehulled the grains, something that is normally done just prior to cooking. The pericarp is a maternally derived cell layer and its color is determined by the maternal genotype rather than the genotype of the fertilized embryo. Therefore, all seeds on any one plant are the same color. By dehulling only a few grains on any plant, it would be apparent which plants carried red or white seeds, allowing for selection. Third, the rc mutation is recessive, and in a highly self-pollinated species such as O. sativa, this means that the trait breeds true; seeds with white pericarp will produce offspring having grains with white pericarp for generations to come.

Why Did Farmers Prefer White Rice?

Almost all wild plants have pigmented seed coats, skin, and flesh, yet humans seem to prefer white starchy staples, and selection against pigmentation is a common theme in the evolutionary history of crop plants. The reasons our ancestors selected against pigmentation in these staple foods are not entirely clear, although it is tempting to speculate.

The novelty of a different color rice may have contributed to its popularity. White rice has the very practical advantage that it is far easier to detect and eliminate insects and pathogens against a light background than a darker red background. Alternatively, differences in cooking properties between red and white rice may have resulted in the preference for white. Colored rice has a harder seed coat than white rice. This means longer cooking times as well as more time spent finding fuel with which to cook. In many cultures, the rice hull and bran layers were removed by pounding. As red rice has a harder seed coat than white, it requires more work to remove the bran layer from reds. No matter what the reason for selecting against red pericarp in rice, selection was strong enough to catalyze a major change from red to white grains across a vast geographic area in a surprisingly short period of time.

Despite the lack of human records, the history of domestication is written in the genomes of the plants and animals that have been changed by our choices. As we begin to decipher these stories, we gain new insight into our own history as part of the domestication process. In rice, at least two important and unlinked domestication alleles, rc (discussed in this paper) and the shattering allele, sh4 (located on Chromosome 4) [ 25], are now known to have arisen only once in evolution and to have been introgressed across the indica-japonica divide, with almost complete prevalence in modern forms of cultivated rice. Other loci affecting domestication traits, such as the Rc-s allele for white pericarp (this paper) or the sh1 allele for non-shattering [ 26], confer alternative or complementary alleles for the same domestication phenotypes, but their distribution is restricted to a specific subpopulation. The fact that cultivars belonging to the deeply divergent indica and japonica subspecies have both common and distinct domestication genes suggests that the process of domestication in O. sativa involves complex patterns of subpopulation isolation and convergence, underwritten by a rich tapestry of cross-talk among ancient farmers. Whether a majority of key domestication alleles are subpopulation-specific or are shared among all Asian rice varieties, and whether the shared alleles arose predominantly within a single subspecies or subpopulation, is not yet known. The answer to this question will illuminate the paths traveled by ancient peoples and their innovations in the rice-growing world.

Materials and Methods

PCR and primer development.

List of primers in Table S4 All primers have an annealing temperature of 55 °C. For the introgression study, primers were designed to amplify small (4-30 bp) insertion/deletion events that were originally detected computationally by aligning the sequences from cv Nipponbare (japonica) and cv 9311 (indica) [ 27, 28]. As 9311 is a white indica, it carries a japonica introgression around rc. In order to identify indels between ancestral indica and japonicas in this region, short stretches of sequence from the red indica cv Mudgo were aligned with Nipponbare. To determine the pericarp color, ten seeds per plant were dehulled and then visually inspected.

Phylogenetic tree construction.

Rc DNA sequenced using overlapping primers. SNPs were confirmed with multiple reads through independent amplicons. Amplified products were sequenced at the Cornell Biotechnology Resource Center. Alignment adjusted manually in MacClade 4.05 ( http://macclade.org/macclade.html). A heuristic parsimony search was conducted on the data matrix, excluding gapped characters, with TBR branch swapping. This identified a single tree at 53 steps, with a consistency index of 0.98 (0.97 without autapomorphies) and a retention index of 0.99. Bootstrap support was estimated using 100 replications of the same search strategy.

Supporting Information

Accession Numbers

The sequences discussed in this paper were assigned National Center for Biotechnology Information (NCBI) GenBank ( http://www.ncbi.nlm.nih.gov/gquery/gquery.fcgi) accession numbers DQ885795-DQ885823 and DQ902346-DQ902352.

Acknowledgments

We are grateful to Gurdev Khush, Martha Hamblin, and Gael Pressoir for critical reading of the manuscript; Jeff Doyle for help constructing the phylogenetic trees in Figure 3; Gael Pressoir for suggesting the plotting of F ST values in Figure 4; Dr. Sutrisno and Fatimah Suwardjo of the Indonesian Center for Agricultural Biotechnology and Genetic Resources Research and Development in Bogor, Indonesia for providing DNA samples of the Indonesian rice varieties; the International Rice Germplasm Collection (IRGC) in the Philippines for access to diverse rice germplasm; Melissa Fitzgerald for helpful discussions concerning reason for rice preferences; and Lois Swales for help with formatting.

Author Contributions

MTS and SRM wrote the paper. SHW and MTS analyzed data. MTS, MJT, YGC, and YJP performed research. MTS, SRM, and CDB designed research.

References

1. 1. Oka HI (1988) Indica-Japonica differentiation of rice cultivars. Origin of cultivated rice. Amsterdam: Elsevier. pp. 141-179.

2. 8. Londo JP, Chiang Y, Hung K, Chiang T, Schaal BA (2006) Phylogeography of Asian wild rice, Oryza rufipogon, reveals multiple independent domestications of cultivated rice, Oryza sativa. Proc Natl Acad Sci U S A 103: 9578-9583.

3. 10. Furukawa T, Maekawa M, Oki T, Suda I, Iida S, et al. (2007) The Rc and Rd genes are involved in proanthocyanidin synthesis in rice pericarp. Plant J 49: 91-102.

4. 16. Thomson MJ, Tai TH, McClung AM, Lai XH, Hinga ME, et al. (2003) Mapping quantitative trait loci for yield, yield components and morphological traits in an advanced backcross population between Oryza rufipogon and the Oryza sativa cultivar Jefferson. Theor Appl Genet 107: 479-493.

5. 22. Matsuo T, K Hoshikawa (1993) Science of the rice plant. Tokyo: Food and Agriculture Policy Research Center. 686 p.

6. 27. Yu J, Hu S, Wang J, Wong GK, Li S, et al. (2002) A draft sequence of the rice genome (Oryza sativa L. ssp. indica). Science 296: 79-92.

7. 28. Goff SA, Ricke D, Lan TH, Presting G, Wang R, et al. (2002) A draft sequence of the rice genome (Oryza sativa L. ssp. japonica). Science 296: 92-100.

# data[].attributes.attachment_link.title
Global Dissemination of a Single Mutation Conferring White Pericarp in Rice

# data[].attributes.attachment_link.description
Author SummaryUnderstanding the history and origin of genetic mutations that have changed wild plants into crops can help us understand the history of the people who cultivated these plants. Rice is one of the oldest crops grown in Asia and it contains two different subspecies that are believed to have been domesticated in different locations by different people.

# data[].attributes.body
The Rc and Rd genes are involved in proanthocyanidin synthesis in rice pericarp

Please put a special emphasize on complementation analysis

Presenters:
Gillian Carney, Aoife McCooey, Claire McGee

# data[].attributes.attachment_link.description
Enter your email address below. If your address has been previously registered, you will receive an email with instructions on how to reset your password. If you don't receive an email, you should register as a new user

# data[].attributes.body
Caught red-handed: Rc encodes a basic helix-loop-helix protein conditioning red pericarp in rice

Please put a special emphasize on genetic mapping

Presenters:Annie Walshe, Stella Tattan, Robyn Healy

# data[].attributes.attachment_link.readable_text
Rough Mapping of QTLs and the Rc Mutant

Previous QTL mapping in this laboratory identified a single, significant QTL associated with red grain ( rg7.1) on chromosome 7 ( Figure 1B, i). This QTL was identified in two independent BC2 populations derived from crosses between an accession of O. rufipogon (IRGC-105491) from Malaysia and, in one case, a U.S. tropical japonica cultivar, Jefferson, and in the other case, a widely planted tropical indica cultivar, IR64. The log of the odds scores associated with the rg7.1 QTL peaks in these two populations were 99 and 33, respectively, and the QTL was detected in multiple environments ( Septiningsih et al., 2003). The peak of both QTLs corresponded to the previously mapped position of the mutant locus, brown pericarp, Rc ( Kinoshita, 1998). All of the BC2F1 plants had red seeds, indicating that the rg7.1 locus is dominant for red color, with the dominant allele donated by the O. rufipogon parent. Using the cv Jefferson/ O. rufipogon population, rg7.1 encompassed a 5.1-centimorgan (cM) region that represented ∼7.2 Mb straddling the border of the centromere on chromosome 7 ( Figure 1B, i). The genetic/physical distance in this region averages 1.4 Mb/cM, much above the genome average of 200 to 250 kb/cM, as expected for a pericentromeric region ( Zhao et al., 2002; Wu et al., 2003). An investigation of the genotype-phenotype relationship in 285 BC2F2 families demonstrated that all 18 families with red grain contained the O. rufipogon allele at either RM125 or the adjacent marker, RG30, suggesting that the gene underlying rg7.1 lay between these two markers.

Fine-Mapping of rg7.1

To fine-map the gene, 1410 BC2F3 plants were genotyped using markers flanking the QTL. The 72 recombinant plants were genotyped using six markers within the QTL region to locate the recombination break points more precisely, and the seed color of each recombinant line was recorded. Nineteen new simple-sequence repeat (SSR) and insertion/deletion (indel) markers were developed to help define break points across the region (see Supplemental Table 1 online). To narrow the region responsible for rg7.1 in each successive generation, we determined which segments of DNA from the red donor parent, O. rufipogon, were shared in all red-seeded progeny and eliminated from further consideration the O. rufipogon segments that appeared in white-seeded progeny. Because the trait is dominant, both the heterozygous and the homozygous O. rufipogon classes had identical phenotypes and therefore were grouped together during fine-mapping.

In the final BC2F6 generation, 4000 plants were genotyped and three classes of informative recombinants were identified, which narrowed the rg7.1 QTL to an 18.5-kb region ( Figure 1B, iii). Class 1 consisted of a single plant with a break point defined by markers RID13 and RM21197, located between the first and second exons of the gene LOC_Os07g11020.1, as illustrated in panel iv of Figure 1B. Upstream of RID13, these plants inherited O. rufipogon DNA, and downstream of RM21197, they were homozygous for cv Jefferson alleles, with the break point delimited to the region between the two markers. Because these plants had white seeds, we concluded that rg7.1 could not be located in the region upstream of RID13, which was heterozygous for O. rufipogon DNA. Recombinant class 2 consisted of 22 individuals having recombination break points between RID14 and RID15. This 14-kb region was highly repetitive, precluding the development of additional markers to help resolve the precise position of each break point. All plants in this class had red pericarp, and because they were all homozygous for cv Jefferson DNA downstream of RID15, we were able to eliminate that region as the location of rg7.1. Class 3 consisted of a single red-seeded recombinant plant that contained an intragenic break point between RM21197 and RID12. It had cv Jefferson DNA upstream of RM21197, allowing us to eliminate the region between RID13 and RM21197 from further consideration. We thus defined an 18.5-kb target region for rg7.1 bracketed by RM21197 and RID15 ( Figure 1B, iii).

Physical/Genetic Distance

The evaluation of 4000 plants for recombination allowed us to identify recombinational hot and cold spots within the 401-kb region. The most recombinagenic interval was a 6-kb region that contained no annotated gene models, defined by markers RID16 and RID17 ( Figure 2). Within this interval, the physical/genetic distance was 4 kb/cM, a recombination rate significantly higher than the average across the 401-kb region (P < 0.001). The 4-kb interval just downstream (RID17 to RID18) also lacked gene models but had no detectable recombination. The second most recombinagenic region (12.3 kb/cM; P < 0.001) was the 14-kb interval between RID14 and RID15 that contained two TEs. Again, this region was juxtaposed with a 2.3-kb region upstream in which there was no detectable recombination. The large interval between RM21177 and RM21194 (∼300 kb), with 50 annotated gene models only 3 of which are TEs, had a recombination rate significantly lower than the average for the region (P < 0.001), even after recombination rates from the two hot spots were excluded from the average ( Figure 2).

Positional and Functional Candidate Genes

One non-TE gene was detected within the 18.5-kb target region. LOC_Os07g11020.1 is a single-copy gene 668 amino acids in length and containing a predicted bHLH domain. This domain is common among transcription factors known to regulate pigment synthesis. As illustrated in Figure 1B, panel iv, a recombination break point between RM21197 and RID12 in recombinant class 3 had eliminated the promoter region and the first two exons of the bHLH gene as the source of the functional nucleotide polymorphism, leaving only exons 3 to 7 within the 18.5-kb target region. In addition to the bHLH protein, two putative transposon proteins, LOC_Os07g11040.1 and LOC_Os07g11030.1, both of the CACTA type, En/Spm subclass, were also present within the 18.5-kb target region. These were eliminated from further consideration based on three lines of evidence. First, the sequences of both of these proteins have >100 BLAST hits in the rice genome, with >98% identity over the entire length of the sequences. Mutations in either of these highly repeated genes would have no phenotypic consequences, because many other copies would remain functional in the genome. Second, domain analysis showed that LOC_Os07g11040.1 contains a transposase domain and LOC_Os07g11030.1 contains a proteinase domain, neither of which has been found in any of the regulatory or biosynthetic proteins responsible for anthocyanin/proanthocyanidin pigmentation in plants. Third, although TEs can be responsible for phenotypic changes if they are inserted within functional genes, neither of the two TEs within the rg7.1 QTL region shows any evidence of having disrupted any other genes. Therefore, based on our positional analysis and in silico functional interpretation, we postulated that the bHLH gene was responsible for red pericarp at the rg7.1 locus. We proceeded to test this hypothesis using three additional lines of evidence: sequence comparison of parental lines, sequence comparison of an allelic series, and expression analysis of the bHLH gene.

Sequence Comparison of the bHLH Gene

We sequenced the bHLH locus in both mapping parents, O. sativa (cv Jefferson) and O. rufipogon, to search for sequence changes that could explain the observed change in pericarp color. Having eliminated the promoter of the bHLH gene through recombination, we focused on changes that could affect the protein sequence. Six indels and 22 single-nucleotide polymorphisms were detected across the genomic sequence. We also compared the sequences of the mapping parents with the publicly available cv Nipponbare sequence. The cv Jefferson allele was identical to the Nipponbare sequence, both of which are japonica cultivars having white seeds.

Keeping in mind the possibility that annotation from the Nipponbare sequence might provide a gene model that differed from the dominant allele, we annotated the allele obtained from O. rufipogon and compared it with the gene model available at TIGR ( http://www.tigr.org/tdb/e2k1/osa1). In the Nipponbare annotation, the 3′ end of exon 5 and the 5′ end of exon 6 were truncated relative to the O. rufipogon gene model, predicting an mRNA in Nipponbare that was 513 bp shorter than the mRNA predicted from the O. rufipogon annotation ( Figure 3C). To confirm the accuracy of the different gene models, we amplified a segment of cDNA from both O. rufipogon and cv Jefferson within the only region in which the gene models differed. The cDNA amplicons were both 400 bp, the size expected from the O. rufipogon annotation ( Figure 3A). Sequencing of this amplicon from O. rufipogon confirmed the splice sites predicted from the O. rufipogon annotation. When the polymorphisms between cv Jefferson and O. rufipogon were aligned with the new gene model, 10 of the sequence polymorphisms fell within the coding sequence, and 5 of those are expected to affect the protein sequence ( Figure 4).

To help identify which of the sequence polymorphisms between the parents was responsible for the altered function of the gene, we also sequenced the bHLH locus in H75, an Rc mutant stock belonging to the japonica subspecies ( Figure 1A). H75, like O. rufipogon, carries a functional allele, but it is much more closely related to cv Jefferson than to O. rufipogon. Thus, a sequence comparison between H75 and cv Jefferson was expected to help eliminate some of the nonfunctional polymorphisms detected between the parents in the bHLH gene. We found that the coding sequence of the bHLH allele in H75 was identical to the cv Jefferson sequence except for a 14-bp indel in exon 6 ( Figure 4). This 14-bp sequence was present in the H75 stock as well as in O. rufipogon, but it was deleted in cv Jefferson and cv Nipponbare. The deletion induces a frame shift in the sequence, resulting in two premature stop codons before the end of exon 6. The stop codons truncate the protein before the bHLH domain. Given that this deletion was the only difference between the alleles of LOC_Os07g11020.1 in the H75 mutant stock (pigmented seeds) and the japonica cultivars cv Jefferson and cv Nipponbare (white seeds), that its location in exon 6 is consistent with the recombinational data, and that it would have a clear and important impact on gene function, we conclude that the 14-bp deletion is the only apparent reason for the lack of pigment in the pericarp of cv Jefferson and cv Nipponbare seeds.

To confirm that the bHLH protein underlying rg7.1 is also Rc, we analyzed a second mutant stock for sequence variation in the bHLH gene. The stock, Surjamkuhi, is an indica line that carries a third allele, Rc-s, conditioning light red seed pigmentation. This genetic stock offered independent confirmation of the identity of the Rc gene because different sequence polymorphisms in the same gene would be expected to distinguish the Rc-s, Rc, and rc alleles. The sequence of the bHLH gene in Surjamkuhi differed from the sequence of the japonica cultivars at many sites (as expected for varieties from different subspecies) but differed from the O. rufipogon allele at only four sites (positions 96, 660, 1353, and 1833 to 1844) ( Figure 4). The first two changes proved to be synonymous substitutions. The change at position 1353 consisted of a C-to-A change in exon 6. This single-nucleotide polymorphism was independent of any change seen in previous comparisons and represented a premature stop codon before the bHLH domain, truncating the protein and rendering the effect of the remaining indel immaterial. The fact that the different alleles of Rc show sequence polymorphisms that clearly account for the observed phenotypic differences is consistent with the conclusion that the bHLH protein is the Rc gene.

Expression Profiles of Rc and Biosynthetic Genes in White and Red Rice

To examine the timing and localization of the Rc transcript, we used RT-PCR to amplify mRNA from leaf, young panicle (before fertilization), pericarp of young seeds (at the milk or dough stage of grain filling), and pericarp from mature seeds. The mRNA was collected from both cv Jefferson (white seeds) and O. rufipogon (red seeds) plants. RT-PCR showed no expression of Rc in leaf tissue, as expected for a gene associated with a seed phenotype; however, expression was seen in several stages of panicle development ( Figure 3A). Because the promoter of the bHLH gene had been eliminated as the source of polymorphism based on the recombination data, we anticipated that similar expression levels of Rc would be detected in red and white seeds. Our results confirmed this expectation and further demonstrated that the RNA transcript from cv Jefferson contained the 14-bp deletion predicted from the sequence information ( Figure 3B).

Phylogenetic Comparison

To explore the evolutionary origin of Rc in rice and to identify putative orthologs in other species, we compared the sequence of Rc with other, previously identified bHLH transcription factors involved in anthocyanin and proanthocyanidin regulation as well as proteins with unknown effects that were recovered from BLAST using Rc as the query. The alignable portion of these sequences extended well beyond the bHLH domain (see Supplemental Figure 1 online), indicating that homology was not restricted to a single conserved functional domain. Analyzed using maximum parsimony or Bayesian analyses, these sequences fell into several clades ( Figure 5).

The divergence between sequences of different clades is substantial, making outgroup selection and the position of the root uncertain. Among clades 1 and 2, clades 4 and 5, and within clade 3, further alignment was possible, strengthening our findings that these groups of sequences are more closely related to each other. Therefore, it is likely that the root lies on one of the branches separating the three main groupings (clades 1 and 2, clade 3, and clades 4 and 5) from one another.

Several copies of this type of transcription factor appear to have been present in the ancestor of the monocots and eudicots, as clades 1 and 2 contain both monocot and eudicot sequences ( Figure 5). A third copy, present in the common ancestor of maize and rice, gave rise to clade 3, which shows gene duplication within each species. The paralogs within maize are known to confer tissue specificity of the anthocyanin pigmentation ( Goff et al., 1992). Clades 4 and 5 contain only eudicot sequences ( Figure 5).

It is clear from this analysis that Rc is not closely related to the rice bHLH proteins regulating anthocyanin, because they fall in different clades ( Hu et al., 1996; Sakamoto et al., 2001). These rice anthocyanin regulators are sister to the maize anthocyanin regulators, and they map to homologous locations on rice chromosome 4 and maize chromosomes 2 and 10, respectively. The phylogenetic analysis clusters Rc with In1 from maize. Rc and In1 are located in homologous chromosomal regions on chromosome 7 of both genomes. In1 is within maize bin 7.02, and of the 29 markers within 7.02 that map to the rice genome, 14 of them hit a 3-Mb region on rice chromosome 7, surrounding Rc. The intron-exon structures of In1 and Rc are similar, but not identical. Rc has seven exons, whereas In1 has nine, but in both genes the bHLH domain spans an intron-exon boundary. Although the phylogeny of Rc and In1 and the positions of the 14 homologous genes support orthology, In1 and Rc have different functions. In1 is a negative regulator of anthocyanin synthesis, which gives an intense purple color when mutated, whereas Rc is a positive regulator of proanthocyanidin synthesis, which has no color when the function is lost ( Burr et al., 1996). All of the other genes on the tree whose functions have been established are positive regulators, suggesting that the ancestral function is positive regulation ( Goff et al., 1992; Nesi et al., 2000; Spelt et al., 2000; Sakamoto et al., 2001; Bernhardt et al., 2003; Elomaa et al., 2003).

It is possible that two orthologous genes, Rc and In1, have evolved separate functions since the divergence of the common ancestor of rice and maize. To further test the hypothesis of orthology, we compared the nucleotide divergence (Ks) values between these genes and two other pairs of genes from bin 7.02 and the region around Rc (see Supplemental Figure 2 online). The Ks values for the Rc- In1 comparison are consistent with the expected values of orthologous gene pairs from these taxa. The surrounding genes had Ks values above and below those from the Rc- In1 comparison, showing that Rc and In1 were not changing faster or slower than other genes in the homologous region. Given the homologous positions of these loci and Ks values that are in agreement with expected values for two of the three genes analyzed, there is no evidence that contradicts the orthologous relationship between Rc and In1.

Thus, it appears that the rice and maize bHLH genes associated with pigment production in plants have evolved separate functions over time. The rice genes found in different clades have specialized, becoming part of either the proanthocyanidin or anthocyanin pathway, whereas in maize, the duplicated genes have become positive and negative regulators of the anthocyanin pathway.

Allelic Variation at Rc

The functional nucleotide polymorphism, a 14-bp deletion that knocked out gene function, was sufficient to explain the change in seed color between O. rufipogon and cv Jefferson. It remains to be seen whether independent mutations in this gene or other genes also give rise to the white phenotype in other lineages. We are currently undertaking an association study to determine the predictive power of these mutations for pericarp color in rice.

When comparing the rc and Rc-s alleles, it is not immediately apparent why the 14-bp deletion that frame shifts the bHLH domain should result in no pigment production and why a premature stop codon before the bHLH domain would give an allele conditioning light red color. In petunia and maize, an1 and b1 mutants lacking the bHLH domain are able to promote anthocyanin synthesis, much like the Rc-s allele ( Liu et al., 1998; Spelt et al., 2002). Insertions that cause a frame shift within the bHLH domain of an1 produce null alleles, a parallel to rc in rice ( Spelt et al., 2002).

Functional Parallels

In maize, petunia, and Arabidopsis, the bHLH proteins R/B, AN1, and TT8 interact with a Myb transcriptional activator ( C1, AN2, and TT2, respectively) to promote transcription of the anthocyanin and proanthocyanidin structural genes. In rice, only Rc has been reported as necessary for the production of colored seeds, although Rd is needed for the seed color to be red. The action of Rd is not consistent with that of the Myb transcription, as a functional copy of Rd is not required for the presence of any color in the pericarp. The Rd locus may be a biosynthetic gene that, when nonfunctional, results in the accumulation of a brown proanthocyanidin precursor. The rice gene OsC1, with similarity to the maize Myb gene C1, has been cloned on rice chromosome 6 ( Saitoh et al., 2004). This locus is acknowledged to have a role in purple (anthocyanin) rice pigmentation but has not been shown to play a role in red (proanthocyanidin) rice pigmentation.

It is possible that, unlike the anthocyanin regulatory systems in maize, petunia, and Arabidopsis, the bHLH gene encoded by the Rc locus is sufficient to activate the transcription of the rice structural anthocyanin genes alone. The bHLH genes do function as transcription factors in homodimers or heterodimers in animal systems, so there is no inherent reason that proteins containing this domain could not activate transcription on their own ( Heim et al., 2003). In maize, the PERICARP COLOR1 Myb transcription factor has evolved so that it does not require interaction with a bHLH protein to activate the transcription of genes involved in phlobaphene production ( Grotewold et al., 1994). An analogous change may have occurred in rice to allow the bHLH protein to function without an interacting partner. It is more likely, however, that rice also requires an interaction with a Myb transcription factor, but if both the red and white rices used in the genetic studies contain a functional copy of the Myb locus, it would not segregate and hence would not be detected as a contributor to pericarp color. If this is the case, a plant carrying a null allele at the Myb locus and a functional Rc should have white seeds.

Despite the differences in function between Rc and In1 from maize, the phylogeny of bHLH genes and Ks values support an orthologous relationship. It should be noted that the expected Ks for orthologous sequences (0.65) is also within the range expected of ancient polyploid events in these taxa ( Schlueter et al., 2004). Polyploidy could also have produced duplicated regions in the common ancestor of these taxa. If a paralogous region was lost in both the rice and maize lineages after divergence, then Rc rice and In1 maize could be ancient paralogs but still have the Ks values we found. Regardless of the possibility of ancient paralogy, In1 would still have diverged in function from the other positive regulators of anthocyanin and proanthocyanidin found in clade 1.

Phenotypic Associations

Red pericarp has long been used as marker for the cluster of domestication traits associated with weedy rice, including dormancy and shattering ( Gu et al., 2005). Several studies have placed QTLs for dormancy and shattering in the pericentromeric region of chromosome 7, encompassing the Rc locus. With the cloning of Rc, it is now possible to ask whether this association is the result of linkage or pleiotrophy. Given the reduced rate of recombination within the rg7.1 QTL, it is logical that genes for shattering, dormancy, and pericarp color have simply hitchhiked together in a linkage block. Indeed, fine-mapping of a rice shattering gene in this region has recently shown that this gene is tightly linked to Rc, although it occupies a different position on chromosome 7 (H. Ji, personal communication). It is also possible that Rc acts pleiotropically, as do many of the other bHLH proteins presented in the phylogeny ( Payne et al., 2000; Spelt et al., 2002; Bernhardt et al., 2003; Zhang et al., 2003). Using the recombinant lines generated in this work, we will be able to test these different hypotheses.

Rice and wheat ( Triticum aestivum) are similar in that red pericarp in both species can be eliminated by one locus. Comparative mapping shows no homology between the position of the Rc gene in rice and the R gene controlling red pericarp in wheat. A reverse genetics approach also failed to locate any ESTs from wheat that map to the Rc locus, although this is not surprising, because no rice or maize ESTs have been found for this locus either. Our work confirms that the Rc transcript only amplifies with a high sensitivity Taq polymerase, and this suggests that low transcript abundance may also explain the lack of EST hits in wheat. The R locus in wheat may be orthologous to the Rd gene in rice, given their homologous positions on wheat chromosome 3 and rice chromosome 1. Although the systems look similar phenotypically, molecular genetics analysis suggests that the mutations leading to white pericarp occurred at different points in the pathway.

Recombinational Analysis

The rg7.1 locus was originally mapped to a 7.2-Mb region that included the centromere on rice chromosome 7. Given the low frequency of recombination across this region, it was not clear whether positional cloning would be feasible. This study demonstrates that even in regions that are recombinationally repressed, map-based gene isolation offers a viable approach.

The Rc/rg7.1 locus is an area rich in gypsy retrotransposons and other repetitive elements ( www.gramene.org). Several previous studies have noted significant repression of recombination associated with the abundance of TEs and other repetitive sequences ( Fu et al., 2002; Wu et al., 2003; Shah and Hassan, 2005). Therefore, we expected that the TEs themselves might contribute to the recombinational repression. However, fine-mapping demonstrated that within the 401-kb region analyzed in the BC2F6 generation, one of the regions with a significantly increased recombination rate was TE-rich.

It has been shown that recombinational break points occur most often in regions with the greatest sequence similarity and that indels (including TEs) decreased the rate of crossovers more than single-nucleotide polymorphisms ( Puget et al., 2002). The TE components of varieties may be different, and the lack of similarity would be one explanation for the low recombination in TE-rich regions. A recent study in humans showed that the retrotransposons THE1A and THEA1B are overrepresented in recombinational hot spots and that the motif CCTCCCT, found associated with the retrotransposons, has a role in hot spot determination ( Myers et al., 2005). In our study, one of the TEs in a region with enhanced recombination near Rc contains a CCTCCCT motif. This offers an alternative explanation for the significant increase in recombination in this particular TE-rich region.

Implications

The cloning of Rc will make possible new methods of fighting weedy rice infestations in rice paddies. Red rice is a noxious weed that is currently responsible for losses of as much as $50 million per year in the United States ( Gealy et al., 2002). It is a perfect mimic of elite varieties, as the red pericarp is not visible until after harvest, when the grains are dehulled. Furthermore, the close association between red pericarp, seed shattering, and dormancy allows it to persist in fields for years despite vigorous attempts to remove it. The fact that the pericarp is maternal tissue, so that its color is dependent on the maternal genotype, means that seeds pollinated by red rice can be white (if the maternal plant carries the rc allele), but plants grown from these seeds will produce red seeds.

An immediate application of the work presented here involves the use of perfect markers that specifically target the 14-bp functional nucleotide polymorphism within the bHLH gene to screen for red rice contamination within certified seed lots. This will also facilitate the use of genes derived from crosses with wild relatives by allowing breeders to conclusively select against progeny carrying Rc, and to do so before the plants set seed.

QTL Population Development and QTL Detection

A BC2F2 population was constructed for QTL mapping using Oryza sativa subsp tropical japonica cv Jefferson as the recurrent parent and a wild accession of Oryza rufipogon (IRGC-105491 from Malaysia) as the donor parent ( Thomson et al., 2003). QTL detection was followed as described by Thomson et al. (2003).

Fine-Mapping of rg7.1

Eighteen BC2F2 families with red grain but with low levels of dormancy and shattering were used to fine-map the QTL. These families represent a subset of the red-grained phenotypes in the BC2F2 population. These segregating families were grown in 50-mm-wide × 178-mm-deep plastic pots in the Guterman Greenhouse at Cornell University. DNA was extracted using the Matrix Mill method ( Paris and Carter, 2000). The microsatellite and indel markers were amplified using standard PCR protocols, run on 4% polyacrylamide electrophoresis gels, and silver-stained as described ( Panaud et al., 1995). Seeds were harvested, and 10 to 15 seeds/plants were dehulled to determine seed color. Color was scored as either present or absent. In each generation, we determined which segments of DNA from the red donor parent, O. rufipogon, were present in all red-seeded progeny and discounted the O. rufipogon segments that appeared in white-seeded progeny as possible gene positions. Only plants heterozygous for the region containing rg7.1 were planted for the next generation of screening. A total of 5922 individuals from BC2F3 to BC2F6 were genotyped.

Molecular Marker Development

The required density of molecular markers in the target region was achieved using previously published SSRs ( McCouch et al., 2002) as well as SSR and indel markers developed as part of this study. New markers were designed from the publicly available rice genome sequence ( http://www.tigr.org/tigr-scripts/osa1_web/gbrowse/rice; http://143.48.220.116/genome_browser/index.html) using the SSRIT tool ( http://143.48.220.116/db/searches/ssrtool). Indel markers were designed by sequencing areas from both cv Jefferson and O. rufipogon and aligning the sequences using the SeqMan program of DNAstar (GeneCodes) to identify indels. Primer sequences, map positions, and amplified lengths of the newly developed SSRs and indel markers used in this study are listed in Supplemental Table 1 online.

Recombination Rate Statistics

The recombination rate for each interval was compared with the average using likelihood ratio test statistics and P values for the test that compares the two-parameter model in which the fragment has its own unique recombination rate (while everything else is uniform) with the model in which everything is uniform. The likelihood of each recombination rate = e λλxi/ xi! where xi = the number of recombinants in interval i, λ = r(length of sequence)(number of plants), and r is the per nucleotide recombination rate and is estimated as (total number of recombinants)/(total length of region)/(number of plants). The P values were corrected for multiple tests using a Bonferroni correction. The RM21177 to RM21194 interval was significantly lower than the average, even after the recombination rates from the two hot spots with P < 0.001 were excluded from the average.

Sequence Analysis of Red and White Cultivars

Overlapping primers between 700 and 900 bp were designed to cover the HLH gene and were used to amplify DNA with Pfu polymerase (Invitrogen). Amplified products were sequenced at the Cornell Biotechnology Resource Center.

RT-PCR

mRNA was collected from cv Jefferson and O. rufipogon panicles from the following stages and tissues: mature leaves, spikelets before pollination, dehulled seeds from the milk/dough stage of seed filling minus the starchy endosperm, and pericarp and seed coat scraped from the hardened endosperm of mature dried seeds. Total RNA was extracted with the RNAeasy miniplant RNA extraction kit (Qiagen). Total RNA (1 μg) was converted into cDNA with reverse transcriptase (Invitrogen), according to the manufacturer's instructions. As a control, water was used in place of reverse transcriptase for the reaction. The reaction was diluted threefold, and 2 μL of cDNA was used for amplification. A total of 10 μL of PCR products was loaded on 1% agarose gels. Primers for RT-PCR were designed to span introns. BLAST searches were used to ensure that primers were specific to the candidate gene of interest. The primers used were as follows: actin F, 5′-CGTCCTCTCTCTGTATGCCAG-3′; actin R, 5′-CTGGTACCCTCATCAGGCAT-3′; Rc F, 5′-GCCTTGTCACTCTTGGCATT-3′; Rc R, 5′-GGTTGGCACTGAAATCACCT-3′; Rc short F, 5′-CAGGCACCACACAGAGAATG-3′; Rc short R, 5′-CTCCTCTCTTTCAGCACATGG-3′.

Sequence Alignment and Tree Building

DNA sequences were aligned in BioEdit ( Hall, 1999) with the assistance of ClustalW ( Thompson et al., 1994) and bl2seq ( Tatusova and Madden, 1999) to identify similar regions, with extensive manual adjustment. Alignment was undertaken to initially preserve the reading frame of each sequence. Some gaps were later introduced that disrupted small sections of the reading frame to minimize the total number of events (substitutions, insertions, or deletions) required to transform one sequence to another. Only some regions of the sequences could be confidently aligned across all taxa, including the bHLH domain (∼1 kb). Further alignment was accomplished among sequences within the three main clades found in our analyses, with gaps for the remaining sequences (∼500 bp). Gaps were treated as missing in all analyses. Indels were not coded as additional characters.

Analyses

The aligned matrix was examined for the presence of secondary signal(s) using reverse successive weighting, implemented in RSW1.1a ( Trueman, 1998). Using a parsimony bootstrap cutoff of 70% and 1000 fast bootstrap replicates, no evidence of secondary signal was found. Further analyses were conducted with equal character weighting using parsimony implemented in PAUP* ( Swofford, 1998) as well as Bayesian analysis using HKY and GTR+G models implemented in Mr Bayes ( Huelsenbeck and Ronquist, 2001) version 3.1.1 to examine the sensitivity of the resulting phylogenetic inference to some model and method variation. Of 1449 aligned positions, 810 (55%) were parsimony-informative. A single most parsimonious tree of 3629 steps was found after a heuristic search using 100 random addition sequence replicates holding a maximum of 100,000 trees (maximum not reached); Bayesian analyses for both models used 10 chains and were run for 5 million generations, saving a tree every 1000 generations. Two replicate runs were compared with check convergence using the likelihood versus generation plot. Full mixing was evident from almost the first few generations, so an arbitrary 100 burnin trees were discarded before producing the consensus topologies and clade posterior probabilities. The parsimony tree differed slightly from the Bayesian majority rule tree ( Figure 5). The parsimony-preferred positions of taxa are shown by arrows where they connect to other branches.

Synonymous changes per synonymous site (Ks) were calculated for Rc rice and In1 maize ( Zea mays) sequences along with sequences from two flanking genes from both species, No Apical Meristem putative protein and putative cellulose synthase sequences. We used the method of Yang and Nielsen (2000) implemented in PAML ( Yang, 2000; Yang and Nielsen, 2000) to test whether these values were consistent with that expected for orthologous sequences, around Ks = 0.65, assuming a synonymous substitution rate of 0.65 × 10 −9 ( Gaut et al., 1996) and the age of divergence of rice and maize at 50 million years ago ( Gaut, 2002). Two approaches were considered when deciding how different Ks values could be from the expected values and not be significantly different. In the first approach, we compared two standard errors above and below the estimated Ks for each gene with the expected value of 0.65. This represents an approximation of the 95% confidence interval around each estimate. The other approach we used was to treat the expected value as the mean of a distribution of Ks values, and therefore an estimate for any one pair of genes. Extending a similar approach that has been used in comparisons of genes duplicated by polyploidy ( Zhang et al., 2003), we might expect the variation for 90% of genes duplicated by the same event (here speciation) to show 2.6-fold variation in Ks. Therefore, a range of 0.361 to 0.939 around 0.65 can be expected for most genes.

Accession Numbers

Sequence data for the genomic sequence of Rc can be found in the GenBank/EMBL data libraries under the following accession numbers: H75, DQ204735; cv Jefferson, DQ204736; O. rufipogon, DQ204737; Surjamkuhi, DQ204738. The cDNA fragment amplified from O. rufipogon has accession number DQ315482.

Supplemental Data

The following materials are available in the online version of this article.

Supplemental Table 1. Primers Developed for Fine-Mapping.

Supplemental Figure 1. Sequence Alignments Used in Phylogenetic Analysis.

Supplemental Figure 2. Ks Ranges for Rc-In1 Surrounding Gene Pairs.

View Abstract

# data[].attributes.attachment_link.title
Caught Red-Handed: Rc Encodes a Basic Helix-Loop-Helix Protein Conditioning Red Pericarp in Rice

# data[].attributes.attachment_link.description
Rc is a domestication-related gene required for red pericarp in rice ( Oryza sativa ). The red grain color is ubiquitous among the wild ancestors of O. sativa , in which it is closely associated with seed shattering and dormancy.

# data[].attributes.body
The molecular genetics of crop domestication
Presenters: David Prendergast, Niamh Tuohy, Jemma Wilson

# data[].attributes.attachment_link.title
The Molecular Genetics of Crop Domestication

# data[].attributes.attachment_link.description
Ten thousand years ago human societies around the globe began to transition from hunting and gathering to agriculture. By 4000 years ago, ancient peop...