Dissection of the genetic basis of genotype-by-environment interactions for grain yield and main agronomic traits in Iranian bread wheat landraces and cultivars

Understanding the genetic basis of performance stability is essential to maintain productivity, especially under severe conditions. In the present study, 268 Iranian bread wheat landraces and cultivars were evaluated in four well-watered and two rain-fed conditions for different traits. According to breeding programs, cultivars were in a group with a high mean and stability in terms of GY, GN, and SW traits, while in terms of PH, they had a low mean and high stability. The stability of cultivars and landraces was related to dynamic and static stability, respectively. The highest number of marker pairs and lowest LD decay distance in both cultivars and landraces was observed on the B genome. Population structure differentiated indigenous cultivars and landraces, and the GWAS results for each were almost different despite the commonalities. Chromosomes 1B, 3B, 7B, 2A, and 4A had markers with pleiotropic effects on the stability of different traits. Due to two rain-fed environments, the Gene Ontology (GO) confirmed the accuracy of the results. The identified markers in this study can be helpful in breeding high-performance and stable genotypes and future breeding programs such as fine mapping and cloning.


Results
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). Accord- Table 1. Mean, standard deviation (SD), broad sense heritability (H 2 ), and combined analysis of variance based on studied traits in 286 Iranian wheat landraces and cultivars and 6 environments. *, ** and *** are significant at the probability level of 5%, 1% and 0.1%, respectively. 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 r 2 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.

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 W i 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 HMR-PGV 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 b i 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). 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 W i . 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   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 W i index. Most b i -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 W i 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 b i , 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 biosyn- www.nature.com/scientificreports/ thesis of secondary metabolites, carotenoid biosynthesis, fatty acid elongation and ubiquinone, and another terpenoid-quinone biosynthesis (Table 5).

Discussion
The high significance of GEI for the studied traits was expected in this study, which accords with the previous reports 9 . Similar to this study, different monthly rainfall in MET studies, which also has drought stress, is one of the main reasons for GEI 10,33 . For some traits, the effect of GEI in cultivars was less than in landraces. This result is due to breeding programs and the small number of samples in the cultivars compared to the landraces, leading to fewer effects of GEI. Severe GEI caused low heritability of traits, especially in GY. In general, heritability and repeatability for complex traits such as GY are low compared to PH 9,28,34 . High yield and stability of wheat cultivars were expected since new wheat genotypes tolerate adverse environmental conditions such as drought www.nature.com/scientificreports/ stress 35,36 . On the other hand, breeding programs have improved wheat adaptation throughout a century 25 and continued to provide adapted wheat germplasm 37 . The genotypes in the fourth group in each of the traits, which included unstable low-yield genotypes, mainly consisted of landraces. Also, landraces had the highest percentage of genotypes selected by the multi-trait stability index. Most likely, the lack of specific selection for high yield and priority of yield stability during wheat domestication has led to such a result 38 . However, severe genetic heterogeneity in the Iranian wheat landraces and the application of some early breeding processes 39 have put landraces in desirable groups in terms of yield and stability. In this context, additional assessments with a large number of locations are needed to fully explain GEI patterns. The concepts of static and dynamic stability can be clearly distinguished based on bi and Wi indices in GY. The genotypes of group I, i.e., most cultivars, had dynamic stability. In contrast, although unstable in terms of dynamic concept, the fourth group, including the landraces, had static stability due to the low values for bi. The static concept is associated with low GY 8 . The genotypes of the second group, which included a small number of cultivars and landraces, were unstable in terms of both concepts despite their high yield and had good adaptability according to WAASBY and HMRPGV indices. We found a distinction between the concepts of stability and adaptability in other traits, especially PH. However, the literature paid scant attention to such a distinction between cultivars and landraces in terms of stability.
We found that the studied SNPs covered the wheat genome well. The number of SNPs based on the new wheat reference genome was higher in the B genome and lower in the D genome. There also seemed to be a direct relationship between marker density and chromosome size, and such a frequency of SNPs results from the evolutionary process of wheat. This conclusion was reported by Alipour et al. 39 in Chinese Spring and W7984 reference genomes. Other similar results were confirmed by Mourad et al. 40 and Edae et al. 41 . The difference of r 2 in cultivars, landraces, and different chromosomes, in addition to the evolutionary process, indicates the effect of breeding programs 42 . In this regard, comparing landraces and cultivars of wheat in China and Pakistan showed that the distances of LD decays in the landraces were less than cultivars. On the other hand, LD decays in genome A was slower than that of B 43 . Given that the landraces are genetically heterogeneous and are collected from areas with different climates, we expected that their heterozygosity would be high. Environmental factors affect genetic diversity and the structure pattern of plant populations 44 . Therefore, the high level of gene diversity in the studied population can be attributed to the geographical diversity of collection sites, differences in growth habit, etc. These factors led us to observe two subpopulations that separated cultivars well from the landraces. Moreover, the breeding programs and improved accessions are the reasons for such a separation. Iranian wheat genotypes have been categorized into two subpopulations in the previous studies 40 . The mean PIC value for all genotypes was 0.27, which is a good value for the bi-allelic marker 39,45 , and given their good distribution throughout the genome, they can be used to understand the genetic basis of GEI control. www.nature.com/scientificreports/ Genome-wide association studies capture the genetic loci linked to significant variation for traits of interest in a vast collection of wild relative populations, breeding cultivars, and landraces 46,47 . It is also an important tool for selecting high-yield genotypes in a group of environments 33   www.nature.com/scientificreports/ including those that were not mapped to any chromosome. The number of MTAs identified in all genotypes was higher for GY and its indices in the D genome compared to the B genome, and for PH, it was higher than in the A and B genomes. In a study, GWAS for GY was run in each environment due to the presence of GEI, and the results showed that the D genome had the highest number of SNPs 33 . This suggests that the role of the D genome in wheat adaptability should be further addressed 48 . The greatest number of significant MTAs were identified on chromosome 6B in both cultivars and landraces datasets, while the least numbers were detected on chromosomes 5D and 4D in cultivars and landraces datasets, respectively. Acuña-Galindo et al. 49 also found two meta-QTL for adaptation to drought stress on chromosome 6B in wheat. A recent study also reported a major grain yield QTL on chromosome 6B and fifteen haplotype blocks associated with two stability indices, including Lin and Binn's superiority index and Eberhart and Russell's coefficient on chromosomes 1A, 4A, 4B, 5B, 6B, 7A, 7B, and 7D 29 . In addition, genomic regions associated with grain yield and yield stability on chromosomes 2B, 3A, 4A, 5B, 7A,  www.nature.com/scientificreports/ and 7B were identified in CIMMYT's spring bread wheat 50 . Considering all genotypes, we located about 44% and 24% of the markers associated with the mean GY on chromosomes 6A and 7D, respectively. Interestingly, GO results showed that one of these markers is in the coding region of proteins that regulate transcription. Previous reports indicate that chromosome 6A contains GY and TGW-related locus in MET data that harbored a TaGW2-6A gene and that other genes influence its expression 51 . Chromosome 7D is of great importance in explaining GY phenotypic variation 28 . Muhu-Din Ahmed et al. 52 identified MTAs for GY on chromosomes 1A, 3A, 4A, 1B, 4B, 6B, 7B, 5D, and 7D under both well-watered and water-deficit conditions. Several studies also demonstrated MTAs for GY in various wheat panels analyzed thorough GWAS on chromosomes 2B, 3A, 3D, 5B, 7A and 7B 53 , 1A, 2D, 3A, 7B and 7D 1 , and 1B 54 under different water regimes. The marker locus on 4B in GY under water stress conditions was also associated with this trait in the Pakistani wheat population 55 . Similarly, in genome-wide association mapping, Edae et al. 56 reported MTAs for GY on chromosomes 4A, 1B, 5B, and 2B of spring wheat association panel under contrasting moisture regimes. Moreover, Lozada et al. 57 found MTAs for GY on chromosomes 5A, 1B, 2B, and 4B in a diverse panel of 239 wheat genotypes evaluated across two growing seasons using SNP markers. Tadesse et al. 58 reported GY-related MTAs on 1B in 120 elite hexaploid wheat genotypes, which were evaluated under rain-fed and irrigated conditions for a genome-wide study.
The multi-trait loci controlling performance and stability were located on chromosomes 1B, 3B, and 7B. Furthermore, chromosomes 2A and 4A in all three cultivars, landraces, and the sum of these two had multi-trait control loci. All chromosomes, except for chromosome 3B, were reported in a similar study 9 . In another study, chromosomes 3B and 2B, 3A, 4A, 5B, 7A, and 7B were associated with wheat yield stability coefficient 50 . Major QTLs with pleiotropic effects on chromosomes 3B and 7B have also been confirmed 59 . One study concluded that a specific combination of photoperiod genes increases the yield stability of durum wheat 14 . Also, the best allelic combination using stepwise regression in markers identified by genome-wide association mapping (GWAM) can lead to increased stability and yield in wheat 50 . Therefore, it is possible to say that yield stability is controlled by genes with pleiotropic effects. However, as the experiment was performed in the same place and under different conditions, the correlation between grain yield in different environments may be a reason to observe common SNPs. In this regard, the lack of correlation between the environments resulted in no common SNP for the GWAS performed in 9 environments 33 . Although several common MTAs were identified in for GY, GN, SW, and PH traits and different stability indices, these traits are not exclusive and independent. Thus, it is possible to select both traits and stability indices in Iranian wheat cultivars and landraces since most significant MTAs (almost 90%) were not common among the trait values and stability indices. Lozada and Carter 9 identified 12 SNP loci linked to both trait value and stability parameters in Pacific Northwest winter wheat. Two major effect SNP markers of Tdurum_contig61410_542 (1B) and BS00022542_51 (7B), were associated with grain yield and yield stability indices. The common MTAs between different traits and yield stability coefficient have already been reported 50 . The low number of MTAs identified for PH is probably due to the fact that this trait is controlled by a small number of genes compared to other traits. However, the above results for yield and its components show that they are controlled by several genes that interact with each other and the environment. Stability-associated genes can also be stress-responsive genes 50 . Therefore, GO results could be well described, given that two of the six environments are under rain-fed conditions. Proteins phosphorylation, especially in wheat grains, play an important role in drought stress 60 . Jasmonic acid biosynthetic modulates drought stress in wheat 61 . Markers related to mean GY and SW were annotated with antioxidant activity. Reducing the effects of drought stress by such activity with various enzymes in wheat was demonstrated by previous researchers 62 . The Synthesis of fatty  www.nature.com/scientificreports/ acids is useful in counteracting the drought stress in oats 63 . Transmembrane transport, DNA-binding transcription factor activity, DNA repair, and peptidase activity were other examples that were annotated and possibly involved in response to drought stress. These results are similar to the previous reports 64 . Earlier efforts have been made to interpret GWAS results and understand GEI using gene annotation 33 . KOBAS is a useful tool for genome annotation 65 . It has been shown that ubiquinone and other terpenoid-quinone biosynthesis are metabolic pathways of response to drought stress in plants 66 . In addition, carotenoid biosynthesis is involved as one of the KEGG pathways in drought stress tolerance 67 . Such an important role for the biosynthesis of secondary metabolites has been proven 68 .

Conclusions
In the current study, GWAS was performed for some important agronomic traits and different static and dynamic stability indices based on those traits were calculated in a diverse panel of 268 Iranian wheat cultivars and landraces. The highest number of marker pairs and lowest LD decay distance in both cultivars and landraces was observed on the B genome, whereas the D genome had the least number of marker pairs and most significant LD decay distance. A total of 846, 653, and 1023 significant MTAs were identified for the traits and their related stability indices in cultivars, landraces, and total genotypes datasets, respectively. The chromosomes 6B and 4D had the highest and lowest number of MTAs, respectively. The multi-trait loci controlling mean traits and stability were located on chromosomes 1B, 3B, and 7B, and GO results for highly significant MTAs almost confirmed the accuracy of the identified markers. The identified markers in this study could provide valuable genetic resources to initiate marker-assisted selection, fine mapping, and cloning the underlying genes and QTLs.

Methods
Plant materials and field evaluation. A set of 268 Iranian bread wheat genotypes, including 180 landraces and 88 cultivars, were studied in six environments (Supplementary Table 1 Table 2). Trials were planted in early November and harvested in July of the next year. The experiments were performed on the research farm of the University of Tehran with latitudes of 50.58 E and 35.56 N and 1112.5 m above sea level in a randomized complete block design with two replications. The dimensions of the plots consisted of four lines with a length of 1 m (80 × 100 cm). The distance was 20 and 5 cm between and within the rows. Plant height (PH, cm), grain number per spike (GN), spike weight (SW, g), grain yield per plant (GY, g plant -1 ) were traits that were measured based on ten randomly selected samples from each plot. Plant height was recorded from ground level to tip of the spike, excluding awns, at maturity stage. After harvesting, all spikes were hand-threshed to determine the GY, SW, and GN. Then, stability parameters (Table 6) of each trait were calculated using 'agricolae' 69 , 'ammistability' 18 , and 'metan' 70 packages in the R and STABILI-TYSOFT online programs 71 . Broad sense heritability of traits was calculated using the following equation: where σ 2 g and σ 2 ge are the variance due to genotype, and genotype-by-environment interaction, respectively. σ 2 ε is the residual variance, and e and r are the number of environments and replications, respectively 72 .
Genotyping. The development and genetic material studied was previously described based on genotyping by sequencing of a GBS library for the Iranian wheat samples have been by Alipour et al. 39 . In brief, sequence reads were first trimmed to 64 bp and were grouped into sequence tags. Then, SNPs were identified using internal alignment allowing for mismatch up to 3 bp. The UNEAK (Universal Network-Enabled Analysis Kit) GBS pipeline was used for SNPs calling, where reads with a low-quality score (< 15) were discarded. Imputation was performed in BEAGLE v3.3.2 73 using w7984 reference genome 74 . Finally, SNPs with heterozygotes < 10%, and minor allele frequency > 5% were used for further analysis.

Genome-wide association study. Both general linear model (GLM) and mixed linear model (MLM)
were employed to obtain the unbiased estimation of marker effects using TASSEL 5.0 75 software and GAPIT R-package 76 . The results of GLM was adjusted using the first three principal components (PCA) and population structure (Q) and MLM was corrected using kinship-matrix with the first three principal components (PCA + K) and population structure (Q + K). Results of all approaches from both TASSEL and GAPIT were evaluated based on the Q-Q plot and significance of associated loci using t-tests. In general, the results of the MLM approach of the first three principal components and kinship-matrix (PCA + K) obtained from GAPIT provided a more robust control of confounding effects. We, therefore, only reported the results MLM obtained from GAPIT. In the MLM model, individuals are considered random effects, and the relatedness among individuals is conveyed through a kinship matrix. A threshold of -log 10 (p) > 3 was used to state statistically significant MTAs 77,78 . Confidence intervals (CIs) for MTAs were calculated for each chromosome using the linkage disequilibrium (LD) decay. Circular Manhattan plots were performed using the CMplot R-package 79 . Gene annotation. Sequences surrounding all significantly associated SNPs were obtained from the blast tools in EnsemblPlants database (http:// plants. ensem bl. org/ index. html) to assess gene annotation using Gramene (http:// www. grame ne. org/) by aligning them to the IWGSC RefSeq v1.0 annotation (https:// wheaturgi. versa illes. inra. fr/ Seq-Repos itory/ Annot ations). After aligning SNPs sequences to the reference genome, we selected overlapping genes with the highest identity percentage and blast score for further processing. The gene ontology of each selected gene, including molecular function and biological process, was extracted from the H 2 = σ 2 g /(σ 2 g + (σ 2 ge /e) + (σ 2 ε /er)) Table 6. Description of the stability statistics studied. g: number of genotypes, e: number of environments, X ij : average yield of genotype i in environment j, X i· : average yield of the genotype i, X ·j : average yield of the environment j, X ·· : grand average yield, Ij: Environmental index, which is the deviation of the average of all genotypes in a specific place from the total average, PC: interaction principal components, SSPC n : sum of squares of the nth PC, GV ij : genotypic value of genotype i in environment j, u j : represents the mean of environment j, g i and ge ij are the BLUP values of genotype i and the interaction between genotype i and or environment j, respectively, θ n : the percentage sum of squares explained by the nth principal component interaction effect, RASV: rank of AMMI stability value, RY: rank of the mean yield of genotypes across environments, μ j : general mean for each environment j, IPCA gn is the score of the genotype g in the nth IPCA, and EP n is the amount of the variance explained by the nth IPCA. θ Y and θ S are the weights for yield and stability, respectively. rG g and rW g are the rescaled values of the gth genotype for yield and WAASB, respectively, G g and W g are the yield and WAASB values for gth genotype, respectively.