Population structure and LD analysis
A total of 120 cultivars and 8,632 SNPs with known physical positions were retained for further analysis. The number of SNPs within 1 MB across each chromosome is shown in Figure S1. Moreover, the number of SNPs across each chromosome is shown in Table S1. These SNPs were distributed on all 21 wheat chromosomes, with 3,260, 4,310, and 1,062 SNPs in the A, B, and D subgenomes, respectively.
The population structure of the natural population was based on unlinked SNP markers (Cormier et al. 2014). According to the trend of the CV error value curve from the results of the admixture software, the K value may be four or five (Fig. 2a). PCA classified the population into four major groups, although there were many admixtures (Fig. 2b). The neighbor-joining tree divided the population into four major groups (Fig. 2c). According to the planting area of each cultivar and the four subgroups, blue parts (A) of the population represent most of the Shandong area cultivars registered in the 1990s, purple parts (B) of the population represent most of the Shandong area cultivars registered in the 2000s, red parts (C) of the population represent most of the Hebei area cultivars, and green parts (D) of the population represent most of the Henan area cultivars. Most of the cultivars from the same planting area clustered together, which inferred that cultivars from the same province were more closely related than those from different provinces (Fig. 2d). Overall, based on the above three methods and the planting area of each cultivar, it was appropriate to classify the association population into four subgroups.
LD was estimated by calculating the r2 among all possible pairs of markers for each of the 21 chromosomes. The obtained r2 values were then plotted against the physical distance for each of the three genomes separately (Figure S2) and across the whole genome. When the cutoff threshold for r2 was defined as 0.1, the LD decay distance for the entire genome was greater than 15 MB.
Phenotypic trait analysis among subpopulations under different conditions
We analyzed each phenotypic trait among subpopulations under different conditions. The performance of all traits in the LN and CK environments is displayed in Fig. 3. From Fig. 3, we found that the same traits of subgroups were different. By comparing the same traits under the CK and LN environments within a subgroup (Table S2), we found that the PH and PH components were both greater in the CK treatment than in the LN treatment, except for the SECILT and THILT of the C subgroup. In the ANOVA results, a significant P value of σe2 existed among the different environments (Table S3), which illustrated that the environment affected the traits. According to the comparison of all traits in the CK and LN environments, the PH, SL, FORILT, FIFILT, and LRS showed significant differences (P < 0.001). The results indicated that the decrease in PH in the LN treatment was mainly caused by decreases in the length of FORILT, FIFILT, and SL in the LN treatment.
Using the ANOVA results and formula of broad sense heritability, we calculated h2 for each trait. The h2 of the tested traits ranged from 0.668 (SIXILT) to 0.97 (PH). For four traits (PH, SECILT, THILT, FORILT), h2 was greater than 0.9, whereas h2 for FIRILT, LRS, PL, FIFILT, and SL was 0.884, 0.875, 0.862, 0.810, 0.799 and 0.668, respectively (Table S3).
The average data based on each trait were used to calculate Pearson's correlation coefficients. PH was positively correlated with SECITL, THILT, and FORILT (Figure S3); SL was negatively correlated with FORILT and FIFILT (Figure S3). PL showed a strong positive correlation with FIRILT (Figure S3), and the LRS showed no correlation with other traits (Figure S3), which indicated that LRS had no significant relationship with the length of internodes (Figure S3).
GWAS of plant height-related traits under low nitrogen and normal environments
To eliminate the influence of the kinship and structure of the population, we performed GWAS for all traits with the MLM method under CK and LN environments, and significant SNPs (P < 1e-3), repetitive SNPs in multiple environments (all, E > 3; CK, E = 2; LN, E = 2) and SNPs located within the LD decay distance (15 M) were defined as repetitive QTLs. Manhattan plots and QQ plots for all traits are shown in Fig. 4 and Figure S4. All detected SNPs and QTLs in the CK and LN environments are displayed in Table S4. GWAS analysis showed that a total of 86 QTLs were detected in the LN and CK environments.
Plant height and plant height components
A total of 12 QTLs were detected after analysis of PH (Fig. 4; Table S4), including 7 QTLs in the CK environment and 6 QTLs in the LN environment, which were mainly distributed on chromosomes 1B, 3A, 3B, 3D, 4B, 6B, 7A, and 7B. Two QTLs, QPH.sdau-3A.1 (142 MB) and QPH.sdau-4B (555–570 MB), were detected in both the CK and LN environments, explaining 8.0-10.6% and 9.2–10% of the phenotypic variation, respectively (Fig. 4; Table S4, PH). Two QTLs, QPH.sdau-3B.2 (814 MB) and QPH.sdau-7A.1 (34–49 MB), were detected in the 2017CK and 2018CK environments and explained 9.5–10% and 8.0-8.3% of the phenotypic variation, respectively. No repetitive QTLs were detected with PH in the 2017LN and 2018LN environments.
A total of 10 QTLs mainly distributed on chromosomes 2A, 2B, 3A, 3B, 3D, 5A, 5B, and 7B were associated with SL, including 6 QTLs in the CK environment and 4 QTLs in the LN environment (Fig. 4; Table S4, SL). Among all QTLs, QSL.sdau-3D (602 MB) in 3D was detected in the 2017CK and 2018CK environments, while no SNP was detected in the LN environment. The phenotypic variation explained by QSL.sdau-3D ranged from 7.9–11%. For FIRILT, 19 QTLs were identified on chromosomes 1A, 1B, 1D, 2A, 3A, 3B, 4B, 5A, 6A, 6B, 7A and 7B, including 14 QTLs in the CK environment and 8 QTLs in the LN environment (Fig. 4; Table S4, FIRILT). QFIRILT.sdau-3A.2 (503 MB) and QFIRILT.sdau-4B (555–570 MB) were both detected in the 2017CK, 2018CK, 2017LN, and 2018LN environments, explaining 9.8–12.8% and 9.9–13.7% of the phenotypic variation, respectively. QFIRILT.sdau-7B.2 (606–608 MB) was detected in the 2017CK and 2018CK environments and explained 10.9–12.3% of the phenotypic variation. A total of 6 QTLs for SECITL were identified on chromosomes 1B, 3A, 4A, 4B, 6B, and 7B, including 2 QTLs in the CK environment and 4 QTLs in the LN environment (Fig. 4; Table S4, SECITL). Among all QTLs, the QTL QSECILT.sdau-4B (535–570 MB) was detected in the 2017LN and 2018LN environments. The phenotypic variation of QSECILT.sdau-4B ranged from 8.6–9.9%. For THILT, 11 QTLs were identified on chromosomes 2A, 2D, 3A, 3B, 4B, 6B, 7B and 7D, including 7 QTLs in the CK environment and 6 QTLs in the LN environment (Fig. 4; Table S4, THITL). A total of 13 QTLs for FORILT were identified on chromosomes 1B, 3A, 3B, 5B, 6B, 7A, 7B, and 7D, including 7 QTLs in the CK environment and 7 QTLs in the LN environment (Table S4, FORILT). For FIFILT, 14 QTLs were identified on chromosomes 2A, 3A, 3B, 5A, 5B, 6A, 6B, 6D, 7B and 7D, including 11 QTLs in the CK environment and 3 QTLs in the LN environment (Fig. 4; Table S4, FIFILT). From SECITL to FIFILT, no repetitive QTLs were detected in the CK or LN environments. Plant height in wheat is the spike length plus the length of all internodes above the ground, while SIXILT existed in a small number of cultivars; therefore, so it was not suitable to perform GWAS with SIXILT in the present study.
Plant height-related traits
Although PL and LRS are not components of PH, they are associated with plant morphology, which is closely related to PH. A total of 13 QTLs were detected with PL on chromosomes 2A, 3A, 3B, 4A, 4B, 5B, 6A and 7B, including 7 QTLs in the CK environment and 9 QTLs in the LN environment (Fig. 4; Table S4, PL); however, none of the QTLs were repetitively deactivated. A total of 7 QTLs for LRS were detected on chromosomes 2B, 3A, 3B, 4A, and 6A, including 5 QTLs in the CK environment and 3 QTLs in the LN environment (Fig. 4; Table S4, LRS). Among them, QLRS.sdau-4A (630 MB) was detected in the 2017CK, 2018CK, 2017LN, and 2018LN environments, explaining 11.2–12.5% of the phenotypic variation. QLRS.sdau-2B (42 MB) and QLRS.sdau-3B (657–664) were detected in the 2017LN and 2018LN environments, while QLRS.sdau-3A.1 (487 MB), QLRS.sdau-3A.3 (513 MB) and QLRS.sdau-6A (614 MB) were detected in the 2017CK and 2018CK environments. QLRS.sdau-2B and QLRS.sdau-3B explained 12.1–12.5% and 12.1–16.8% of the phenotypic variation, respectively (Table S4, LRS). QLRS.sdau-3A.1, QLRS.sdau-3A.3 and QLRS.sdau-6A explained phenotypic variations of 10.8–11.4%, 11.0-11.3%, and 11.5–12.1%, respectively (Table S4, LRS).
GWAS of plant height-related traits in response to nitrogen stress
According to the comparison of the same traits in the CK and LN environments, different nitrogen contents in the soil caused the difference in FORILT, FIVILT, SL, PH, and LRS. Tdelt_aver was measured as the trait that responded to nitrogen stress. Based on GWAS and Tdelt_aver, QTLs that responded to nitrogen stress were detected (Fig. 5; Table 1). A total of 18 QTLs that responded to nitrogen stress were detected.
In total, three QTLs for Tdelt_aver_FORILT responding to nitrogen stress were detected on 3A, 6A, and 6D (Fig. 5; Table 1). Qdelt_aver_FORILT.sdau-3A (657 MB), Qdelt_aver_FORILT.sdau-6A (18 MB) and Qdelt_aver_FORILT.sdau-6D (17 MB) explained 12.1–12.5% and 12.1–16.8% of the phenotypic variation, respectively (Table 1). A total of 11 QTLs detected for Tdelt_aver_FIFILT in response to nitrogen stress were detected on 1A, 1B, 1D, 2A, 2B, 3A, 3B, 4A, 5B and 7B (Fig. 5; Table 1). Among them, Qdelt_aver_FIFILT.sdau-3A.2 (648–658 MB) had the lowest P-value and explained 36.3% of the phenotypic variation, which was considered a main-effect QTL (Table 1). Qdelt_aver_SL.sdau-3A (705 MB) was detected for Tdelt_aver_SL in response to nitrogen stress, which explained 11.3% of the phenotypic variation (Fig. 5; Table 1). Qdelt_aver_PH-5B (701 MB) was detected for Tdelt_aver_PH in response to nitrogen stress, which explained 9.5% of the phenotypic variation (Fig. 5; Table 1). Two QTLs, Qdelt_aver_LRS.sdau-3B (150 MB) and Qdelt_aver_LRS.sdau-5B (520 MB), for Tdelt_aver_LRS responded to nitrogen stress and explained 12.0% and 10.4% of the phenotypic variation, respectively (Fig. 5; Table 1).