Zhiwei Wng, Chenyng Ho, Jing Zho, Chng Li, Chengzhi Jio, Wei Xi,Jin Hou, Tin Li, Hongxi Liu, Xueyong Zhng,*

a Key Laboratory of Crop Gene Resources and Germplasm Enhancement, Ministry of Agriculture and Rural Affairs/The National Key Facility for Crop Gene Resources and Genetic Improvement/Institute of Crop Sciences, Chinese Academy of Agricultural Sciences,Beijing 100081, China

b Novogene Bioinformatics Institute, Beijing 100083, China

Keywords:Wheat Polyploidization Differentiation Asymmetric recombination Haplotype block

ABSTRACT Common wheat(Triticum aestivum L.)is one of the most important crops because it provides about 20% of the total calories for humans. T. aestivum is an excellent modern species for studying concerted evolution of sub-genomes in polyploid species, because of its large chromosome size and three well-known genome donors.Establishment of common wheat genome reference sequence and development of high-density SNP chips provide an excellent foundation to answer questions of wheat evolution and breeding at the genomic level.By genotyping more than 600 accessions of common wheat and their diploid and tetraploid ancestors using a Wheat660K SNP array, we found dramatic genome changes due to tetraploidization and hexaploidization, in contrast to weaker influences of domestication and breeding on them. Further, since common wheat was introduced in China in 1500 BCE,Chinese landraces formed two subgroups(T.aestivum-L1 and T.aestivum-L2)with considerably diverse geographic distributions and agronomic traits.T.aestivum-L2,mainly distributed in central and east China is found to have more but smaller oval grains with early maturity characteristics.We found that variation and selection in intergenic regions of the A and B sub-genomes dominated this differentiation, in which chromosomes 7A and 3B took the leading roles due to the existence of putative genes related to defense responses and environmental adaption in the highly differentiated regions. Large haplotype blocks were detected on 3B (232.6-398.3 Mb) and 7A (211.7-272.9 Mb) in the landraces, forming two distinct haplotypes, respectively. We discovered that artificial crosses in breeding promoted recombination in the whole genome, however, this recombination and differentiation was highly asymmetric among the three sub-genomes in homoeologous regions. In addition, we found that the wide use of European and northern American cultivars in breeding at early era,led dramatic changes in Chinese wheat genome,whereas,the recent breeding functioned to optimize it.This study will provide the insight for reconsideration of wheat evolution and breeding, and a new strategy for parent selection in breeding.

1. Introduction

Common wheat (Triticum aestivum L.) is one of the most important cereal crops, which provides approximately 20%of the total calories for human intake [1]. Wheat was originated and domesticated in the Fertile Crescent region of central Asia. Tetra- and hexa-ploidization events caused reconstitute and epigenetic remodeling of wheat genome, making it an excellent modern species for studying concerted evolution of sub-genomes in polyploid species [2]. Recent technological advances in genomics have accelerated the generation of high-quality genome sequences of wheat and its progenitor species [3-7]. These genomic resources provide an excellent basis for development of high-density SNP chips to: (i) address how polyploidization and genomic changes occurred naturally; (ii) underpin the concerted evolution of sub-genomes; and (iii) detect driving elements for rearrangement and differentiation of sub-genomes after polyploidization. To date, the wheat 9K [8], 15K [9], 35K [10], 50K [11],90K [12], 280K [13], 660K [14], and 820K [15]SNP arrays,etc., have been developed and widely used in linkage map construction and QTL mapping [16-18], genetic diversity and phylogenetic analysis [19,20], genome-wide association[21,22], and selection signature and haplotype block survey[23,24]. More recently, a review article made a comprehensive comparison of seven widely used high-throughput wheat SNP arrays (e.g. SNP number, distribution, density, and associated genes), and suggested the Wheat660K SNP array with great potential for targeted genotyping and marker-assisted selection in common wheat [14].

Wheat was introduced into China about 1500 BCE. After long time spreading, cultivation and farmer selection, China was recognized as a wheat secondary diversity center with rich genetic variation because of vast difference in geographic,climate, cultivation ecosystem and cooking culture from west to east, and south to north China. In addition, Chinese wheat breeding program has achieved great success for feeding the rapidly escalating population. Since 1949, landraces in wheat production were gradually replaced by modern cultivars that were bred by crosses. More than 3000 new varieties of common wheat have been released over the past few decades[25], providing a platform to investigate genetic and evolutionary mechanisms of cultivars. In the current study, the Wheat660K SNP array was used to genotype 475 common wheat accessions, together with a set of tetraploid wheat,and T. urartu, Ae. speltoides, Ae. tauschii collections. With the advent of genome sequences of Chinese Spring (T. aestivum),we intend to gain a glimpse of genome evolution throughout farmer selection and modern breeding that highlights a new era of breeding to feed the increasing global population.

2. Materials and methods

2.1. Plant materials

A panel of 475 common wheat collections, comprising 158 Chinese landraces (CL), 297 modern Chinese cultivars (MCC),and 20 introduced modern cultivars(IMC)was used in present study. More than half of the materials were from Chinese wheat mini core collection, and the core contains 1% of Chinese wheat initial accessions but represents more than 70%of the total genetic diversity [26].

A set of 130 diploid wheats,including eight T.urartu,seven T. boeoticum, 15 T. monococcum, 41 Ae. speltoides, 18 Ae. longissima, seven Ae. sharonensis, and 34 Ae. tauschii, as well as 44 tetraploid wheats,including seven T.dicoccoides,15 T.dicoccum,16 T. durum and six T. turgidum, were also used to analyze the genetic effect of domestication, improvement, and polyploidization on genome-wide differentiation during wheat evolution. Detailed information of each studied plant material is listed in Tables S1, S2.

2.2. Phenotypic assessment and statistical analyses

Common wheat collections were planted at the Chinese Academy of Agricultural Sciences (CAAS) Xinxiang Experimental Station, Henan (113.5°E, 35.2°N) in four growing seasons (2014-2017) and once at the CAAS-Shunyi Experiment Station, Beijing (116.3°E, 40.0°N) in 2017, and they were named 2014XX, 2015XX, 2016XX, 2017XX, and 2017SY, respectively.Each accession was planted in a 2 m four-row plot with 25 cm between rows, and 40 seeds per row. Field management followed local practices. Ten plants from the middle of each plot were measured to phenotype eleven traits, including heading date (HD, day), flowering date (FD, day), spike length(SL, cm), spikelet number per spike (SN), plant height (PH, cm),kernel number per spike (KN), effective tiller number (ETN),1000-kernel weight (TKW, g), kernel length (KL, mm), kernel width (KW, mm), and kernel thickness (KT, mm).

Statistical analysis, including a two-tailed Student’s t-test was conducted in R software (https://www.r-project.org/) to analyze phenotypic trait differences between two subgroups of Chinese landraces classified based on SNP markers in the entire genome, as well as between modern Chinese cultivars released in different decades. The mean value of each trait was estimated by the best linear unbiased prediction (BLUP)method [27,28].

2.3. DNA extraction

Genomic DNA of each material was extracted from fresh leaves of 10-day-old seedlings(100 mg)in a single plant using DNA quick Plant System by Tiangen Biotech (Beijing) Co., Ltd(www.tiangen.com) according to the manufacturer’s instructions (Tiangen). The DNA with a final diluted concentration of 40 ng μL-1was used for genotyping.

2.4. Genotyping by a Wheat660K SNP array

The recently developed high density Wheat660K SNP array,also known as the Affymetrix Axiom Wheat660, contained 630,517 SNPs covering the entire genome[14].The Wheat660K was used to conduct genotyping by CapitalBio Corporation(http://www.capitalbio.com),and data capture was performed using the Affymetrix Analysis Suite system according to the manufacturer’s instructions (https://downloads.thermofisher.com/Axiom_Analysis_Suite_v_4.0.1_User_Guide.pdf)(Fig.S1a).Keeping in the varied ploidy level of different wheat germplasms, all materials were first grouped into diploid (A,B, and D sub-genomes), tetraploid, and hexaploid to capture data with DQC >0.6 and QC call rate>70%, with the penalties of four in both diploid and tetraploid, and 10 in hexaploid.Together with manually filtering (Fig. S1b), the final available number of SNPs in T.uratu(AA),Ae.speltoides(SS),Ae. tauschii(DD),tetraploid wheat(AABB),and common wheat(AABBDD)were 218,246, 210,872, 130,396, 467,351, and 582,381, respectively. All or part of genotype data will be available upon reasonable request.

2.5. Genetic structure evaluation

All accessions of natural populations were evaluated for population structure using 184,121, 202,934, and 62,305 SNPs(MAF > 0.05, Miss < 0.2) in wheat sub-genomes A, B, and D,respectively, which were used to infer three neighborjoining trees in the MEGA software with 1000 bootstrap replicates under the p-distance model[29].The phylogenetic trees were constructed using ITOL online software [30]with whole markers, genic markers, and intergenic markers. Division of markers to a corresponding chromosome was referenced using the Chinese Spring v1.0 Genome[7].Furthermore,principal component analysis (PCA) was performed to reveal the genetic relationship in the tested genotypes based on the aforementioned three subsets of markers from different genomic regions using the smart-pca program of EIGENSOFT v4.2 software [31], and the first three eigenvectors were plotted in a three-dimensional figure with R.

2.6. Selective sweep detection and gene flow calculation

Two parameters, π (nucleotide diversity) and π-ratio (ratio of nucleotide diversity),were used to depict the diversity pattern among subgroups. π measures the diversity of sub-genomes or one chromosome in each subgroup [32]. π-ratio is developed from π to estimate the reduction of diversity for a subgroup with respect to a control population [33]. In this study,the π-ratio was used to detect selective sweeps between subgroups in the R package popgenome [34]using a 30-SNP sliding window with a step size of 10-SNPs. Sliding windows with the highest 5%π-ratio values across different subgroups were picked as significant selection signals. Additionally, the R package CMplot (http://cran.uvigo.es/web/packages/CMplot/index.html) was used to depict not only selective sweeps, but also Manhattan plots.

FST, the genetic differentiation coefficient, was applied to reflect population differentiation on the basis of variances in allele frequency between subgroups [35], and was calculated by the R package popgenome [34]. Modern Chinese cultivars were divided into six different 10 year periods based on their released time i.e. 1950s, 1960s, 1970s, 1980s, 1990s and 2000s. Gene flow from Chinese landraces and introduced modern cultivars to each group in different breeding periods is calculated using the ABBA-BABA test (D statistic) [36]by detecting differences in allele sharing between two populations (P1 and P2) with a third population (P3). Here,Thinopyrum elongatum and Hordeum vulgare were used as multiple outgroup (O) [37,38]. The radar map of gene flow is drawn with the Excel software.

2.7. Identification of haplotypes in Chinese landraces

The identity scores (IS) [39]were calculated to visualize the shared haplotypes in Chinese Spring within 1-Mb windows.The IS was calculated using the following formula:

In which, Ds represents the difference between sample genotype and Chinese Spring in a single locus. If it was identical to CS,Ds was given value of 0;if it was different from the CS,Ds was given value of 1;and if it was heterozygous,Ds was given value of 0.5.

To determine the boundary of the haplotype blocks, we calculated r2of all SNPs pairwise in chromosomes 3B and 7A using PLINK [40], then the mean r2of each window(2 Mb)was calculated using a custom perl script.The window was defined as a haplotype block when its mean r2was greater than 0.5 in Chinese landraces.

2.8. KEGG pathway analysis

For candidate genes in the chromosomal regions with the highest differentiation on chromosomes 3B and 7A,they were enriched with rice pathways (http://www.kegg.jp) by KOBAS v2.0 [41], respectively. The filtering criteria of P < 0.05 and FDR < 0.05 were given for the analysis of KEGG enrichment analysis.

2.9. Recombination rate and asymmetry among subgenomes

Based on high-quality SNPs(MAF>0.05,Miss<0.2),the recombination rate of different subgroups was calculated using the software FastEPRR [42], with the sliding window of 10-SNPs and a step size of 2-SNPs, according to the formula ρ=4Ne×r.Here,Ne represents effective population,r means recombination, and ρ has a high positive correlation with r.Then, the collinear regions among A, B, and D sub-genomes were assigned and drawn based on the conserved protein sequence, using the software blast-n [43]and MCScanX [44].Collinear analysis data were obtained from conserved protein sequences of Chinese Spring reference genome [7]. Finally,the recombination rate among sub-genomes was intercompared in the regions of high collinearity to reveal asymmetry among the three homoeologous chromosomes.

3. Results

3.1. Integration of the D sub-genome broke the stable genetic partnership between the A and B sub-genomes in emmer wheat

A total of 475,995 effective SNPs were obtained in 158 Chinese landraces and 317 modern Chinese cultivars after the original data was filtered, in which 191,240, 208,592, and 76,163 polymorphic SNPs were physically mapped in the A, B, and D sub-genomes, respectively, based on Chinese Spring reference genome [7](Table S3). To reveal genome differentiation among diploid, tetraploid, and hexaploid wheat, we also genotyped the three donor species, including 8 T. urartu, 41 Ae. speltoides, 34 Ae. tauschii, 7 T. dicoccoides, and 36 T. dicoccum accessions. We filtered the genotypic data firstly and then used SNP markers shared among common wheat, tetraploid wheat, and the three diploid donor species (Table S4; Fig. S2)to reveal genome changes.

Significant variations occurred between the T. dicoccoides and their donor genomes of T. urartu and Ae. speltoides. The average FSTvalues reached 0.553 and 0.519 for the A and B sub-genomes, respectively. In total, 6.95% and 7.92% of 30-SNP windows reached 0.8, respectively. However, mean FSTvalues were only 0.160 for A sub-genome and 0.125 for B sub-genome between T. dicoccum and T. dicoccoides. This indicates that dramatic changes occurred between the two subgenomes and their potential donors.But domestication of tetraploids produced a much weaker effect on the sub-genomes in comparison with tetraploidization (Fig. 1; Tables S5, S6).

About 7000-8000 years ago, the domesticated T. dicoccum naturally crossed with Ae. tauschii and created the hexaploid wheat T.aestivum L.,also known as hexaploidization.FSTanalysis indicated that hexaploidization not only induced differentiation of the D sub-genome from its donor species, but also affected the stable genetic partnership between the A and B sub-genomes formed in emmer in the long-term evolution. Interestingly, a new round of genetic differentiation occurred in the A and B sub-genomes of common wheat(landrace)in comparison with that of T.dicoccum.For example,at 30.53% and 30.12% of the 30-SNP windows, the FSTvalues exceeded the cutoff of 0.5; at 2.09% and 0.69% of the 30-SNP windows, FSTvalues even reached 0.8 in the A and B subgenomes,respectively.As expected,D sub-genome differentiation from its donor species was much stronger than in the A and B sub-genomes from those in T. dicoccum: 36.17% of window FSTvalues between D sub-genome and its donor species exceeded the cutoff value of 0.8 (Fig. 1; Tables S5, S6).

Further, the π-ratio between species clearly showed that the distribution of selective sweeps was uneven among subgenomes,even within single chromosome during wheat polyploidy. For example, in tetraploidization, strong selective sweeps occurred on 3A (223-372 Mb), 3B (203-252 Mb), and 5B (496-519 Mb),while strong selective sweeps were detected on 4A (626-701 Mb), 2B (250-350 Mb), and 4D (306-354 Mb) in common wheat (Fig. S3). In domestication from wild emmer wheat to cultivars (T. dicoccum), strong purification selection mainly occurred on 5A,6A,7A,3B,5B,and 6B.In the transformation of common wheat from landraces to modern cultivars, the A and B sub-genomes experienced stronger selection than the D sub-genome, especially in 1A, 4A, 5A,2B, 5B, 6B, and 7B (Fig. S3).

3.2. Very strong geographic and genetic differentiation occurred in Chinese wheat landraces

Cluster analysis based on SNP markers in different subgenomes, even in the genic, and intergenic regions showed that Chinese landraces were clearly clustered in two subgroups, T. aestivum-L1 and T. aestivum-L2, and SNPs from the intergenic regions on B and A sub-genomes mainly drove the distinct differentiation between the two lineages(Figs. 2, 3a; Fig. S4). Further, chromosomes 3B and 7A are mainly responsible for genetic differentiation in the landraces,which was clearly indicated by the FSTvalues between T. aestivum-L1 and T. aestivum-L2 (Fig. 3b). At the critical regions on 3B and 7A,FSTreached 0.839 and 0.614,respectively(Fig. 3c, d). These regions potentially harbor candidate for grain size and length based on multi-environmental BLUP data (Fig. 3e, f). The re-cluster results based just on the SNPs within the critical genomic regions on 3B and 7A separately remained unchanged compared with SNPs in the entire genome, further supporting the leading role of 3B and 7A in genetic differentiation of T. aestivum-L2 from T. aestivum-L1(Fig. S5).

Further,haplotype block map was constructed using polymorphic SNPs on chromosomes 3B and 7A(Fig.4a,c).Notably,an obvious block was existed in the two chromosomal regions leading to highly differentiated T.aestivum-L1 and T.aestivum-L2, and the existence of two different haplotypes in single block (232.6-398.3 Mb on 3B and 211.7-272.9 Mb on 7A) in the two subgroups of Chinese landraces, further demonstrated that the two genomic regions played the crucial role in differentiation of the two landrace subgroups. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of these two high-differentiation genomic regions showed the putative involvement of genes related to phenylpropanoid biosynthesis,porphyrin and chlorophyll metabolism,arginine and proline metabolism and plant-pathogen interaction,which are related to defense responses and environmental adaption (Fig. 4b, d; Table S7).

Accessions in T. aestivum-L1 were mainly distributed in northwestern China, whereas those in T. aestivum-L2 were mainly from central to eastern China (Fig. 3a). The most distinct agronomic trait between the two subgroups was grain size, i.e. the T. aestivum-L2 accessions usually had smaller grain size than the T. aestivum-L1, which was achieved by reduction in grain length (Fig. S6). In addition, T. aestivum-L2 accessions have shorter duration cropping seasons, shorter spike length, and fewer spikelets than T. aestivum-L1 (Fig. S6;Table S8). The πL1/πL2and πL2/πL1ratios also illustrated that genes controlling grain size, such as CWI-5D [45], MOC1 [46],and TFE-7A [47], etc., were genetically involved in seed size differentiation (Fig. 3g, h). The π ratio also indicated that T.aestivum-L2 has experienced stronger selective pressure than T. aestivum-L1.

3.3. Early breeding in China led dramatic genome change while recent breeding subsequently optimized it

Fig. 1 - Both tetraploidization and hexaploidization induced much stronger diversification than either domestication of tetraploids or modern breeding of common wheat. Diversification index (FST) between species, the modern cultivars, and landraces of common wheat, which were estimated based on the effective SNPs in the A-, B-, and D- sub-genomes,respectively. The color dots means FST for each sliding window of 30-SNP markers. For each of the 21 chromosomes, FST values are given in the grey shaded area at the top of each chromosome.The blue and red dashed lines are the cutoffs for 0.5 and 0.8, respectively.

First,genetic diversity(π value)in the different breeding periods indicated that modern Chinese cultivars (MCC) in 1970s had the highest diversity, followed by the ones in 1980s(Fig. 5a). Further, based on the breeding periods, gene flow between MCC,Chinese landraces(CL)and introduced modern cultivars (IMC) were calculated, respectively, it showed a higher gene flow from CL to MCC released before 1970s, but from IMC to MCC released since 1980s (Fig. 5b). To highlight alterations at the chromosomal level more clearly, we simply compared genetic differentiation of the cultivars released in the 1950s, 1980s, and 2000s. The agronomic traits in multiple environments clearly indicated that the heading date of MCC has become earlier, while plant height has been significantly reduced. The plant height reduction from the 1950s to the 1980s was extremely significant,which matches with the first global green revolution. However, the length, thickness, and width of grains (1950s-2000s) have continuously increased,which led to the increased 1000-kernel weight and grain number per spike(Fig.5c,d;Table S9;Fig.S7).As a result,a significant increase in grain yield occurred in the past 70 years[25].

When we tested the genome differentiation of cultivars released in the three aforementioned periods, we found significant differentiation of cultivars released in the 1980s with those released in the 1950s on all 21 chromosomes.However,differentiation of cultivars released in the 2000s compared with those in the 1980s was relatively weaker(Fig. 5e, f). This indicates that wide use of European and American cultivars as parent to cross with local cultivars in breeding from 1940s to 1980s, led dramatic change of genome and subsequent breeding program has further optimized the genome of the cultivars (Fig. 5g, h).

3.4. Breeding has promoted effective recombination with an asymmetric distribution among three homoeologous chromosomes in wheat

Artificial crossing has promoted effective recombination,because along all of the 21 chromosomes, the recombination frequency in the modern cultivars is higher than that in the landraces (Fig. 6a, c; Fig. S8). The most surprising finding was an asymmetric sub-genome distribution of recombination frequency in some homologous regions. In other words,usually one sub-genome has a higher recombination frequency than the other two at the homologous regions(Fig. 6a, c), such as 218.9-267.2 Mb on 1A and 249.7-377.5 Mb on 2B,etc.(Table S10).The πL/πMratio also showed an obvious asymmetry among the three homoeologous chromosomes(Fig. 6b, d; Fig. S9). This might mean that in differentiation between the modern cultivars from the landraces, one genome contributed more than the other two at the homologous regions among the three sub-genomes, such as 396.4-447.7 Mb on 1A and 419.4-453.6 Mb on 2B, etc. (Table S11).Notably,it could be concluded that A and B sub-genomes carried more selective sweeps than the D sub-genome. These results showed that an asymmetric distribution of recombination frequency and selective sweeps was present among the three sub-genomes and homoeologous chromosomes in common wheat.

Fig. 2 - Cluster analysis of wheat accessions based on all markers, genic markers, and intergenic markers on three subgenomes. Light blue lines are diploid accessions (A sub-genome), orange lines are diploid accessions(B sub-genome), black lines are diploid accessions (D sub-genome), green lines are emmer wheat, dark blue lines are landraces, dark red lines are modern cultivars and yellow lines are introduced cultivars.The intergenic SNPs of B sub-genome divided Chinese landraces into two subgroups.

4. Discussion

4.1. Mechanisms for integration of different genomes to form new polyploid species

Polyploid organisms contain multiple sets of chromosomes.Genomic studies of polyploid organisms illustrate that the evolutionary flexibility contributed by polyploidy has reshaped the genomes of most eukaryotes [48], which was also demonstrated by strong genetic differentiation between ploidy in the current study (Fig. 1). Genome reshuffling is likely to be a crucially important source of genetic diversity in polyploid species, especially in the earlier stages of polyploidization [49-51]. Loss of duplicated genes and new functionalization are common mechanisms in polyploid evolution[52,53].A high proportion(84%)of transposable elements(TEs,transposons)may be used to contribute to wheat improvement.TEs provide a source of allelic regulatory variation in genes that are implicated in stress adaptability, and might facilitate rapid genome restructuring after polyploidization, because TEs can jump from one of the subgenomes to the others [54-56]. The cluster images based on intergenic SNPs are quite similar to those based on whole genome SNPs (Fig. 2). This result suggested the influence of repetitive DNA sequences in ecological differentiation [7,57].

Fig. 3 - Very strong geographic and genetic differentiation occurred and 3B and 7A took the leading roles for genetic differentiation in Chinese landraces. (a) T. aestivum-L1 subgroup is mainly distributed in north-western and north-eastern China,but T.aestivum-L2 is distributed in central-eastern China.(b)FST values indicate that 3B and 7A lead the differentiation of T. aestivum-L2 from T. aestivum-L1. (c & d) The FST values reached 0.514 and 0.318 on whole chromosomes 3B and 7A,respectively;at the critical regions they reached 0.839 and 0.614 on 3B and 7A,respectively,which are much higher than the whole genome differentiation between T.aestivum-L1 and T.aestivum-L2(0.21).On 3B and 7A,as well as their critical regions,the differentiation indexes between the modern Chinese cultivar (MCC) and T. aestivum-L2 are higher than those for the whole genome. (e & f) The T. aestivum-L2 usually has shorter and smaller grains than T. aestivum-L1. The grain length and 1000-kernel weight difference was significant between these two subgroups in all environments.(g&h)The πL1/πL2 and πL2/πL1 clearly indicate that genes related to grain size are selected in this geographic and genetic differentiation. *, P < 0.05, **,P < 0.01, ***, P < 0.001, ns means not significant.

4.2. Contrasting selection criteria for grain size in ancient farmer selection and modern breeding

Clustering based on the SNPs in an entire genome and subgenomes clearly indicated that A and B sub-genomes, and especially their intergenic regions, genetically drive differentiation of T. aestivum-L2 from T. aestivum-L1 (Fig. 2). The two subgroups of landraces have divergent haplotypes in the critical genomic region of 3B and 7A (Fig. 4a, c). Appearance of a few T. aestivum-L2 collections in north China also supported their evolution from T. aestivum-L1. The most interesting thing is that accessions of T. aestivum-L2 usually flowers and mature earlier and thus have smaller grains compared to those in T.aestivum-L1(Fig.3;Fig.S6).The selection for earlier heading date per se should be related to the cultivation system, i.e. central-eastern China usually has two crop seasons each year, which was consistent with the priority of early maturing character in modern breeding (Fig. 5c). The priority of producing smaller grains might be related to the cooking style in ancient times.Ancient Chinese used to eat all cereals as full grains rather than milling wheat into flour [58,59]. In contrast to ancient trends, modern-day breeding has prioritized the selection of larger grains since 1950s (Fig. 5d), due to technological development in the milling industry as well as ample electricity. At present, more than 90% of wheat grains are used to make steamed bread, pan-bread, noodles,and dumplings. Strong progressive selection for larger grain size is also related to the rapid development of irrigation systems and the fertilizer industry after 1950s. In addition, the subgroup T. aestivum-L2 is an insufficiently exploited gene pool for breeding, because it is genetically distinct from the modern cultivars (Fig. 3c, d).

Fig. 4 - Haplotype blocks and KEGG enrichment analysis on chromosomes 3B and 7A. (a & c) A large haplotype block was detected in the high-differentiation chromosomal regions in Fig. 3b in Chinese landraces. Obviously, the two subgroups of landraces,T.aestivum-L1 and T.aestivum-L2,carried almost different haplotypes in the two strongly differentiated regions on chromosomes 3B and 7A. Each row is a genotype of Chinese landraces, and each column is a haplotype. Alleles that are identical to or different from the ones of Chinese Spring are indicated by blue and red bars, respectively. And the missing allele is labelled with white color. (b & d) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the key genomic regions, that is 232.6-398.3 Mb on chromosome 3B and 211.7-272.9 Mb on 7A.

4.3. Asymmetric sub-genome evolution and selection occurred in polyploid crops

Asymmetric selection to subgenomes was frequently detected in polyploid crops,such as in upland cotton(Gossypium hirsutum), canola (Brassica napus) and wheat. This has become an important hotspot with advances in crop genomics. For examples. In the new synthetic canola (Brassica napus),rapid gene loss and expression divergence were found in the Anand Cnsub-genomes, and selection for lower erucic acid and high-oil accelerated the loss of glucosinolate genes,while preserving expansion of oil biosynthesis genes [49-53]. Further, resequencing 588 canola accessions uncovered that during the first stage of improvement(FSI),A subgenome specific selection promoted the stress tolerance and oil accumulation,while asymmetrical selection in the C sub-genomes improving developmental traits [60]. In upland cotton, A-subgenome was selected specifically for fiber quality, but D sub-genome for biotic and abiotic stress tolerance[61].Moreover, resequencing of 352 wild and domesticated cotton accessions revealed asymmetric selection for long fibers in the A sub-genome, but for fiber color in the D sub-genome[62].

In wheat, through wild wheat domestication QTL mapping, Peng [63]found that there were more QTLs controlling domestication syndromes in the A sub-genome than in the B sub-genome in evolution of tetraploid wheat from wild emmer to cultivated emmer. Furthermore, in tetraploid wheat, the B sub-genome showed higher plasticity [64]. In common wheat, the D sub-genome showed as the most conserved and consistent with its donor genome [65]even if the dramatic genome change caused by hexaploidization was inevitable (Fig. 1). Another possible reason for the A and B sub-genomes having more selective sweeps than the D subgenome was frequent genome fragment introgressions from wild emmer to common wheat, which was clearly revealed by recent genome resequencing of wheat and its relatives[66,67]. Importantly, the second step hybridization between tetraploid wheat and Ae. tauschii happened about 7000-10,000 years ago [68], and the history of selection on the D sub-genome was much shorter than that on A and B subgenomes. Considering the advances in various omics techniques with reducing cost, it would be of interest to find out most of the functional mechanisms of asymmetric subgenome evolution and selection in bread wheat.

Fig. 5 - Agronomic trait difference and genome change caused by crosses and selection in modern breeding. (a) Genetic diversity π value (10-6, in green line) of modern Chinese cultivars (MCC) in different breeding periods. (b) Gene flow from Chinese landraces(CL,in blue line)and introduced modern cultivars(IMC,in red line)to MCC released from 1950s to 2000s.(c&d)MCC released in the 1950s,1980s,and 2000s showed large differences in heading date and grain size.(e&f)At genome level, high differentiation was detected between cultivars from the 1980s and the 1950s; differentiation between cultivars from the 2000s and 1980s is not so obvious. (g & h) π ratio (π1950s/1980s and π1980s/2000s) indicates much stronger selection occurred in breeding from 1950s-1980s than 1980s-2000s. The two periods targeted different chromosomes or genomic regions of the same chromosome.

4.4. Wheat breeding in the post-genomics era

Selective breeding leads to reduced diversity in most crops[69].However, artificial crossing in wheat breeding has promoted recombination and gene exchange among gene pools, which results in breaking the bottle-necking phenomenon to some extent[24].This might be related to the self-pollination character of wheat. Post 1960s in China, modern breeding selection resulted in significant differentiation at the genome level(Fig. S9) and development by the state of the irrigation and the fertilizer industry, which promoted a rapid Green Revolution of modern cultivars during the 1970s-1980s. Another important reason causing increased genetic diversity in modern cultivars is the wide use of genetic resources from abroad,especially those from the former USSR, Italy, and the USA,which has been illustrated with multi-functional genes[70].

Fig. 6 - Asymmetric distribution of recombination and genetic diversity differentiation among the three sub-genomes in landraces and modern cultivars. The modern cultivars have a higher recombination frequency than the landraces. On the grey regions anchored by homologous genes,one or two genomes usually have higher RF than the others(a&c).One of the three homology regions usually leads the differentiation of modern cultivars from the landraces (b & d).

Since 1980s, Chinese wheat yield potential has continued to increase, even if differentiation at the genome level is not very significant between cultivars released in the 1980s and 2000s (Fig. 5f). This means that modern-day breeding has optimized the wheat genome.Therefore, it is crucial to search out cultivars with the most optimized genome based on the variations of essential haplotype blocks [71], further improve the block types [72], and then use them as recurrent parents in crosses with exotic collections conserved in gene banks,which will be a good strategy to bridge the gene banks with modern breeding in the post genomic era for wheat improvement [73,74].

Declaration of competing interest

Authors declare that there are no conflicts of interest.

Acknowledgements

We are grateful to Professors Rudi Apples from University of Melbourne and Long Mao from ICS-CAAS for their critical discussion of this manuscript. This study was supported by the National Key Research and Development Program of China(2016YFD0100302),the CAAS Program(Y2017PT39),and Jiangsu Collaborative Innovation Center for Modern Crop Production.

Appendix A. Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2020.08.006.