Shizhou Yu,Jianhui Wu,Mng Wang,Wiming Shi,Guangmin Xia,Jizng Jia,Zhnshng Kang, Djun Han

aState Key Laboratory of Crop Stress Biology for Arid Areas,Northwest A&F University,Yangling 712100,Shaanxi,China

bMolecular Genetics Key Laboratory of China Tobacco,Guizhou Academy of Tobacco Science,Guiyang 550081,Guizhou,China

cState Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008, Jiangsu,China

dKey Laboratory of Plant Cell Engineering and Germplasm Innovation, Ministry of Education, School of Life Sciences, Shandong University,Jinan 250100,Shandong,China

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

Keywords:Genome-wide association study Linkage disequilibrium Salt tolerance Haplotype tracing Triticum aestivum

A B S T R A C T Most modern wheat cultivars were selected on the basis of yield-related indices measured under optimal fertilizer and irrigation inputs. With climate change, land degradation and salinity caused by sea water encroachment, wheat is increasingly subjected to environmental stress. Moreover, expanding urbanization increasingly encroaches upon prime agricultural land in countries like China, and alternative cropping areas must be found.Some of these areas have moderate constraining factors, such as salinity. Therefore, it is important to investigate whether current genetic materials and breeding procedures are maintaining adequate variability to address future problems caused by abiotic stress.In this study, a panel of 307 wheat accessions, including local landraces, exotic cultivars used in Chinese breeding programs and Chinese cultivars released during different periods since 1940, were subjected to a genome-wide association study to dissect the genetic basis of salinity tolerance. Both marker-based and pedigree-based kinship analyses revealed that favorable haplotypes were introduced in some exotic cultivars as well as a limited number of Chinese landraces from the 1940s. However, improvements in salinity tolerance during modern breeding are not as obvious as that of yield. To broaden genetic diversity for increasing salt tolerance, there is a need to refocus attention on local landraces that have high degrees of salinity tolerance and carry rare favorable alleles that have not been exploited in breeding.

1. Introduction

It is predicted that the world population will reach 10 billion by 2050. To meet food demands, production of bread wheat(Triticum aestivum), a staple crop providing 20% of current calorie requirements, must increase by 60% [1]. This is a challenging task, as not only must wheat yields be increased at an unprecedented rate, it must be done by simultaneously addressing problems associated with climate change, land degradation and declining water availability and quality[2].

Genomic buffering in bread wheat allowed by its polyploid nature improves the chances of retaining natural variation.A dynamic and plastic genome allowed wheat to adapt to diverse environments and enabled global spread [3]. Hence,valuable genes for abiotic stress tolerance that are likely conserved in wheat gene banks need to be exploited. After decades of domestication and extensive selection in breeding,genetic variation in bread wheat has been eroded [4], making it more vulnerable to environmental stress. In modern breeding systems, and particularly in developing countries,yield has been the most important goal, and most modern wheat cultivars were selected based on yield-related indices measured under conditions of high nutrition and optimum field management [5]. Whereas such breeding systems were beneficial for maintenance of yield performance, they possibly caused unexpected and even negative effects on abiotic stress tolerance [6]. Therefore, it is worthwhile to compare tolerance characteristics and abiotic stress-related genetic diversity among landraces/early cultivars and modern cultivars, in order to make a better judgement on the current breeding system, and possibly gain clues for future improvement of abiotic stress tolerance in bread wheat[7].

Due to the complexity of the hexaploid wheat genome, it was previously impossible to perform certain traitassociated genetic diversity comparisons at the whole genome level [8]. Next-generation sequencing (NGS) technologies have radically changed that landscape.NGS enables high throughput discovery of DNA variants and facilitates faster development of new generation markers based on single nucleotide polymorphisms (SNP) [9]. With updated whole wheat genome sequencing and development of highdensity SNP chips, massive amounts of data on genetic variation are becoming more easily and more cheaply procurable. Several high-density and ultra-high-density wheat SNP arrays have been constructed for genetic studies and breeding [10].

These achievements greatly increased the number of genome-wide association studies (GWAS) in wheat. GWAS, a powerful tool for dissecting complex genetic variation and for estimating their effects on phenotype, has been widely used in crops such as rice and maize [11–13]. Various researchers have used high-density SNP arrays in GWAS to identify QTL for complex agronomic traits[14,15],such as yield[16,17]and biotic/abiotic stress tolerance [18–20]. Genetic diversity analyses of post-domestication evidence of selection and evolutionary histories have also been made [21–25]. All of these studies provide ways to trace the origins of beneficial alleles in breeding populations developed over time. However, few investigations have attempted to analyze changes in abiotic stress-related genetic diversity over the period of modern breeding; Thai is 1950–2019.

Bread wheat cultivation reached China about 4500 years ago, and then spread throughout the country except for subtropical southern areas [25]. Landraces, that evolved in different agro-ecological zones and presumably underwent diverse environmental adaptations, should be valuable genetic resources for modern breeding programs. As cross breeding in the west began more than half a century earlier than in China, many western bread wheat cultivars introduced from the 1930s formed another genetic resource during the early stages of Chinese cross breeding [26–28]. As the largest developing country with approximately 20% of the world population,yield was the top priority for Chinese wheat breeders. Fortunately, a large database of pedigrees and agronomic information regarding Chinese wheat cultivars is available and easily accessed [28]. Therefore, cultivars generated at different times are ideal for study of the effects of past selection for high yield on abiotic stress tolerance.

Soil salinity/alkalinity is one of the most serious stresses threatening crop productivity, especially in arid and semiarid regions worldwide [29]. Thus far, about 800 million hectares of arable land globally or more than 6% of the total have become saline or are endangered by salinity[30].China is one of the countries with the largest saline-alkaline areas,ranking third among the top 10 most affected countries.The saline-alkaline lands are distributed in the northwest,northeast, north, and coastal areas, and include 17 wheatgrowing provinces and more than 30 million hectares(http://www.stats.gov.cn/). With rapid urbanization on previously prime agricultural land, other areas less favorable for crop production are being considered for agriculture and the challenges they pose must be addressed [31]. Many studies attempted to address the problem of lack of salinity tolerance in wheat by genetic studies [32,33] QTL mapping [34–36] or association analysis [18,37–39]. However, few studies provided critical evaluations of the current Chinese breeding system based on comprehensive comparisons of genetic and phenotypic diversity within a broad representation of Chinese wheat germplasm.

Here, a collection of 307 wheat accessions, including local landraces, important exotic introductions and Chinese cultivars generated at different periods since 1940 was obtained from the China Agriculture Research System (CARS). Seed germination indices obtained under salinity stress were analyzed in a GWAS using the 660 K SNP array.The objectives of the study were to: 1) dissect and trace the genetic basis of salinity tolerance in a Chinese wheat panel,2)investigate the impact of the current breeding system with sole emphasis on yield on salinity tolerance,and 3)discuss strategies for future improvement of abiotic stress tolerance in wheat.

2. Materials and methods.

2.1. Plant materials

The bread wheat accessions used for tests in response to salt stress included 31 local landraces, 126 exotic cultivars(most of them past parents in Chinese breeding programs)and 150 cross-derived cultivars from the main wheat-growing areas in China (Fig. 1-a; Table S1). Among them, 16 landraces and 27 cultivars were from the Chinese mini core collection,that represents much of the genetic diversity (~70%) in 23,090 accessions held in the Chinese national collection [40,41]. The other materials were selected from more than 5000 accessions based on pre-diversity assessment using KASP markers or wheat 660 K SNP genotypic data for all 21 chromosomes. Some landraces were parents used in breeding modern cultivars.The modern cultivars were from different regions and released at different times. A second panel of 999 accessions with available genotypic data was used to validate the origin of favorable haplotypes.

Approximately 30 g of seed of each accession were planted in a 3 m2plot in a field nursery at Yangling (34°16′N(xiāo), 108°4′E)in early October 2016, and harvested in June 2017. After threshing, the seeds were fumigated and stored for 3 months in the granary of the College of Plant Protection. During the wheat growing season, we measured agronomic traits such as heading date, anthesis date, plant height, spike length,kernels per spike, thousand grain weight (TGW), and grain yield per plot. Cultivars Early Premium and Xiaoyan 22 were used as non-tolerant and tolerant controls throughout the study [42,43].

2.2. Phenotyping

In order to reduce variability associated with seed quality, we selected seeds with sound appearance and treated them before treatment as follows: a) soaked in 10% sodium hypochlorite (NaClO) solution for 10 min followed by three or four washes with water; b) soaked in 100 mmol L?1potassium nitrate (KNO3) solution for 4 h for osmotic adjustment,followed by at least four washes with water; and c) dried with an air-blower for 1 day at 25 °C [44].

After pre-processing, 30 seeds of each genotype were germinated in Petri dishes (90 mm diameter) with 20 mL of 0, 0.2, 0.4, 0.6, 0.8, 1.0 salt solution (standard concentration of seawater (SCS) in which almost all land plants cannot germinate) (Table S2) [45,46]. All seeds were incubated in a greenhouse at (25 ± 2) °C with a photoperiod of 14 h light(22,000 lx) and 10 h darkness at (20 ± 2) °C. Percentages of germinated seeds were recorded after 7 days. The tests were conducted three times during September and October 2017.

2.3. Salt tolerance index (STI)

The salt tolerance of a genotype at the germination stage was assumed to be a constant (θ) determined by genetic factors[47,48]. Wheat seeds were germinated in a dilution series of salt solutions (defined as variable, x) and germination rate was set as y. Since other environmental factors were treated as constant germination rate was the function:

The curve of y was a function of x. By integrating x from 0(pure water) to 1 (SCS). (The composition of seawater was formulated on the basis of Kester et al. [45] enabling Y to be estimated as:

Y is also a constant(Y ∝θ),thus if we obtain Y,we have evaluated θ.Here,we used the salt tolerance index(STI)to represented the actual level of salt tolerance(θ),namely STI ∈[0,1].

Applying an integral approach,STI was represented by the area under the differential concentration curve (AUDCC)calculated as follows:

2.4. Genotyping

The GWAS wheat panel was genotyped with the Affymetrix wheat 660 K SNP array(containing 630,514 SNP probes)[49]at Capital Bio Corporation, Beijing (http://www.capitalbiotech.com). SNP genotype calling and allele clustering was processed with the polyploid version of the Affymetrix Genotyping Console (GTC) software. To ensure genotyping accuracy, SNP markers with minor allele frequencies (MAF) <0.05, missing data >10%, and Hardy-Weinberg equilibria (HWE) >0.01, were excluded from further analysis. To obtain the physical positions of SNPs the software Nucleotide-Nucleotide BLAST(Basic Local Alignment Search Tool) 2.2.29+ was used to map the probe sequences to the IWGSC RefSeqv1.0 reference [50].The blast results were filtered with the parameters: ratio(AL/QL) ≥0.95; identity ≥0.95; E_value <1e?10. Polymorphism information content (PIC) values were calculated for each SNP marker using the formulawhere Piwas the proportion of the population carrying the ith allele[51].

2.5. Population structure

The p-distance was used to construct an NJ phylogenetic tree with 1000 bootstraps using the software MEGA-CC [52].Principal components analysis (PCA) of the population was performed using the software GCTA[53].Population structure was evaluated with 10,718 randomly selected SNPs from remaining SNPs using STRUCTURE V2.3.4 [54]. The model was applied without use of prior population information and the K-value representing the number of subgroups was set from 2 to 20; and five independent runs for each K were performed, with burn-in length set to 20,000 and iterations set to 10,000. The most likely number of sub-populations was determined using the delta K (ΔK) method based on the rate of change in LnP (D) between successive K-values, as previously described [55].

2.6. Estimation of linkage disequilibrium

Genome-wide linkage disequilibrium (LD) analysis for the A,B, and D genomes was performed on the 307 accessions using the software PLINK. LD was estimated by pairwise comparisons among filtered SNP markers using squared allele frequency correlations (r2). The parameters for calculating r2were set to –r2 –ld-window-kb 30,000 –ld-window 1000 –ldwindow-r2 0. For each subgenome, we calculated r2between pairs of SNP loci within a distance of 30,000 kb. For LD decay analysis, 30,000 kb LD regions covering SNP pairs were divided into bins of 100 kb, and r2values were plotted against distance in Mb. A LOESS curve was drawn to fit the data using seconddegree locally weighted scatter plot smoothing in the R program. The genetic distance corresponding to LD ≤0.1 was considered as the critical distance spanning a QTL.

2.7. Genome-wide association analysis

GWAS was conducted using a univariate linear mixed model in GEMMA software [56]. The suggestive threshold for P-values(P = 1/Ne) was calculated based on the modified Bonferroni correction; Ne represented the effective number of independent SNPs calculated using Genetic Type I Error Calculator software [57]. Our results showed suggestive Pvalue thresholds ranging from 1.29 × 10?4to 8.6 × 10?4for each chromosome (Table S3); we therefore considered 1 × 10?5as the criterion for genome-wide significance. Significant markers from GWAS were visualized using Manhattan plots,and important P-value distributions were visualized by quantile-quantile (QQ) plots, both drawn by the qqman package in R 3.0.3 (http://www.r-project.org/). The phenotypic variance explained by significant SNPs was evaluated by GCTA software [53].

2.8. Pedigree information

Pedigree information was sourced from publicly available breeders' records, genebank passport information (http://genbank.vurv.cz/ewdb/), information associated with variety release in China (http://www.cgris.net/cgris_english.html), the book titled “Chinese Wheat Improvement and Pedigree Analysis” [28], and with permission, from breeders' private records.Pedigrees were visualized using Helium version 1.18.03.15.

2.9. Candidate gene identification and expression

Based on IWGSC RefSeq v1.0 gene annotations (https://wheaturgi.versailles.inra.fr/Seq-Repository/Annotations), high confidence (HC) genes located within the LD block around significant SNPs were used for candidate gene analysis.Potential candidate genes or homologs in other plants involved in abiotic stress response were examined for gene expression. We downloaded transcripts per million (TPM)gene expression values for previously mapped RNA-seq samples on website http://202.194.139.32/ [58,59]. Transformed TPM values were visualized using the heatmap function in the NMF package in R v3.6.1.

Total RNA extraction and cDNA synthesis were performed following Wang et al.[73].Primers of TraesCS1A01G306300.1 for the qRT-PCR assay were F-5′-ACAGAGCCTTTGCTTAGTATGA-3′and R-5′-TCCACCGTCTTCTTGCTC-3′. TaEF1-α (M90077) was chosen as the reference gene for transcriptional profiling of wheat samples.Estimates of transcript abundance were based on four technical replicates made from each of three biological replicates.

3. Results

3.1. Genetic diversity, population structure and linkage disequilibrium

After removal of low quality SNP markers (MAF <0.05 and missing data >0.1), the remaining 431,673 SNPs were implemented in construction of a phylogenetic tree, and in principal components (PCA) and linkage disequilibrium anal-yses (Fig.S1). Population structure was also determined using STRUCTURE V2.3.4 with 10,718 randomly selected SNPs from the cleaned SNPs covering the whole genome and distributed evenly on all chromosomes. According to the ΔK method, a break in the slope was observed at K = 7; hence, the 307 accessions were divided into seven sub-populations (Sp1–Sp7)(Figs. 1-d, S2). To investigate phylogenetic relationships among the 307 accessions, genetic distances between clusters were computed from data generated with Wheat660K SNP markers, and neighbor-joining (NJ) trees were constructed (Fig.1-b). The NJ trees clustered the accessions into seven groups consistent with the sub-population number and result of principal components analysis (PCA) (Fig. 1-c). Sp1, Sp2, Sp3,Sp4, Sp5, Sp6, and Sp7 contained 11, 28, 34, 108, 54, 44, and 28 accessions, respectively (Table S4). Sp1 consisted mainly of old exotic wheat cultivars from East Africa and the Fertile Crescent; Sp2 mainly included exotic wheat cultivars cultivated during the past 70 years; Sp3 included mainly Chinese landraces; Sp4, Sp5 and Sp6 comprised mixed accessions from crosses between landraces and introduced varieties; Sp7 mainly included derivatives of crosses between modern Chinese cultivars.

The extent of LD was estimated using 173,066, 196,915, and 54,295 SNP loci from the A, B and D sub-genomes, respectively. LD analysis was assessed based on 717,701,068 pairwise comparisons from 424,276 SNPs, and pairwise LD was estimated using the squared-allele frequency correlation (r2).Plots of LD estimation (r2) as a function of physical distance(Mb) indicated clear LD decay with increasing physical distance (Fig. 1-e). In this study, a critical value for significance of r2was selected as 0.1, and the point in the curve that intercepted the critical r2was determined as the average LD decay. The average LD decay distance for the whole genome was approximately 3.2 Mb; the lowest LD decay distance of 2.0 Mb was in the D genome and about 2.7 and 4.8 Mb in the A and B genomes, respectively (Fig. 1-e).

3.2. Phenotypic data analysis

Xiaoyan 22 displayed higher salt tolerance with an STI of 0.9146 compared to Early Premium at 0.1256, suggesting that the salt treatments were appropriate and the results were reliable (Table S4). There was no significant variation among replicates with correlation coefficients of 0.98 (P < 0.01) (Fig.S3). The frequency distributions of STI values were approximately normal indicating that the responses were quantitative (Fig. S3). These results indicated that QTL for salt tolerance had large effects on germination percentage.

There was significant variation in salt tolerance (ANOVA P< 2e?16) among the 307 accessions (Fig. 2; Table S4). The detected phenotypic variation was confirmed by the values for range, standard deviation, and coefficient of variation(Table S4). Compared to the water control, all salinity concentrations significantly reduced germination. Based on the phenotypic heatmap, accessions were classified into two major groups (G1–G2) (Fig. 2-a): G1 showed a relatively stable germination potential under different salinity concentrations;b) G2 showed low germination rates at all salt concentrations or displayed a sharp drop with increasing of concentration. An STI value of 0.6 was considered the boundary between salinity tolerance and susceptibility (Fig. 2-b). We also compared STI values between subpopulations; there were significant differences in STI values between Sp1 and Sp3, Sp5 and Sp3, Sp6 and Sp3, however, there were almost no differences among the other subpopulations (Fig. 2-c). The proportions of accessions with superior salt tolerance in different subpopulations were consistent with the above results (Fig. 2-d). Based on the heat map for STI in individual accessions, there was a cluster of accessions highly sensitive to salt in Sp3; otherwise most accessions in Sp3 had low STI values (Fig. 2-c; Table S4).

3.3. Genome-wide association analysis

A total of 402,176 SNPs with average SNP density of 0.49 Mb were used in GWAS to detect QTL for salt tolerance. GWAS was performed using the univariate linear mixed model in GEMMA. One hundred and seventeen significant SNP traitassociations were detected with a suggestive threshold Pvalue <1.0e?4(Fig. 3; Table S5). Of these, 102 SNPs were clustered in three main genomic regions on chromosomes 1A(72), 3B (10) and 6B (20), respectively. Combining the 660 K SNP genetic and physical maps [60], most SNPs on 1A (map length 281.69 cM) were within the interval 152.08–154.05 cM corresponding to a physical interval of 2.85 Mb. SNPs on 3B (map length 475.88 cM) were in the region of 83.96 cM corresponding to a physical interval of 0.66 Mb. Most SNPs on 6B (map length 216.39 cM) were within the interval 42.15–43.85 cM corresponding to a physical interval of 0.39 Mb (Table S5). These QTL were designated as QSt.nwafu-1A, QSt.nwafu-3B and QSt.nwafu-6B and they contributed 8.93%, 5.23% and 1.98% of the phenotypic variation,respectively(Table 1).Many QTL for salt tolerance were identified in previous studies and were distributed on chromosomes 1AL, 2AS, 2BL, 2DS, 3AS, 3BS,4BL, 5BL, 6AL, 6BL, and 7AL. All of these QTL were associated with Na+exclusion,reduced Na+uptake,regulation of K+and/or Na+transport,or Ca2+and Mg2+accumulation[18,34–36,61].Compared with a previously reported QTL based on the physical positions of Chinese Spring reference genome, QSt.nwafu-6B was likely to be novel, whereas QSt.nwafu-1A and QSt.nwafu-3B were probably the same as QTL identified in Xiaoyan 54 [39]. This was consistent with pedigree analysis indicating that Xiaoyan 54 was derived from Xiaoyan 6.

Fig.2– Phenotypic clusters and genetic population structure and their relationships.(a)Similarity trees and heat-maps showing phenotypic (left)and genotypic relationships(right)of the traits.The phenotypic tree was constructed from the correlation matrices. (b)Two major clusters classified in the phenotypic tree (salt tolerant and salt sensitive) and boxplots of their STIs.(c)Accessions with STIs>0.6 for different sub-populations in the diversity panel.(d)Boxplot of STIs for different sub-populations in the diversity panel.

To obtain more genetic evidence on the relationship of the three QTL with salinity tolerance,we mapped these three QTL regions to the wheat reference genome (Chinese Spring,IWGSC RefSeq V1.0) and found a total of 53 high-confidence genes (Table S6). These genes were enriched in categories encoding transcription factors, calcium-dependent protein kinases, antioxidative enzymes, and osmoprotectants that have proved to be important in salinity tolerance [32].Moreover, when checking these genes in the bread wheat RNA-seq data [58], 35 genes were up- or down-regulated by salt stress (Fig. S4-a). One of the candidate genes for QSt.nwafu-1A, TraesCS1A01G306300.1, encoded a bZIP transcription factor annotated as an abscisic acid-insensitive 5-like protein (TaABI5-like) and its transcription was strongly affected by salt stress.Moreover, TaABI5-like expressed highly in seeds as well as young wheat seedings (Fig. S4-b). ABI5 plays a significant role in salt stress response in multiple plant species, particularly at the seed germination and young seedling development stages[62–64]. These results suggested that our GWAS results were reliable. Using a quantitative real time PCR (qRT-PCR) assay, we found that TraesCS1A01G306300.1 was up-regulated following treatment of wheat seeds and young seedlings with a 0.2 salt solution (Fig. S4-c). Moreover, expression in the tolerant cultivar Xiaoyan 22 was much higher than in the intolerant cultivar Bima 4.

3.4. Haplotype effects on salt tolerance and its confirmation by independent validation

Evaluation of local LD revealed that salt tolerance-associated SNPs were in haplotype blocks that we designated as Hap-ST-1A, Hap-ST-3B and Hap-ST-6B, with each having 72, 10,and 20 SNPs, respectively (Tables S5). The strongest associations exceeding or approximating a P-value < 1.0e?5(AX-111135532, AX-108794146, and AX-111029292) were located at 500.84 Mb in chromosome 1A, 70.60 Mb in 3B and 124.09 Mb in 6B (Fig. 3; Table S5). For Hap-ST-1A, there were four haplotype blocks, i.e. 1A-Block1 to 1A-Block4 based on combinations of SNP loci. Among them, 1A-Block1 represented the inferior haplotype (Hap-ST-1a); 1A-Block2 and 1ABlock3 represented favorable haplotype named Hap-ST-1A-1;and 1A-Block4 was a second favorable haplotype named Hap-ST-1A-2. There were two haplotype blocks for Hap-ST-3B; an inferior 3B-Block1 (Hap-ST-3b) and the favorable 3B-Block2(Hap-ST-3B). There were three haplotype blocks for Hap-ST-6B; 6B-Block1 was inferior (Hap-ST-6b) whereas 6B-Block2 and 6B-Block3 were favorable (Hap-ST-6B). To investigate the QTL effects, random accessions were divided into ten genotypic groups (designated as Hap_1a + 3b + 6b, Hap_1a + 3b+ 6B, Hap_1a + 3B + 6B, Hap_1A-1 + 3b + 6b, Hap_1A-1 + 3B + 6b,Hap_1A-1 + 3b + 6B, Hap_1A-1 + 3B + 6B, Hap_1A-2 + 3b + 6b,Hap_1A-2 + 3b + 6B, and Hap_1a + 3b + 6b +?) based on combinations of haplotypes for the three QTL (Table S7). The mean STI value (0.322) was significantly lower in accessions that carried all three adverse haplotypes (1a + 3b + 6b) than in accessions carrying three, two, or one favorable haplotypes(Fig. 4). Pairwise interactions between these loci showed significant additive effects (P = 0.05).

In a validation experiment to test the effect of the favorable haplotype combination Hap-ST-1A-1/Hap-ST-3B/Hap-ST-6B, 48 lines had the combination. A comparison of STI values of 20 random lines from this group with 20 Hap_1a + 3b + 6b accessions (Table S8) showed a highly significant difference (P = 0.0001) (Fig. S5-a, b, d). A pot experiment with the seeds sown in soil confirmed that Hap_1A-1 + 3B + 6B accessions had higher rates of germination and superior seedling growth than Hap_1a + 3b + 6b accessions when salt solution was added to the soil (Fig. S5-c).

3.5. Salt tolerant haplotypes traced by kinship and pedigree analysis

The proportions of different haplotypes were calculated for both panels of accessions. There was no difference in haplotype distribution frequencies between the two populations, indicating that the smaller selected population also represented the total genetic diversity (Fig. S6-a). Hap-ST-6B,the most frequent haplotype, was present in more than 75% of accessions whereas Hap-ST-3B, the least frequent haplotype,was present in less than 15% of accessions. Hap-ST-1A-1/2 were present in about 50% of accessions (Fig. S6-a). To trace the origins of favorable haplotypes in Chinese cultivars, we performed a combined kinship and pedigree analysis. A haplotype network analysis and comparison with population substructure revealed that carriers of favorable haplotypes were mainly from two sources. Hap-ST-1A-1 and Hap-ST-3B were mainly derived from introduced accessions whereas Hap-ST-1A-2 was common among Chinese landraces (Figs. 5 and 6). As indicated above, Hap-ST-6B was very common.Some Chinese landraces and their derivatives, such as Baihuamai, Youzimai, Beijing 10, Jimai 1 and Zaosui 30,harbored none of the favorable haplotypes, but displayed superior salt tolerance, indicating the likelihood of further rare genetic variants (Fig. 4; Table S2).

3.6. Salinity tolerance and genetic diversity in modern cultivars

To determine the effect of modern breeding on salinity tolerance, we divided the 307 accessions into different groups based on time of release. As expected, wheat yields increased with time; however, there was almost no difference in STI value among the different groups (Figs. 6, S7-b). Likewise,there was no difference in STI value when the accessions were divided into groups based on agro-ecological wheat zones (Fig.S7-c). These results indicated that little progress was made in improving the salinity tolerance of modern cultivars in Chinese modern breeding. A QTL in interval 42.15–43.85 cM on chromosome 6B with a large influence on peduncle length and kernel size in bread wheat cultivars from the Chinese Yellow and Huai Valley was reported [15]. Hap-ST-6B is closely linked with this QTL for yield-related traits, indicating that this genomic region had gone through strong directional selection in breeding presumably causing Hap-ST-6B accessions to be the most frequent haplotype. When accessions with Hap-ST-6B were compared with those lacking it, there were differences in multiple grain characters including seed width and thousand kernel weight (Figs. 6, S8-a, b; Table S9).Similarly, when accessions with STI >0.6 were divided into pre- and post-1980 groups the proportion of accessions with Hap-ST-1A-1 was higher in the latter group (Figs. 6, S6-b).Further comparative analyses of grain characters between accessions with Hap-ST-1A-1 and those with Hap-ST-1A-2,also showed differences in thousand kernel weight (Table S9).These results suggested that most of the Hap-ST-1A-1 and Hap-ST-6B carriers were selected for superior yield components, and that the favorable salt tolerance haplotypes were carried along by linkage(Figs.6,S8-c).

Fig. 3 – GWAS results for salt tolerance index (STI). (a, b) Q-Q and Manhattan plots of SNPs associated with STI. (c, d, e)Manhattan plots of SNP clusters on chromosomes 1A, 3B, and 6B and corresponding linkage disequilibrium (LD) r2 patterns. Dashed horizontal line depicts significance threshold level. SNPs beyond threshold are marked with light green dots. SNP names and corresponding local LD r2 value patterns are in the lower part of the graph. Numbers within the diamonds of the triangular LD matrix are r2 values multiplied by 100. (f, g, h) Haplotype genotypes in the candidate regions of Qst.nwafu-1A, Qst. nwafu-3B and Qst.nwafu-6B.

The origin of Hap-ST-1A-1 in Chinese cultivars was traced to the same backbone parents using combined kinship and pedigree analysis (Figs. 5 and 6). We calculated PIC values for different subpopulations and compared the genetic diversity between modern cultivar groups based on population structure using 52,303 markers. Sub-groups Sp4, Sp5 and Sp6 had relatively high mean PIC values (0.31, 0.30, 0.32, respectively)whereas Sp7 had a lower value (0.28). Compared with earlier cross-bred cultivars(Sp4,Sp5,and Sp6),the most recent group(Sp7) exhibited an even lower degree of diversity (0.28),agreeing with similar comparisons of genetic variation reported by others (Fig. 6) [26,27,65]. All of the above results provided evidence that modern cultivars are becoming homogenized as they frequently share common parents.

4. Discussion

4.1. Genetic basis of salinity tolerance in a Chinese modern wheat cultivar panel

Fig.4– Effects of combined salt tolerance(ST)haplotype variants Hap-ST-1A,Hap-ST-3B,and Hap-ST-6B on STI values from 307 wheat accessions.Box plots(quartiles are boxes,medians are continuous lines,means are crosses,whiskers extend to the farthest points that are not outliers,and outliers are black dots)for STI values with the identified haplotype and their combinations.

Wheat in China has a history that dates back more than 4000 years [26,28]. Initial wheat improvement was based on selection from landraces. The practice of cross-breeding and selection was introduced from the 1940s. This involved crossing between landraces and introduced cultivars [26]. In contrast, wheat improvement over the last 30 years was heavily based on crosses between elite varieties [26,28].Chinese landraces and introduced cultivars are the basis of improvement in tolerance to abiotic stresses. The present study used kinship and pedigree analyses to show that favorable haplotypes Hap-ST-1A-1 and Hap-ST-3B were derived from exotic wheat cultivars such as St1472/506, St2422/464,Ardito,Funo,and CIMMYT-derived lines(Fig.5;Table S1).These cultivars were widely adopted locally and/or were parents in cross-breeding [28]. Most of their derivatives, such as those prefixed as Yangmai, Yumai, and Zhoumai, carried favorable haplotypes or haplotype combination Hap_1A-1 + 3B + 6B, and showed high STI values. The carriers of the adverse Hap_1a + 3b + 6b combination were mostly derivatives of American accessions Quality and Early Premium,including Shijiazhuang 54, Beijing 8, and Jinan 9. By contrast,Hap-ST-1A-2 was inherited from Chinese landraces such as Yanda 1817, Jiyumai, Baihuomai, Sanyuehuang, Pingyuan 50,Jiangdongmen, and Youzimai (Fig. 5; Table S1). These landraces were widely grown before 1950 and some of them became leading parents in modern breeding.

Surprisingly, kinship analysis revealed that the genetic backgrounds of modern cultivars are more distinct from local landraces(Fig.2-a),indicating that landraces are treasure troves to broaden the genetic diversity of modern cultivars, a conclusion also reached by Hao et al. [26]. More importantly,some Chinese landraces and derivatives had the adverse haplotype 1a + 3b + 6b, but displayed superior salt tolerance,indicating that they were a source of rare favorable genetic variation (Fig. 4; Table S7). More attention should be given to those genotypes by breeders in order to capture the rare alleles.

4.2. Effects of the current yield-pursuing breeding system on salinity stress tolerance

Fig.5–Salt tolerance haplotypes traced by kinship and pedigree analysis in Chinese cultivars.Red orange and blue represent that the wheat accessions carry the favorable and adverse alleles for salt tolerance in this study,respectively;yellow means unknown.

Fig.6– Impact of selection for high yield on salinity tolerance in Chinese wheat breeding.The favorable salt tolerance haplotypes Hap-ST-1A-1 and Hap-ST-1A-2 were from exotic germplasms and local landraces. The percentage of accessions with Hap-ST-1A-1 was increased whereas Hap-ST-1A-2 was decreased after 1980 along with differences in thousand kernel weight(TKW).Hap-ST-6B is closely linked with a QTL for yield-related traits.The increasing frequencies of favorable haplotypes for salinity tolerance were due to linkage with yield traits in modern breeding.

Table 1–SNP markers significantly associated with salt tolerance index.

In order to meet future food production demands for the growing Chinese population, there is a need to produce cultivars adapted to more marginal land, such as the Bohai Gulf area in northeastern China where salinity is a major problem. It is important to investigate whether current breeding materials and procedures are maintaining adequate variability to adapt to this area.In this study,we identified the effects of the current breeding system on salt tolerance in leading Chinese wheat cultivars. Firstly, strong selective breeding for yield apparently has not enhanced the levels of salt tolerance in modern cultivars(Figs.6,S7).This was largely caused by selection and certification for high yield under optimal field conditions. Similar studies also showed that selective breeding under optimal fertilizer and irrigation inputs has resulted in diminished root mass[66,67].Secondly,there was no systematic selection for salt tolerance in breeding, and any favorable haplotypes for salt tolerance were retained by serendipitous linkage with yield traits (Figs.6, S8). Selection reduces variation at target loci but leads to hitchhiking effects around the selected loci;this raises linkage disequilibrium, reduces variability and causes a skewed distribution of allele frequencies at linked neutral sites [68].For example, selection for early heading in European wheat cultivars diminished genetic variation for root growth [67].The third impact of the current yield-pursuing breeding .system is the genetic bottleneck effect. Pedigree analysis in this study indicated that the frequency of favorable haplotypes for salt tolerance linked with yield traits trend towards homogenization due to common parents (Fig. 5; Table S1).There is increasing concern that breeding has drastically narrowed the overall genetic diversity in elite germplasm to an extent that there is insufficient variation in tolerance/resistance to specific stresses[26,69].In this study,derivatives of crosses between modern Chinese cultivars (Sp7) displayed lower molecular diversity (PIC, 0.28) than accessions derived from crosses between landraces and introduced varieties(Sp4, Sp5, and Sp6 with PIC values of 0.31, 0.30, 0.32,respectively). Obviously, this is not beneficial for future wheat breeding, especially for the enhancement of salt tolerance in modern wheat cultivars for the Bohai Gulf area.

4.3.Future targets and strategies in breeding salinity-tolerant wheat cultivars in China

Our data and the practical reality of the current Chinese breeding system place constraints on improvement of salinity tolerance in Chinese bread wheat cultivars.The first constraint is narrowing genetic diversity as most elite parents used in breeding, even by different breeding agencies in different regions, tend to share common relatives; the second is that most newly generated cultivars are assessed under optimal field conditions with sole emphasis on yield potential. Therefore,it is necessary to broaden genetic diversity by introducing useful alleles from external gene pools[4].Exotic germplasm is a valuable resource that contributed significantly to salinity tolerance in Chinese cultivars following its introduction from the 1930s[70].However,after almost a century of utilization,the favorable haplotypes introduced from exotic germplasm are already carried by most Chinese modern cultivars, and thus cannot make further contributions to improving salinity tolerance (Fig. S6-b). By contrast, favorable haplotypes carried by Chinese landraces, for example Hap-ST-1A-2, are present in only a limited number of modern cultivars (Table S7) and should be more widely incorporated in breeding programs. It should be noted that some landraces had a high degree of salinity tolerance but carried none of favorable haplotypes detected in this study (Fig. 4; Table S7). Moreover, since many of these landraces are already adapted to the growing environments in China [14,25], it could be worthwhile to use such germplasm for improvement of abiotic stress tolerance. Another important gene pool is the diploid and tetraploid ancestral species, which have been regularly used as a backbone in European elite varieties [71], and therefore make a large contribution to modern wheat genetic diversity. In the evolutionary history of wheat, particularly in the process of polyploidization and domestication, only a limited representation of the diversity of the progenitor species was transferred to the polyploid descendants [4]. Therefore, many favorable alleles could be locked in the ancestral gene pools. Although our study provides no information in this respect, the introduction of the salinity tolerance gene Nax2 from Triticum monococcum into a commercial wheat genotype greatly improved its tolerance [33].Additionally, species belonging to the tertiary gene pool comprise a number of grasses within the tribe Triticeae,harboring a treasure trove of potential alleles for stress tolerance. For example, tall wheatgrass (Thinopyrum ponticum),a halophytic relative of wheat, was used in intergeneric hybridization and contributed useful chromosome segments to modern wheat. For example, Shanrong No.3, an introgression cultivar with superior salinity tolerance, was allegedly generated by somatic hybridization between common wheat and tall wheatgrass [33,72,73]. Another successful outcome using tall wheatgrass was a number of Xiaoyan cultivars. In this study,Xiaoyan 6 and Xiaoyan 22 exhibit salinity tolerance, at least at the seedling germination stage.

It is understandable that yield along with resistance to a few key diseases is the most important trait for food stability.Under climate change and land degradation, it could be necessary to expand breeding objectives by greater emphasis on abiotic stress tolerance. One potential pipeline is to develop more practical and efficient phenotyping platforms to assess the tolerance levels of newly generated germplasm.Recent strides in high-throughput phenotyping technology will become increasingly important [74]. Moreover, as more favorable haplotypes and molecular markers associated with superior stress tolerance, such as those identified in this study, marker-assisted selection and ultimately genomic selection will enable this to occur [75]. Looking forward,more versatile avenues such gene editing and modification,might be applied when the candidate genes underlying stress tolerance are identified.

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

Author contributions

Shizhou Yu conducted the experiments, analyzed the data,and wrote the manuscript. Jianhui Wu and Meng Wang assisted in analyzing the data and wrote the manuscript.Weiming Shi, Guangmin Xia, and Jizeng Jia revised the manuscript. Dejun Han and Zhensheng Kang conceived and directed the project and revised the manuscript. All authors read and approved the final manuscript.

Declaration of competing interest

The authors declare no conflicts of interest.

Acknowledgments

The authors are grateful to Prof.R.A.McIntosh,Plant Breeding Institute, University of Sydney, for language editing and proofreading during the preparation of this manuscript. This study was financially supported by the National Youth Foundation of China (31901494, 31601306, and 31901869), the National Natural Science Foundation of China (31971890).Meng Wang was also supported by Young Elite Scientists Sponsorship Program of China Association for Science and Technology (2017QNRC001), and the Natural Science Fund of Jiangsu Province,China(BK20161092).