Analysis of the bread wheat genome using whole-genome shotgun sequencing.

Bread wheat (Triticum aestivum) is a globally important crop, accounting for 20 per cent of the calories consumed by humans. Major efforts are underway worldwide to increase wheat production by extending genetic diversity and analysing key traits, and genomic resources can accelerate progress. But so far the very large size and polyploid complexity of the bread wheat genome have been substantial barriers to genome analysis. Here we report the sequencing of its large, 17-gigabase-pair, hexaploid genome using 454 pyrosequencing, and comparison of this with the sequences of diploid ancestral and progenitor genomes. We identified between 94,000 and 96,000 genes, and assigned two-thirds to the three component genomes (A, B and D) of hexaploid wheat. High-resolution synteny maps identified many small disruptions to conserved gene order. We show that the hexaploid genome is highly dynamic, with significant loss of gene family members on polyploidization and domestication, and an abundance of gene fragments. Several classes of genes involved in energy harvesting, metabolism and growth are among expanded gene families that could be associated with crop productivity. Our analyses, coupled with the identification of extensive genetic variation, provide a resource for accelerating gene discovery and improving this major crop.

Introduction

With a global output of 681million tonnes in 2011[1], bread wheat (Triticum aestivum) accounts for 20% of human consumption of calories[2] and is an important source of protein, vitamins and minerals. It originated from hybridization between cultivated tetraploid emmer wheat (AABB, T. dicoccoides) and diploid goat grass (DD, Aegilops tauschii) approximately 8,000 years ago[3]. Bread wheat cultivation and domestication has been directly associated with the spread of agriculture and settled societies such that it is now one of the most widely cultivated crops due to its high yields and nutritional and processing qualities. The three diploid progenitor genomes, AA from T. urartu, BB from an unknown, possible Sitopsis section species (to which Ae. speltoides belongs), and DD from Ae. tauschii, radiated from a common Triticeae ancestor between 2.5-4.5 MYA (million years ago) and AABB tetraploids arose less than 0.5 MYA[4,5]. Nucleotide diversity in the AABB and DD genomes is substantially reduced compared to ancestral populations, indicating a major diversity bottleneck upon the transition to cultivated lines[6].

Grass genomes exhibit extensive long-range conservation of gene order[7-9]. Nevertheless, they are highly dynamic due to the activities of repeats contributing to tremendous variation in genome size[10], changes in local gene order, and to pseudogene formation, particularly in larger genomes such as maize[11] and wheat[12]. From analysis of BAC contigs on chromosome 3B, the 17 Gb genome was estimated to be composed of approximately 80% repeats, primarily retroelements, with a gene density between 1/87-1/184kb[13]. Despite both the substantial knowledge gained of the wheat genome from these studies and the central importance of the wheat crop, a comprehensive genome-wide analysis of gene content has yet to be conducted due to its large size, repeat content and polyploid complexity.

We have analysed low-coverage, long-read (454) shotgun sequence of the hexaploid wheat genome using gene sequences from diverse grasses. From this we created assemblies of wheat genes in an orthologous gene family framework, utilized diploid wheat relatives to classify homoeologous relationships, and defined a genome wide catalog of single nucleotide polymorphisms (SNPs) in the A, B and D genomes. These analyses provide a new foundation for genetic and genomic analysis of this key crop.

Sequence analysis

The wheat variety Chinese Spring (CS42) was selected for sequencing because of its wide use in genome studies[14,15]. Purified nuclear DNA was sequenced using Roche 454 pyrosequencing technology (GS FLX Titanium and GS FLX+ platforms) to generate 85 Gb of sequence (220m reads), approximately 5× coverage based on an estimated genome size of 17 Gb. Supplementary Table 1 shows 79% of the reads had matches to the TREP repeat database, and most hit retrotransposons, consistent with previous studies[13]. To identify A-, B- and D-genome derived gene assemblies in the hexaploid sequences, we utilized Illumina sequence assemblies of T. monococcum, related to the A genome donor, Ae. speltoidies cDNA assemblies, and 454 sequences from the D genome donor Aegilops tauschii respectively. The SOLiD platform was used to generate additional sequence of CS42 and three commercial wheat varieties to increase the accuracy of homoeologous SNP identification. Datasets are summarized in Table 1a and Supplementary Table 2, and SNP identification methods are described in Supplementary Section 5.2.

Table 1a

Summary of sequence sources used for analysis.

Genome	Platform	Size of dataset	Reference
T. aestivum (CS42) genomicDNA	454 GS FLXTitanium/454GS FLX+	85Gb	EBI Study:ERP000319
T. aestivum(CS42) genomic DNA fromsorted chr 1A, 1B, 1D	454 GS FLXTitanium	1A: 287Mb1B: 392Mb1D: 375Mb	Wicker et al.,2011[20]
T. aestivum (CS42, Avalon,Rialto, Savannah) genomicDNA	SOLiD3/SOLiD 4	15.2 billion reads	EBI Study:ERP001493
T. aestivum (CS42) cDNA	454 GS FLXTitanium/454GS FLX+	1.6Gb	EBI Study:ERP001415
T. monococcum genomicDNA	IlluminaGAIIx/HiSeq	A/B/D sequences: 3.7Gb A/B/D SNPs: 401Gb	NCBI archive:SRP004490.3
Ae. speltoides cDNA	[Pre-assembleddata]	151Mb	(M. Trick & I.Bancroft,unpublished)
Ae. tauschii genomic DNA	454 GS FLXTitanium	12.8Gb	Luo et al., 2012,submitted
Ae. tauschii genomic DNA	SOLiD 4	80-100x	(J. Dvorak,unpublished)

Sequence assembly

An Orthologous group Assembly (OA) (Supplementary Table 3) was created by clustering 454 reads by sequence similarity to orthologous grass gene sequences, and separate assembly of the clusters at high stringency using Newbler (Supplementary section 2). The orthologous genes were derived from rice[16], sorghum[8], Brachypodium,[9] and barley full-length cDNAs, by OrthoMCL[17] clustering. This generated 20,496 Orthologous Groups (OG) (Supplementary Table 4 and Supplementary Figure 1). The gene model with highest similarity to wheat (termed the “OG Representative”) was selected from each OG by stringent BLASTX comparison to a Low Copy-number Genome assembly (LCG) made by filtering out repetitive sequences and assembling the remaining low-copy sequences de novo (Supplementary Table 3). Nearly 90% of the metabolic genes in Arabidopsis matched OG Representatives. The 20,051 OG Representatives matched 92% of publicly available wheat full-length cDNAs[18] and 78.7% of harvESTs (Supplementary Figure 2), indicating they represent nearly all wheat genes.

We optimized parameters for wheat gene assembly using MetaSim[19] to generate simulated 5× 454 reads from the allotetraploid maize genome and from a triplicated rice gene set, with the introduction of sequence variation (Supplementary Section 2.7). Similar degrees of coverage over the OG Representatives were seen for the simulated datasets and wheat 454 reads (Figures 1A). Rice reads followed the same depth distribution as the wheat reads (Figure 1B), suggesting they are a reasonable representation of hexaploid sequences. Maize reads covered their OG Representatives to a median depth of approximately 5, consistent with 5× coverage.

Figure 1

Coverage of OG Representatives by wheat 454 sequence reads and simulated 454 reads from rice and maize

a. Coverage of OG Representatives by repeat-masked wheat 454 sequence reads (black line), wheat LCG (black dashed line), the OA (blue line), together with rice genes (red line) and maize simulated reads (green line).

b. Median coverage depth over protein coding regions of OG Representatives (N terminus = 0; C = 100). The colour coding is the same as in panel 1A, except simulated hexaploid reads from rice (red line) were used.

c. The distribution of wheat gene copy numbers from the OA.

Simulated maize and triplicated rice 454 reads were used to optimize assembly parameters. Assembly at 99% minimum sequence identity (mi) using 40bp overlap length predicted gene family sizes most accurately (Supplementary Figures 3-6). Wheat 454 reads were pre-processed (Supplementary Table 5) and assembled using 99% mi (Supplementary Tables 6 and 7) to create the OA. Figure 1B shows that the depth of coverage of the OA assembly followed a similar pattern to maize, consistent with multiple gene copies. In contrast, the low depth coverage by the LCG assembly suggested gene family numbers were collapsed. The number of wheat assemblies for each OG Representative was calculated to determine gene copy numbers (Supplementary Table 7). Figure 1C shows that most OG Representatives had between 1-5 distinctive wheat gene assemblies, with a peak of 2 genes.

The A, B and Ae. tauschii (D) genome gene sets[13,20,21] have been estimated to be 28,000, 38,000 and 36,000 respectively. We estimated the number of genes in the hexaploid wheat genome to range between 94,000-96,000 (Supplementary Section 2.10). This is reasonably consistent with estimates based on wheat chromosome sequences[13]. Comparing our transcriptome assembly (Supplementary Sections 2.8 and 2.9) and wheat harvEST to the wheat OG Representatives showed that 76% and 65% respectively were expressed under the conditions used for RNA isolation. Similar results were found in barley[22], rice[16] and maize[23], indicating the assemblies are bona fide wheat genes.

We defined the overall extent of gene conservation between wheat and the most closely related sequenced pooid grass Brachypodium distachyon[9,24]. Figure 2 track 1 shows that there is a high degree of overlap between the gene sets of Brachypodium and wheat, but with regions of lower conservation, for example on Brachypodium chromosomes 1 and 4. Syntenic maps of the Brachypodium genome and the A, B and D chromosome groups were created by integrating high-density wheat EST-based markers[25] with Brachypodium genes (Figure 2 tracks 5, 6 and 7 respectively). Supplementary Figure 7 shows the A, B and D genome markers separately. Syntenic alignments were readily identifiable and conformed to the predicted major patterns[9,26]. We identified many insertions and/or translocations of blocks of genes within the overall conserved patterns of gene order, including the major rearrangement on chromosome 4A as shown on Brachypodium chromosome 1[20]. Lower marker density on the D genome is evident in track 7. The higher resolution genetic map identified a new syntenic alignment of Triticeae group 5 to Brachypodium chromosome 3 genes.

Figure 2

Alignment of wheat 454 reads, SNPs and genetic maps to the Brachypodium distachyon genome

The inner circle represent gene order on the 5 Brachypodium chromosomes. Track 1 illustrates conservation between wheat 454 reads and Brachypodium genes, shown as a window of genes present in wheat. Tracks 2-4 show SNP density (the mean number of SNPs per gene in a window of 20 genes) in the A (track 2), B (track 3) and D (track 4) genomes of wheat. Tracks 5-7 display wheat synteny with Brachypodium for the A (track 5) B (track 6) and D (track 7) genomes. Genetic markers[25](shown in darker colours) were colour-coded by wheat chromosome. Gaps between markers were filled in to show synteny (lighter colours).

Genome change in polyploid wheat

We determined the influence of polyploidy on gene content in hexaploid wheat by defining the sizes of gene families in hexaploid wheat and the diploid progenitor Ae.tauschii from the copy number of genes for each OG Representative, which were then paired with the gene family size of the OG Representative in sequenced diploid grasses (Supplementary Section 2.6). The mean family size was 1.4 members. Supplementary Figure 8 shows relationships between wheat and diploid orthologous gene family across the full scale of orthologous gene family sizes. This approach accurately reconstructed gene family sizes in simulated maize and “hexaploid” rice genomes (Figures 3A and 3B), although larger gene family sizes tended to be under-estimated. Figures 3C and 3D show the relationships between Ae. tauschii and wheat genes. Single member gene families in hexaploid wheat and Ae. tauschii were maintained to a similar extent as those seen in sequenced diploid grasses, consistent with Southern blot analyses of single copy genes[27]. Using the D genome as a diploid reference, we calculated the ratio of gene family sizes in hexaploid:diploid Triticeae as between 2.5-2.7:1, derived from the geometric mean and the slopes of the blue and red lines in Figure 3E respectively. Comparing the expected ratio of 3:1 for the hexaploid indicated the loss of between 10,000-16,000 genes in hexaploid wheat compared to the three diploid progenitors (Supplementary section 2.10). This is consistent with earlier studies of gene loss in newly synthesized wheat polyploids [28] and the erosion of genetic diversity during wheat domestication[6].

Figure 3

Gene family sizes in orthologous assemblies of hexaploid wheat, Ae. tauschii, simulated maize and hexaploid rice

The boxes and whiskers contain 50% and 90% of the OA genes respectively. The box colours indicate the number of genes in diploid gene families of different sizes (x axis). The black lines represent expected gene family sizes; the red line is the gene family size determined from the OA, derived by polynomial regression fit. Only gene families with up to ten members are shown.

a. Maize gene family sizes predicted from orthologous assembly of simulated 454 reads.

b. Rice gene family sizes predicted from orthologous assembly of simulated 454 reads derived from triplicated rice genes.

c. Ae. tauschii gene family sizes obtained from orthologous assembly of repeat-masked 454 reads. Expanded gene families are shown as green dots.

d. Wheat gene family sizes in the OA.

e. Amalgamation of wheat and Ae.tauschii gene copy numbers. The black line shows the expected gene copy numbers for wheat and Ae. tauschii respectively. The red line shows the regression fit for wheat, and the blue line for Ae.tauschii. The grey zone between these lines estimates the extent of gene loss in hexaploid wheat.

Despite this overall trend of gene family size reduction, gene families with fewer or greater than expected members were identified in Ae. tauschii and hexaploid wheat, shown as green dots (expanded) and brown dots (contracted) in Figures 3C and 3D, respectively. Supplementary Tables 10, 11 and 12 show the over- and under-represented functional categories of proteins. Most of the over-represented categories in expanded gene families are common to wheat and Ae.tauschii: these include ribosome proteins, components of photosystem II, storage proteins, transposon-related proteins, cytochrome P450s, NB-ARC domain proteins involved in defense responses, pollen allergen-related proteins, and F-box proteins. Five of the 11 families encoding hydrogen ion transmembrane transporters were significantly more numerous in Ae. tauschii compared to wheat. Analysis of gene families (Supplementary Figure 9) showed they encode different subunits of ATPases. We speculate that they may provide proton gradients to support Na+ exclusion in Ae. tauschii[29] and the accumulation of minerals in other Aegilops species[30].

Pseudogene analysis

Several classes of plant DNA transposons[31,32] and retroelements[33] create and amplify gene fragments, disrupt genes and create pseudogenes, which can influence gene expression through epigenetic mechanisms[34]. We identified a set of almost 233,000 gene fragments that mapped to the same regions of their OG Representatives, forming “stacks” that were sufficiently divergent not to assemble into their cognate gene assemblies (Figure 4A). Two classes were identified; those containing Pfam domains, and those aligning with non-Pfam domains of OG Representatives. Nearly 30% of the OG Representatives had associated gene fragments (Supplementary Table 13) that most frequently covered between 10-12% of their length (Figure 4B). Figure 4C shows that the alignment identities of gene fragments against their OG Representatives were substantially lower than the identities for cognate regions within wheat gene assemblies. Supplementary Figure 10 shows the distribution of “stacks” along genes and Ka/Ks analyses. Pfam domains found in “stacks” were enriched for Zn-finger motifs in Mutator transposons (Supplementary Table 14), consistent with their role in pseudogene formation[31]. F-box, protein kinase and NB-ARC domains, found in the most rapidly evolving gene families in plants[9,35], are also over-represented.

Figure 4

Pseudogene identification and analysis

a. Visualization of an OG Representative and associated wheat sequences. The top track shows the hit count profile of mapped 454 reads. The lower tracks shows sub-assemblies of three wheat genes and a stacked region of gene fragments. Read depth is represented by the heat map.

b. The distribution shows coverage of the OG Representative by Pfam-containing gene fragments and pseudogenes. The blue and red lines represent stacks with and without protein domains, respectively.

c. The distribution shows protein identity between sub-assemblies forming stacks of gene fragments. The blue and red lines represent stacks with and without protein domains, respectively, and the black line represents sub-assemblies forming genes.

Determining homoeologous relationships of gene assemblies

We classified gene assemblies as A-, B- or D-derived according to sequence similarity to Illumina sequence assemblies from T. monococcum, cDNA assemblies from Ae. Speltoides, and 454 reads from Ae. tauschii respectively by applying a Support Vector Machine-learning approach (SVM, Supplementary Section 5, Supplementary Figures 11 and 12, Supplementary Tables 15-18). Supplementary Figure 13 shows 66% of the gene assemblies were classified with high overall precision (>70%) and recall into 28.3% A genome, 29.2% B genome and 33.8% D genome. The remaining 34% with low classification probabilities are likely to be very similar homoeologs. Comparison to a subset of A, B and D genome SNPs confirmed 72% of A genome classifications and 85% of D genome classifications (Figure 2 and Supplementary Table 19). Discrimination of putative B genome genes was only ≈60%, possibly due to the use of cDNA sequences for classification when most of the informative sequence polymorphisms are intronic, and due to uncertainty about the ancestry of the B genome[5]. The set of 132,552 SNPs allocated to the A, B and D genomes is displayed using Brachypodium as a template in Figure 2, tracks 2-4.

There were no significant differences in the distribution of GOSlim molecular function categories in the A, B and D genes (Supplementary Figure 14), indicating that at this level of functional categorization there is no biased gene loss[36] in one of the genomes. Nevertheless, analysis of GOSlim terms associated with stop codons in A, B and D gene assemblies showed a strong tendency to retain functional copies of genes encoding transcription factors in all three genomes (Supplementary Figure 15), similar to the preferential retention of these genes in Arabidopsis genome duplications[37]. This indicates that genome-specific transcriptional regulatory networks tend to be maintained in wheat.

Conclusions

Using whole genome 454 sequencing, we assembled gene sequences representing an essentially complete gene set, and a significant number were assigned to the A, B or D genomes. Although the assemblies are fragmentary, they form a powerful framework for identifying genes, accelerating further genome sequencing, and facilitating genome-scale analyses. The identification of over 132,000 SNP in A, B and D genes facilitates QTL analysis and association studies of traits. Comparison with the sequences of diploid progenitors and relatives showed pronounced reductions in the size of large gene families in wheat despite the relatively recent formation of the hexaploid (Figure 3E), consistent with smaller scale analyses[28,38]. The scale of gene loss in hexaploid wheat compared to maize[36] and Brassica rapa[39] is significantly less, possibly due to its relatively recent origin and the absence of inter-genome recombination[40]. Nevertheless gene loss in wheat could be rapid, as shown in the newly created allopolyploid Tragopogon miscellus[41]. Most functional classes show equal gene loss in the three genomes, but families of transcription factors showed a clear tendency to be retained as functional genes in all three genomes. These may maintain transcriptional networks in each genome and contribute to non-additive gene expression[42] and genome plasticity. In contrast to the overall loss of gene family members, several classes of gene families with predicted roles in defence, nutritional content, energy metabolism and growth have increased sizes in the Triticeae lineage, possibly due to selection during domestication. Major efforts are underway to improve wheat productivity by increasing genetic diversity in breeding material and through genetic analysis of traits[43]. The genomic resources we have created promise to accelerate progress by facilitating the identification of useful variation in genes of wheat landraces and progenitor species, and by providing genomic landmarks to guide the progeny selection. Analysis of complex polygenic traits such as yield and nutrient use efficiency will also be accelerated, contributing to sustainable increases in wheat crop production.

Methods Summary

A single seed descent line of Triticum aestivum landrace “Chinese Spring” was sequenced as it is widely used for cytogenetic analysis[44] and physical mapping[15]. Triticum monococcum accession 4342-96 is a community standard line for TILLING, physical mapping and genetic analysis, and Aegilops tauschii ssp strangulata accession AL8/78, used for physical and genetic mapping, was sequenced using 454 technology48.

Sequence for T. aestivum wheat gene assembly was generated using Roche 454 Pyrosequencing on the GS FLX Titanium and GS FLX+ platforms. Additional sequence read datasets for T. aestivum, T. monococcum and A. tauschii were generated using 3 platforms: Illumina, 454 and SOLiD, to analyse homoeologous sequences and SNPs (a list of all datasets is in Supplementary Table 2). Orthologous Groups were created from rice, sorghum, and Brachypodium distachyon genome sequences and barley full-length cDNA sequences. Wheat gene assemblies were named according to their Orthologous Group Representative, and identified by a seven-digit identifier and their predicted genome viz: Traes_Bradi1g12345_0000001_D and Traes_Sb3g33333_6543210_A. Gene and cDNA assemblies can be searched at http://mips.helmholtz-muenchen.de/plant/wheat/index.jspp. Sequence assemblies are available for download from EBI accession PRJEB217 (OA: CALO01000001-CALO01945079; LCG: CALP010000001-CALP-15321847; cDNA ) and annotated gene sequences and their relationships can be viewed in a Brachypodium synteny-based Ensembl genome browser at: http://plants.ensembl.org/Brachypodium_distachyon.

Full Methods and any associated references are available in the online version of the paper at www.nature.com/nature.

Table 1b

Summary of assembly statistics of the Orthologous Assembly (OA), the LowCopy Genome (LCG) and cDNA assemblies.

	OA 99% mi1	LCG2	cDNA assembly2
Number of sequences	949,279	5,321,847	97,481
Total sequence (bp)	437,512,281	3,800,325,216	93,340,842
Minimum / maximum length (bp)	79 / 7,312	100 / 21,721	100 / 10,382
N10 / N50 / N90 (bp)	766 / 481 / 331	2,234 / 884 / 420	2,707 / 1,325 / 509
Mean length (bp)	460.89	714.10	957.53
GC content (%)	48.25	47.69	47.74

combined set of 454 sequences that cluster and form contigs and 454 sequences remaining singletons

only set of 454 sequences that cluster and form contigs