Genotype-by-environment interaction
The effects of genotype, environment, and GEI were significant at different probability levels for the four traits in the total population and subpopulations (Table 1). Due to drought stress in the study and different rainfall patterns in different years (Fig. 1), such a result was not unexpected. Broad sense heritability was low for GY, moderate for SW and GN, but high for PH. To select the desired genotypes in terms of mean traits and stability, we used different statistics, and the results are presented in Fig. 2. Based on these two criteria, genotypes were divided into approximately four classes: (I) high mean and stable, (II) high mean and unstable, (III) low mean and stable, and (IV) low mean and unstable. In terms of GY, 54, 29, 70, and 115 genotypes were present in these classes, respectively. Cultivars included 37.5%, 9.1%, 36.4%, and 17%, and landraces included 11.7%, 11.7%, 21.1% and 55.5% of the members of these classes, respectively (Fig. 2A). In terms of GN, in class I 43 genotypes (38.6% of cultivars and 5% of landraces), in class II 113 genotypes (42% of cultivars and 42.2% of landraces), in class III 99 genotype (14.8% of cultivars and 47.8% of landraces), and in class IV 13 genotypes (4.5% of cultivars and 5% of landraces) were present (Fig. 2B). According to SW, 135 genotypes (79.5% of cultivars and 36.1% of landraces), 13 genotypes (3.4% of cultivars and 5.5% of landraces), 45 genotypes (6.8% of cultivars and 21.7% of landraces), and 66 genotypes (10.2% of cultivars and 36.7% of landraces), were observed in the mentioned four classes, respectively (Fig. 2C). For plant height, there were 85 genotypes (19.3% of cultivars and 37.8% of landraces) in class I, 40 genotypes (10.2% of cultivars and 17.2% of landraces) in class II, 97 genotypes (62.5% of cultivars and 23.3% of landraces) in class III, and 46 genotypes (8% of cultivars and 21.7% of landraces) in class IV (Fig. 2D). On the other hand, as expected, some indices, especially HMRPGV and WAASBY, were correlated with the mean of the traits and were in the same group (Fig. 2). Since it is difficult to select genotypes with simultaneous stability for all four traits, we calculated the multi-trait stability index based on yield and yield components (Supplementary Fig. 1). The results interestingly showed that 11 cultivars (12.5%) and 29 landraces (16.1%) formed the genotype selected based on this index (Supplementary Table 3).
Average rainfall in different months each year.
Heatmap based on stability indicators for grain yield (A), number of grains (B), spike weight (C), and plant height (D) in Iranian wheat landraces and cultivars. I: high mean and stable, II: high mean and unstable, III: low mean and stable and IV: low mean and unstable. Mean: average trait in all environments, ASV: AMMI stability value, bi: Finlay-Wilkinson regression, Wi: Wricke’s ecovalance measures, MASI: Modified AMMI stability index, YSI: Yield stability index, HMRPGV: harmonic mean of the relative performance of the genetic values, WAASBY: weighted average of absolute scores from the singular value decomposition of the matrix of BLUP for the GEI effects generated by an LMM and response variable. The heatmaps were created using “gplots” package “heatmap.2” function in R32.
Genetic data and population structure
Based on the results, the distribution of SNP markers showed that genome B alone accounted for 50% of the total markers, while genome D had the lowest number of SNP markers by far. In the A, B, and D genomes, chromosomes 7A, 3B, and 2D, respectively, had the highest number of SNPs in cultivars, landraces, and the sum of the two (Table 2). The density (SNP/Mbp) was similar, with the B genome having the highest density of SNPs, especially for chromosomes 6B and 3B. This is more conveniently illustrated in Fig. 3. The average minor allele frequency (MAF) and gene diversity (GD) in cultivars were slightly higher than the landraces. The amount of heterozygosity (HET) of the landraces in each of the chromosomes and consequently the genomes were higher than the cultivars. The polymorphism information content (PIC) in cultivars ranged from 0.240 (4D) to 0.309 (2A) and in landraces from 0.232 (2D) to 0.292 (4A). The mean PIC in cultivars, landraces, and the sum of these two was equal to 0.280, 0.267, and 0.270, respectively (Table 2). On the other hand, the total number of SNP pairs (TNSP) and the number of significant SNP pairs (NSSP) were higher in the B genome (especially on chromosomes 3B, 2B, and 6B) and lower in the D genome (especially on chromosomes 4D, 5D, and 3D). The percentage of NSSP in cultivars ranged from 25.11% (4D) to 58.26% (4A) and in landraces ranged from 26.16% (4B) to 53.27% (4A). The r2 values of cultivars were higher than landraces, especially in B and D genomes. Such a difference in distance (cM) can also be seen in the D genome (Table 3). The results of genetic population structure analysis indicated the existence of two subpopulations (Fig. 4A). The highest value of ΔK was observed at K = 2 (Fig. 4B), and its average log-likelihood value confirmed it (Fig. 4C). One of these subpopulations consisted mainly of cultivars, and the other contained landraces.
Density plot by different chromosomes in total Iranian bread wheat cultivars and landraces.
Barplot (A), the average log-likelihood value (B), and delta K for different numbers of sub-populations (C), in the analysis of population structure using 43,446 SNP markers.
MTAs for mean traits and stability indices
An overview and detailed information of MTAs results are provided in Supplementary Tables 4 and 5. A total of 846, 653, and 1023 significant MTAs were identified for the studied traits and stability indices of cultivars, landraces, and total genotypes, respectively (Fig. 5). Circular Manhattan plots for common regions associated with different traits are plotted (Fig. 6). Ten and 12 markers were related to the mean grain yield of cultivars and landraces, respectively, mainly located in genome A. This number was higher with 55 markers for all genotypes. ASV and MASI statistics had the highest MTAs in the B genome, while for Wi in the D genome (especially chromosome 7D) and the B genome, most markers in landraces were identified on the B and A genomes. There were 22 and 14 significant associations for HMRPGV in cultivars and landraces, respectively. Chromosomes 4A and 2A for cultivars and 3D for landraces were important. WAASBY was significantly associated with 24 and 21 SNPs in cultivars and landraces. These markers were mainly distributed on chromosomes 6B, 2B, and 2D. Although bi for cultivars and landraces had the lowest MTAs in the D genome, this genome (especially its 6D chromosome) contained the highest MTAs considering the total genotypes. Finally, among all the indices, YSI in the cultivars was associated with the highest number of SNPs in the B genome (Fig. 5A).
GWAS results stability indicators for grain yield (A), number of grains (B), spike weight (C), and plant height (D) in Iranian wheat landraces and cultivars.
Circular Manhattan plots to draw common regions associated with grain yield (A), number of grains (B), spike weight (C), and plant height (D) in Iranian wheat landraces and cultivars. Inner to outer circles represents average trait and stability indices including ASV, bi, HMRPGV, MASI WAASBY, Wi, and YSI, respectively. The chromosomes are plotted at the outmost circle where thin dotted blue and red lines indicate significant level at p value < 0.001 (− log10 (p) > 3) and < 0.00001 (− log10 (p) > 5), respectively. Green and red dots indicate genome-wide significantly associated SNPs at p value < 0.001 and < 0.00001 probability level, respectively. Scale between ChrUn and Chr1A indicates − log10 (p) values. Colored boxes outside on the top right side indicate SNP density across the genome where green to red indicates less dense to dense.
For GN, 14 MTAs for the mean and 209 MTAs for the stability parameters were identified in the cultivars, compared to 24 and 171 MTAs for the landraces, respectively. Like GY, more MTAs were identified based on all genotypes. Chromosomes 6B and 2B in cultivars contained the highest markers associated with ASV and MASI, while the SNPs identified for these two indices were low in landraces and scattered on different chromosomes. Genome B, especially chromosome 2B, had the highest QTLs associated with Wi. In total, 14 and 30 MTAs were determined for HMRPGV in cultivars and landraces, respectively, with chromosomes 1A, 4A, and 5A having the highest SNPs in the landraces. The highest number of SNPs associated with WAASBY in cultivars and landraces were located in B and A genome, respectively. The highest number of bi and YSI-related SNPs belonged to the B genome, and chromosomes 6B and 7A in landraces were important for the YSI index (Fig. 5B).
Mean SW in cultivars, landraces, and total genotypes was identified as 31, 19, and 57 MTAs, respectively, mainly located on chromosomes 6A, 3B, 4B, 6B, and 7D. The B genome and then the A genome had a highly significant number of HMRPGV-related SNPs. For WAASBY in the cultivars, 16 and 3 MTAs were identified in the B and A genome, respectively. Although no significant MTAs were observed in genome D cultivars, like the A genome, it contained SNPs related to the WAASBY index in the landraces. In terms of ASV and MASI indices, 6.4 and 2.6 times more MTAs were detected in the cultivars compared to landraces, respectively. More MTAs were observed on chromosomes 6B and 6D in cultivars and 1B, 3D, 6A, 6B, and 7B in landraces for Wi index. Most bi-related SNPs were located on chromosomes 1A, 6B, and 7A in the cultivars and on 1D, 3B, 6B, and 7A in the landraces. Finally, for the SW, like GY, the highest number of MTAs we could see in a genome was the YSI index (Fig. 5C).
Among the traits, the lowest MTAs were observed for PH and its stability indices. Moreover, 13 markers on 1D, 2B, 3B, 5B, 7A, 7B, and 7D chromosomes in cultivars and seven markers on 1A, 1D, 2A, 5B, and 7B chromosomes were associated with mean trait. The markers identified for ASV and MASI were the same in the cultivars and slightly different in the landraces. Such similarity was observed by considering the sum of genotypes, with 3D and 7A chromosomes having a larger number of SNPs. Although genome B had the lowest number of MTAs for Wi in cultivars, it showed the highest association in landraces. Chromosomes 3A, 6D, and 7A in cultivars and chromosomes 4A and 6B in landraces were important for this index. For HMRPGV in cultivars, seven markers were identified on chromosomes 6D, 7A, 3D, and 5B. These numbers were equal to 12 and were distributed on chromosomes 7B, 5D, 5B, 1D, 1A, 6A, and 2B. According to WAASBY, 14 SNPs were identified in cultivars on different chromosomes, including 1A, 1B, 2B, 3B, 3D, 4B, 4D, 5B, 6A, and 7B. In the landraces, 10 SNPs were identified, more than half of which were located on chromosome 1D. The bi, in cultivars on chromosomes 6B and in the landraces on chromosomes 3D and 3D, had the highest number of MTAs. Finally, the number of MTAs detected for YSI in cultivars was three times higher than in landraces (Fig. 5D).
Among the identified markers, 171, 131, and 224 cases in cultivars, landraces, and the sum of these two overlapped with different traits and indices, respectively (Supplementary Table 5). For example, the marker rs65138 in cultivars and rs51479 in total genotypes were associated with the mean of three traits GY, GN, SW, and some of their stability indices and were located on chromosomes 1B and 3B, respectively. One such marker in the landraces was rs58587, which was located on chromosome 7B and was associated only with the stability indices of GY, SW, and PH. Other SNPs with many pleiotropy effects were located on chromosomes 6B, 2A, 2B, 4D, 3B, and 4A in the cultivars. These cases in the landraces included 1A, 2A, 4A, 2D, 6A, 3D, 1B, and 7D. Considering the total genotypes, we found that the SNPs associated with most of the traits and indices were on chromosomes 4B, 4A, 2A, 2B, 7D, 2B, 6A, and 5D (Supplementary Table 5).
Gene ontology
For a closer look, we studied the ontology of highly significant markers (P < 0.0001). Except for PH, some of the identified MTAs were involved in important biological and molecular processes for all traits. These genes were distributed on different chromosomes, including 1A, 1B, 1D, 2D, 3A, 4A, 4B, 6A, 6B, and 7A, with chromosome 4B, 1B, and 7A having the highest number (Table 4). Genes with MTAs mainly encoded proteins wrapped in biological and molecular processes associated with adaptation, including drought stress tolerance. Oxidoreductase activity, DNA-binding transcription factor activity, ATPase-coupled transmembrane transporter activity, protein kinase activity, protein binding, and integral component of the membrane were some of the molecular processes. Some biological processes also included the oxidation–reduction process, regulation of transcription, jasmonic acid biosynthetic process, transmembrane transport, protein phosphorylation, fatty acid biosynthetic process, and DNA repair. The KEGG orthology system was also used to accurately annotate the identified SNPs. The results showed that genes were involved in various pathways such as biosynthesis of secondary metabolites, carotenoid biosynthesis, fatty acid elongation and ubiquinone, and another terpenoid-quinone biosynthesis (Table 5).
