Responses to salinity stress, variations and correlations between phenotypic traits
The phenotypes of 191 wheat accessions were characterized during three crop seasons (2014–2017) in LS and HS environments. Significant differences existed among accessions for all the traits in the two different salinity treatments. Compared with the LS conditions, the mean values for the eight phenotypic traits decreased to different degrees in the HS treatment, indicating more significant inhibitory effects on the growth of wheat accessions under HS conditions (Table 1). Additionally, correlation co-efficient analyses among the three field experiments for the same trait under the same treatment were conducted. All the adult-stage traits had highly significant positive correlations, implying the moderate repeatability for these traits in field experiments (Table S3). Estimations of heritability for all the phenotypic traits showed moderate to high heritability levels. All the phenotypic traits showed higher heritability levels under LS conditions than under HS conditions (Table 1). PH showed the highest heritability levels under LS (H2 = 88.58%) and HS (H2 = 56.30%) conditions. SL and KPS have moderate and relatively stable heritability levels, with H2 values of 56.83% and 60.03% under LS conditions and 46.52% and 47.86% under HS conditions, respectively. Compared with the other traits, YPP and GN showed similar but lower heritability level under LS and HS conditions (Table 1). Pearson’s co-efficient of correlation between traits was calculated based on data averaged across three years under the two salinity conditions (Table 2). PH showed significantly positive correlations with SN, SL, KPS, BM, YPP and GN under both salt-stress conditions, and the correlation coefficients between PH and both BM and GN were comparatively higher than those between PH and other traits in the two salinity environments (Table 2). TKW showed significantly negative correlations with PH, SN, KPS and GN, but was positively correlated with SL and YPP under both salinity treatments (Table 2). In addition, BM and GN showed extremely high and positive correlations with YPP under LS and HS conditions (Table 2). The correlations between yield and its components (SN, KPS and TKW) were significantly positive under both salinity conditions (Table 2). Importantly, the correlation coefficients between YPP and SN were highest under both conditions, while the correlation coefficients between YPP and TKW were lower, compared with those between YPP and KPS (Table 2). Interestingly, compared with the LS treatment, the correlation coefficients between YPP and both SN and KPS slight increased with the HS treatment (Table 2).
Based on the quality preprocessing results provided by PLINK software (http://zzz.bwh.harvard.edu/plink/data.shtml ), 8,451 and 221,878 SNPs were removed for having call rates less than 80% and minor allele frequencies less than 0.05, respectively. No accession was removed for low genotyping (MIND > 0.1). Finally, 191 accessions containing 387,657 SNPs were retained for further analyses. And among these, 322,590 SNP markers were mapped on the Wheat660K consensus maps (Prof.Jizeng Jia, pers.comm.). The numbers of SNPs presented on each wheat chromosome were provided in Table S4. In total, 740,124 variation-containing alleles were detected by 387,657 SNPs among the 191 cultivars (Table S4). The range of the mean allele numbers on the 21 wheat chromosomes was from 1.38 to 3.25. The greatest mean allele numbers were on chromosomes 6D, 2D and 4A, successively. Chromosome 3B had the lowest mean allele number (Table S4).
Genetic structure and linkage disequilibrium
Different methods were used to analyze the population structure. The neighbor-joining method based on Nei’s standard genetic distance  was used to classify the 191 accessions, and it indicated that they were divided into two groups (Fig. 1a). The first group (45) mainly consisted of improved varieties from Shaanxi. The second group (146) was mainly comprised of improved varieties originating from Beijing, Hebei, Henan and Shandong. These two groups were further confirmed by a PCA and population structure analysis. When the number of subpopulations (Ks) was plotted against the ΔK calculated using software Structure, the highest ΔK was observed at K = 2 (Fig. 1b and 1e), indicating that the 191 accessions could be divided into two subgroups. The PCA also showed two distinct clusters (Fig. 1c). According to the phylogenetic tree constructed using the neighbor-joining method, the well-known salt-tolerant wheat accessions, such as Cang 6001, Cang 6002, Cang 6005, Dekang 961, Jimai 32 and Keyi 26, clustered into group 2.
Linkage disequilibrium was calculated using 153,725, 184,755, and 49,177 SNP markers for the A, B, and D sub-genomes, respectively. Based on the critical r2 threshold value at 0.1, the LD decay distance of the whole genome was approximately 5.0 Mb. The highest LD decay distance of 8.5 Mb was found in the A sub-genome, followed by 8.0 Mb in the B sub-genome. The D sub-genome has the lowest LD decay distance of 3.5 Mb (Fig. 1d).
Associations between traits and SNPs
Association analyses were performed between the phenotypic traits at the adult stage and SNP markers. In total, 611 significantly associated signals among 389 SNPs were identified for eight phenotypic traits in the six environments at P < 1.25e-6 (Table S5). The phenotypic variation explanation rate (R2) ranged from 9.14% to 68.81% (Table S5). A total of 152 SNPs showed significantly associations with the same trait in two or more environments, and 44 SNPs showed significantly associations with two or more traits (Table S6). In general, 11 significantly associated loci (a 5 Mb interval containing three or more significant associated SNPs was regarded as a locus) that might affect the growth of wheat under salt-stress conditions were among the 389 SNPs associated with PH, SL, SN, YPP, GN, BM and TKW (Fig. 2).
Up to 141 significant SNPs associated with PH were detected and mainly distributed on chromosomes 2B, 2D, 3D, 4A, 4D, 5A, 6B, 7B and 7D (Table S5). Of these, 79 SNPs were detected in at least two environments under LS conditions and were mainly located on chromosome 5A (Table S6). Significant associations were detected in one environment under HS conditions for eight SNPs related to PH, which were located on chromosome 4A (2), 4D (4) and 5A (2) (Table S6). One SNP locus named QPh-5A1, which contained four SNPs and was located on chromosome 5A at 535.91–536.80 Mb, was significantly associated with PH in 15LS and 15HS, explaining 44.46%–67.67% of the phenotypic variation. More importantly, one SNP of QPh-5A1 was significantly associated with YPP in 15HS. On chromosome 4D, four SNPs located in the interval 16.93–19.29 Mb (Fig. 2), defined as QPh-4D, were associated with PH in 17HS and explained 52.81%-53.62% of the phenotypic variation.
A total of 80 significant SNPs were associated with SL and distributed on chromosomes 1A, 1B, 2A, 2D, 3B, 3D, 5B, 5D, 6D and 7A (Table S5). Among them, nine SNPs on chromosome 1B and three SNPs on 7A were detected in two environments under LS conditions (Table S6). Meanwhile, four SNPs on chromosomes 1A (1), 1B (1), 2A (1) and 2D (1) were detected in two environments under HS conditions (Table S6). The important loci associated with SL were mainly located on chromosomes 1B, 3B, 5B and 7A based on the repeatability scores. Locus QSl-1B, including 11 significant SNPs in 15LS and 17LS, was located at 674.96–675.13 Mb on chromosome 1B (Fig. 2) and explained 16.38%–23.12% of the phenotypic variation. Locus QSl-3B, consisting of five SNPs, was associated with SL in 15HS and located at 62.13–62.14 Mb on chromosome 3B (Fig. 2) and explained 12.18%–12.82% of the phenotypic variation. On chromosome 5B, four SNPs located at 574.06–578.35 Mb (Fig. 2), named as Q-5B1, were associated with the relative value of SL (16HS/16LS) and explained 12.94%–14.29% of the phenotypic variation. Importantly, two SNPs of Q-5B1 were significantly associated with the relative values of BM, GN and YPP. Locus QSl-5B2, containing 21 SNPs, was located in the interval 583.48–584.89 Mb on chromosome 5B (Fig. 2). It was significantly associated with the relative value (16HS/16LS) of SL and could explain 13.03%–15.28% of the phenotypic variation (Table S6). Furthermore, one SNP of QSl-5B2 was significantly associated with the relative values (16HS/16LS) of BM (Table S6). Locus Q-7A, involving 12 SNPs, was mapped at 22.75–22.88Mb on chromosome 7A (Fig. 2), among which ten SNPs were significantly associated with SL in 15LS and 16LS, explaining 20.20%–28.47% of the spike length variation. In total, seven SNPs of Q-7A were significantly associated with YPP in 15LS and 16LS (Table S6). Additionally, two SNPs of Q-7A were significantly associated with GN in 15HS and 16HS (Table S6).
In total, 60 SNPs were significantly associated with GN and were mainly located on chromosomes 2A, 3A, 5B and 7A (Table S5). A comparison of the associated SNPs for SN and GN showed abundantly common SNPs from chromosome 5B that were accumulated at the common locus Q-5B3 (690.08–690.70 Mb) (Fig. 2). At this locus, 34 SNPs were detected in LS environments for GN, and five SNPs were detected under both LS and HS conditions. Notably, the correlation analysis of pairwise traits showed extremely high co-efficient of correlation between GN and SN, with values were 0.859 and 0.891 under LS and HS treatments respectively.
Up to 18 significant SNPs were associated with BM, and they were mainly distributed on chromosomes 4A, 4B, 5A and 5B (Table S5). Of these, nine and three SNPs on chromosome 5B were detected in 15LS and 16HS/16LS, respectively, and could explain 10.04%–13.90% of the phenotypic variation.
Yield per plant
A total of 48 significant SNPs were associated with YPP and were mainly distributed on chromosomes 2D, 3A, 5A, 5B and 7A (Table S5). Of them, six SNPs, forming locus QYpp-3A, were located on chromosome 3A at 694.01–694.96 Mb. The locus was associated with YPP in 16HS, 17LS and 17HS with a phenotypic variation of 11.36%–16.32% (Fig. 2). In addition, four SNPs, composing locus Q-5A2, were located on chromosome 5A at 692.16–692.39 Mb and detected in 15HS, with a phenotypic variation of 10.69%–12.70% (Fig. 2). The SNPs of Q-5A2 were significantly associated with BM in 15HS. The common locus associated with YPP and SL was Q-7A (22.75–22.88 Mb), which was observed on chromosome 7A (Fig. 2).
As many as 68 significant SNPs were associated with SN and were mostly located on chromosomes 1A, 1B, 2A, 2B, 5A, 5B, 5D, 6D and 7B (Table S5). Of these, 19 SNPs on chromosome 5B were detected in at least two environments under LS conditions, while two SNPs on chromosome 5A were detected in at least two environments under HS conditions. Locus Q-5B3 consisting of 36 SNPs was found on chromosome 5B at 690.08–690.70 Mb, and 24 SNPs in this locus were associated with SN in 15LS, 16LS and 17LS (Fig. 2) and could explain 17.29%–32.82% of the phenotypic variation. Moreover, the SNPs of Q-5B3 were almost significantly associated with not only GN in 15LS (34) and 16HS (5), but BM in 15LS (9) (Table S6).
Totally, 31 significant SNPs associated with TKW were mainly distributed on chromosomes 3B, 3D, 4A, 4B and 7B (Table S5). Among these, the largest number of SNPs (up to 25) was from chromosome 4A. These SNPs, forming locus QTkw-4A at 14.68–16.28 Mb, were detected in 15LS and 17LS and could explain 18.60%–35.57% of the phenotypic variation.
Phenotypic effects of alleles in the natural population
A t test was performed to verify the favorable alleles in the natural population based on the significantly associated SNPs, and some loci in this study were identified as possibly affecting the agronomic traits of wheat under salt-stress conditions. For instance, accessions with favorable alleles at loci QSl-1B (674.96–675.13 Mb), QSl-3B (62.13–62.14 Mb) and QSl-5B2 (583.48–584.89 Mb) showed longer SLs (increased by 0.54–0.97 cm with LS and by 0.57–2.44 cm with HS). Favorable alleles of locus QTkw-4A (14.68–16.28 Mb) increased the TKW by 4.06–5.32g with LS and by 2.49–5.95 g with HS. Favorable alleles of locus QPh-4D (16.93–19.29 Mb) increased PH by 3.42–4.57 cm with LS and by 4.21–4.31 cm with HS. Accessions with favorable alleles at locus Q-5B1 (574.06–578.35 Mb) have longer SLs (by 1.06–2.46 cm) and greater BMs (by 4.72–5.45 g), GNs (by 23.80–65.29) and YPPs (by 2.58–8.69 g) with HS. Accessions with favorable alleles at locus Q-5B3 (690.08–690.70 Mb) have greater BMs (by 1.49–4.54 g with LS and by 1.44–1.89 g with HS), SNs (by 0.46–1.50 with LS and by 0.35–0.64 with HS), GNs (by 14.49–60.75 with LS and by 16.16–25.68 with HS) and YPPs (by 2.71–5.73 g with LS and by 2.43–2.88 g with HS). Favorable alleles of locus Q-7A (22.75–22.88 Mb) increased SL and YPP by 0.71–1.28 cm and 2.16–3.39 g with LS, respectively, and increased SL by 1.07–1.11 cm with HS (Table S7).
The frequencies of the alleles of the above key loci ranged from 5.76% to 93.72% (Table S7). An analysis of the distributions of favorable alleles indicated that introduced varieties possessed the highest average number of favorable alleles, followed by improved cultivars. The average number of favorable alleles in landrace was lowest (Fig. S1a). Except accessions from Gansu and Sichuan, which were mostly landraces, accessions from Hebei, Henan and Shaanxi had comparatively higher average numbers of favorable alleles, followed by the cultivars from Shandong, Beijing and Shanxi. However, cultivars from Qinghai had the lowest average number of favorable alleles (Fig. S1b). Thus, it was believed that the introduced varieties and the improved cultivars should be utilized to improve these yield-related traits in salinity tolerant wheat cultivars. In addition, the associated traits of favorable alleles varied for some well-known cultivars. For example, although accessions Kenong 181, Liangxing 99, Zaoyou 504 and Jinhe 9123 had similar SL, Kenong 181 and Liangxing 99 contained the favorable alleles from loci QSl-1B, QSl-5B2 and Q-7A, whereas Zaoyou 504 and Jinhe 9123 had the favorable alleles from loci QSl-1B, QSl-3B1, and Q-7A. Besides, the YPP of landrace Chadianhong and Baiqimai was similar. Chadianhong contained the favorable alleles from loci Q-5A2 and Q-7A, while Baiqimai had the favorable alleles from loci QYpp-3A and Q-5A2.
Development of Kompetitive Allele Specific PCR (KASP) markers for the important loci
In total, 14 KASP markers were designed and used to detect polymorphisms in ZX-RIL and HL-DH populations (Table S9). Eight KASP markers were polymorphic. Of these, four (kasp-5695, kasp-9337, kasp-6491 and kasp-4516) showed polymorphisms in ZX-RIL population, two (kasp-9179 and kasp-4767) were polymorphic in HL-DH population, and two (kasp-8827 and kasp-9508) were polymorphic in both ZX-RIL and HL-DH populations. Based on the genotyping results of each KASP marker, a t test was conducted for corresponding traits in the related DH and RIL populations.
KASP marker kasp-9508, which came from locus Q-5A2, was polymorphic in ZX-RIL as well as HL-DH populations. In HL-DH population, kasp-9508 was significantly correlated with PH, BM and TKW in two environments (LS and HS) and was correlated with YPP under HS conditions (Fig. S2a–d). For the locus Q-7A, markers kasp-9179 and kasp-4767 were polymorphic in HL-DH population. Kasp-9179 was significantly correlated with PH in two environments (LS and HS), with SL under LS conditions, and with BM and YPP under HS conditions (Fig. S2e–h). Kasp-4767 was significantly correlated with PH, BM and YPP under HS conditions and with SL under LS conditions (Fig. S2i–l). Thus, loci Q-5A2 and Q-7A should be stable and pleiotropic loci affecting BM, PH, SL, TKW and YPP under salt-stress conditions.