Dissection of genotype-by-environment interaction and yield stability analysis in barley using AMMI model and stability statistics

Barley is one of the most important cereal crops with considerable tolerance to various environmental stresses, which can maintain its productivity well in marginal croplands. The selection of stable and high-yielding barley genotypes and ideal discriminative locations is an important strategy for the development of new cultivars in tropical climates. Different statistical methods have been developed to dissect the genotype-by-environment interaction effect and investigate the stability of genotypes and select discriminative environments. The main objective of the present study was to identify high-yielding and stable barley genotypes and testing environments located in the tropical regions of Iran using 23 parametric and nonparametric stability statistics. In the present study, the grain yield stability in nineteen barley genotypes was investigated across five different locations over two consecutive years (2018–2020). The additive main effects multiplicative interaction (AMMI) analysis showed that environments (E), genotypes (G) and GE interaction effects were significant for grain yield. Using Spearman’s rank correlation analysis, a pattern map developed simultaneously for assessing relationships between grain yield and stability statistics and clustering of them, which allowed identifying two main groups based on their stability concepts. The biplot rendered using the weighted average of absolute scores (WAASB) and mean grain yield identified superior genotypes in terms of performance and stability. Among test environments, Darab, Gonbad and Zabol showed a high discriminating ability and played the highest contribution in creating GEI. Hence, these regions are suggested as discriminative sites in Iran for the selection of high-yielding and stable barley genotypes. As a conclusion from this research, all stability statistics together identify G10 and G12 as the superior barley genotypes; these genotypes could be released as commercial cultivars in tropical regions of Iran.

Page 2 of 12 Pour-Aboughadareh et al. Bulletin of the National Research Centre (2022) 46:19 heat stress, drought stress or low water availability across growing seasons. Hence, breeders are always testing large numbers of genotypes in different locations for the development of high-yielding varieties to be adapted to various environments (Khalili and Pour-Aboughadareh 2016). Like other cereal crops, grain productivity in barley has a high economic value. This quantitative trait is governed by the additive main effects of environment (E) and genotype (G), as well as by nonadditive effect of G × E interactions (GEI) (Bocianowski et al. 2021). The GEI effect has a key role in breeding programs because it reveals yield changeability of genotypes across different environments which unexplained by individual E and G effects . Nevertheless, this effect is the biggest challenge for breeders, and it can reduce the heritability of grain yield, decline the correlation between phenotypic and genotypic values (Yan and Fregeau-Reid 2018) and finally complicate the selection of superior varieties in terms of yield performance and stability (Ebdon and Gauch 2002). Therefore, considering the GEI effect in breeding programs increases the efficiency of selection processes and helps to identify high-yielding genotypes with both general and specific (Vaezi et al. 2018). In this regard, the multi-environment trials (METs) have an important role to interpret the GE effect and select superior genotypes at the end of the cultivar development pipeline (Gerrano et al. 2020). In this regard, various statistical models and approaches have been developed and numerous studies have used them to study grain yield stability and adaptability of different crops. The statistical methods can be categorized into two main classes. The univariate methods included several parametric and nonparametric stability statistics such as Huehn's (1979) and Nassar and Huehn's (1987) nonparametric measures (S (i) ), Thennarasu's (1995) nonparametric measures (NP (i) ), Kang's (1988) rank-sum, Plaisted and Peterson's (1959) mean variance component (θ i ), Plaisted's (1960) Wricke's (1962) ecovalence (W i 2 ), a joint regression which in turn included three parameters of the regression coefficient (bi; Finlay and Wilkinson 1963), variance in regression deviation (S di 2 ; Eberhart and Russell 1966) and coefficient of determination (R 2 ; Pinthus 1973), Shukla's (1972) stability variance (σ i 2 ), Francis and Kanenberg's (1978) coefficient of variability (CV). These statistics are common to analyze GEI effect and description of phenotypic adaptability and stability of genotypes. The multivariate methods are grouped into two empirical and biological models. Of these, the empirical models are important including the additive main effect and multiplicative interaction model (AMMI; Gauch 1988) and GGE biplot approach . These methods can efficiently interpret the GEI effects due to several concepts such as 'which-one-where' pattern in data sets, selection of ideal genotypes for use in a wide range of environments and discovering of mega-environments, and determine the discriminating power of test environments (Gauch et al. 2008).
Among both types of univariate and multivariate methods, the AMMI model is widely used to analyze of GEI effects in many plant species. This model combines both integration analysis of variance (ANOVA) and principal components analysis (PCA) in a unique model with fixed effects. Moreover, this model has a graphical tool which can evaluate yield performance and stability simultaneously, as well as the discovery of ideal testing environments (Gauch 1988 Sneller et al. 1997) and the absolute value of the relative contribution of IPCAs to the interaction (Za i ; Zali et al. 2012)-are calculated based on the AMMI model. Vaezi et al. (2017) investigated the grain yield stability of barley genotypes using several stability models and finally reported that the AMMI model is more suitable to identify the superior genotypes in terms of yield performance and its stability. In a study conducted by Ahakpaz et al. (2021), it has been demonstrated that the AMMI model is more useful than other statistical methods for identification and selection of superior barley genotypes for the cold zone in Iran. Despite many advantages of AMMI model in explaining the GEI effect, disability of it in dissecting the structure of the linear mixed-effect model (LMM) is the main disadvantage. To circumvent this issue, Olivoto et al. (2019) introduced a novel model-the weighted average of absolute scores from the singular value decomposition of the matrix of BLUP for GEI effects generated by a LMM (WAASB)--to better characterizing ideal genotypes. Indeed, this model combines the features of best linear unbiased prediction (BLUP) models and AMMI model in a unique index. In this regard, Ahakpaz et al. (2021) and Pour-Aboughadareh et al. (2021) revealed that the WAASB index can be used as an efficient index in selecting the high-yield and stable barley genotypes across multi-environment experiments. Due to the fact that barley ranks second and fourth in the Iran and world, as well as the wide climate changes cause most parts of croplands to encounter high temperature and low water availability, hence the main objectives of the present research were to (i) investigate the grain performance of new advanced barley genotypes, which from part of the national barley breeding program at the Iran, (ii) analysis of the GEI for grain yield using a complex of parametric and nonparametric methods and (iii) identification of high-yielding and stable genotypes across different areas in Iran.

Genetic materials and setup experiments
A set of 17 barley genotypes were derived from hybridization between Iranian local cultivars and International genetic materials received from ICARDA along with two check varieties (cv. Auxin (G 1 ) and WB95-3 (G 19 )) were used in this study (Additional file 1: Table S1). The multi-environment trials (METs) were carried out at five tropical regions of Iran during the 2018-2020 cropping seasons. The test environments included Ahvaz, Zabol, Moghan, Gonbad and Darab. More information about test environments is shown in Fig. 1. Sowing and crop managements in all regions were done based on experts' advice. In the planning stage, fertilizer treatments of 32 kg N ha −1 and 100 kg P 2 O 5 ha −1 were applied to the soil at all research stations. In each experiment, barley genotypes were investigated in a randomized complete block design (RCBD) with three replications. Each experimental plot consisted of six rows, 5 m in length, and spaced 20 cm apart. Plant density was determined as 400 seeds per square meter, and seed sowing was done using an experimental plot planter (Wintersteiger, Austria). At harvest time, an experimental combination (Wintersteiger, Austria) was used to harvest plants and finally grain yields were estimated by converting data to tons per hectare.

Statistical data analysis
The grain yield data were pooled across locations and years and subjected to the AMMI model (Gauch 1993). The AMMI analysis was performed in R software 4.0.3 (R Core Team 2020) using package 'metan' (Olivoto and Lucio 2020). Twenty-three stability statistics were calculated, and further genotypes were ranked based on each statistic. All AMMI-based stability statistics were calculated using 'metan' package in R. Moreover, several parametric and nonparametric stability statistics were calculated using a web-based STABILITYSOFT program . The list of the calculated stability statistics is found in Table 1. Principal component analysis (PCA) was computed to detect interrelationships among calculated stability statistics using 'factoextra' package in R. For grouping investigated barley genotypes and measured traits a hierarchical cluster analysis (HCA) was computed based on the Euclidean distances using 'ggdendro' and 'ggplot2' packages in R.

AMMI analysis
The results of AMMI analysis of variance indicated that grain yield was significantly affected by environments (E), genotypes (G) and their interaction (GEI) (P < 0.001) ( Table 2). The environments explained 46.20% of the total variation in grain yield, while genotyping differences justified only 7.01%. The proportion of GEI in explaining variation of yield performance was 14.94%, revealing the magnitude GEI in MET trials. This result confirmed by the fact that the mean yield across the investigated genotypes varied between 1.63 tons ha −1 (genotype G 18 in environment AHZ1) and 7.24 tons ha −1 (genotype G 5 in environment DAB2), whereas the mean grain yield of genotypes ranged from 3.76 tons ha −1 (genotype G 9 ) to 5.15 tons ha −1 (genotype G 5 ) ( Fig. 2A). Moreover, Fig. 2A shows a high changeability in environment productivity so that it varied between 2.88 tons ha −1 (environment AHZ1 in 2018-19 cropping season) and 5.99 tons ha −1 (DAB in 2019-20 cropping season). The interaction's sum of squares of GEI is further divided into five significant interaction principal components (IPCA1-IPCA5). The IPCA1 (26.80%) and IPCA2 (22.305) together accounted for 49.10% of the total variation due to the GEI ( Table 2).
The distribution of test environments and barley genotypes based on mean grain yield and IPC1 scores is found in Fig. 2B. In this biplot, the mean total grain yield (4.56 tons ha −1 ) is shown using a line perpendicular to the horizontal. Half of the genotypes showed a higher mean grain yield than the total, and among them genotype G 5 and followed by G 19 , G 1 , G 2 and G 4 with 5.15, 5.11, 4.88, 4.81 and 4.81 tons ha −1 had the highest mean grain yield, respectively. Genotypes G 6 (3.77 tons ha −1 ), G 9 (3.76 tons ha −1 ) and G 18 (4.01 tons ha −1 ) showed the lowest mean grain yield compared with other genotypes. The environments showed a distribution pattern like genotypes. There was a big difference between AHZ1 and DAB2 environments. The IPC1 scores genotypes for genotypes G 12 , G 6 , G 13 and G 16 were close to zero; hence, these genotypes can represent have high yield stability (Fig. 2B). The position of test environments and genotypes is shown using AMMI2 biplot (Fig. 2C). Based on this biplot, environments and genotypes close to origin of the biplot show the lowest effects on GEI, whereas environments and genotypes with a large distance from origin of the biplot showed the most influential in creating GEI. Genotypes G 1 , G 2 , G 3 , G 4 , G 14 , G 17 and G 19 have narrow adaptation, whereas G6, G11, G12, G13, G16 and G18 genotypes showed broad stability to all environments. Among test environments, DAB1, GON1 and ZAB2 with the longest vectors compared with other environments showed the high discriminating ability and hence played the highest contribution in creating GEI.  Shukla (1972) 7 Coefficient of variance CVi Francis and Kannenberg (1978) 8-11 Huehn's and Nassar and Huehn's nonparametric statistics S (i) Huehn (1990), Nassar and Huehn (1987) 12-15 Thennarasu's nonparametric statistics NP (i) Thennarasu (1995) 16 Kang's rank-sum KR Kang (1988)
The WAASB statistics was used to better characterize ideal genotypes based on both mean grain yield and stability. Figure 2D shows distribution of barley genotypes based on the mean grain yields and WAASB values. The first quadrant included AHZ1, DAB1, ZAB1, ZAB2 and GON2 environments, as well as genotypes G 3 and G 14 . These genotypes showed lower grain yield compared with the average grain yields. In addition, all environments had lower productivity than the average mean; thus, the genotypes and environment belonging to this section have the largest role in GEI. Genotypes G 4 , G 7 and G 19 along with environments DAB2, GON1 and MOG1 were placed in the second quadrant. These genotypes have a higher grain yield than average yield and on the other hand, environments have a big role in GEI. Besides, the environments placed in this quadrant provide aboveaverage production; hence, they deserve special attention to discriminate the high-yielding genotypes. The G 6 , G 8 , G 9 , G 15 , G 17 and G 18 genotypes with lower grain yield than average yields, along with the AHZ2 environment fell in the third quadrant. Environment AHZ2 with low performance showed the minimum discrimination power in GE interaction. Indeed, the WAASB values for these genotypes and environment were minimum. The fourth part of biplot comprised the rest of genotypes. Hence, G 1 , G 2 , G 5 , G 10 , G 11 , G 12 , G 13 and G 16 genotypes with low WAASB values and high performance compared to average yield were identified as the most stable genotypes (Fig. 2D). The G 13 and G 15 showed yield performance equal to average yields and were placed on the border of the third and fourth quadrants. Furthermore, the status AHZ2 environment was like G 13 and G 15 ; hence, it means that the productivity of this environment was equal to the average performances.

Spearman's rank correlation between stability statistics
A heatmap-based Spearman's rank correlation coefficient was used to investigate relationships between mean grain yield (GY) and the calculated stability indices (Fig. 3).

Clustering of stability parameters and genotypes
The calculated stability statistics were classified into two main clusters. The first cluster (CI) further divided into two sub-clusters (CSI and CSII). The CSI included GY, CV,Pi,NP (2) , NP (3) , NP (4) , S (3) , S (6) and KR stability statistics. In the second cluster (CSII), stability statistics S (1) , S (2) , NP (1) , W 2 , σ 2 , S di 2 , θ i ', R 2 , EV classified together in the same sub-sub-cluster, while ASV, WAASB, Za and SIPC statistics placed a separate the second sub-sub-cluster. The second main cluster lonely included θ i stability  (Fig. 3). In the case of investigated genotypes, a hierarchical cluster analysis was computed based on squared Euclidean distance using Ward's method. The grouping pattern obtained from this analysis showed that all barley genotypes were separated into two main clusters (Fig. 4). The first cluster (CI) further was divided into two sub-clusters (SCI-I and SCI-II). The SCI included genotypes G1, G2 and G10 and had a higher average grain yield than the average mean (4.81 tons ha −1 vs. 4.56 tons ha −1 ) and the lowest ranks for almost all nonparametric stability statistics. Of these, G1 and G10 also had acceptable average ranks (ASR) for parametric stability statistics; hence, they can be selected as the high-yielding and stable genotypes ( Table 4). The SCII comprised genotypes G11, G12, G13, G15 and G16. The average grain yield of this sub-cluster was higher than the average yields (4.65 56 tons ha −1 vs. 4.56 56 tons ha −1 ), and all genotypes had the lowest ASR for parametric stability statistics. The second main cluster (CII) was splatted into three sub-clusters. The first sub-cluster (SCII-I) included genotypes G4, G5 and G19 so that its average Fig. 3 Heatmap showing the relationships between grain yield and stability statistics based on Spearman's rank correlation test. GY grain yield, S (i) Huehn's (1979) and Nassar and Huehn's (1987) nonparametric stability statistics, NP (i) Thennarasu's (1995) nonparametric stability statistics, LR Kang's (1988) rank-sum, W i 2 Wricke's (1962) ecovalance, σ i 2 Shukla's (1972) stability variance, S di 2 deviation from regression, bi slope of regression line, θ i Plaisted and Peterson's (1959) mean variance component, θ i ' Plaisted's (1960) GE variance component, Pi Lin and Binns's (1988) superiority index, R 2 Pinthus's (1973) coefficients of determination, ASV AMMI stability value, SIPC sum of the absolute value of the IPCA scores, EV average of the squared eigenvector values, Za the absolute value of the relative contribution of IPCAs to the interaction, WAASB the weighted average of absolute scores. ns, * and ** indicate nonsignificant, and significant at 0.05 and 0.01 probability levels, respectively grain yields were more than the average yields (5.03 tons ha −1 vs. 4.56 tons ha −1 ). Among these genotypes, only G5 had a minimum ASR value for nonparametric stability statistics. The SCII-II included two genotypes G6 and G9 with lower mean grain yield relative to average yields. Remain genotypes were categorized in the SCII-III. The average grain yield of this cluster was lower than average performance (4.38 tons ha −1 vs. 4.56 tons ha −1 ). All genotypes in this sub-cluster had a higher ASR compared to other genotypes, and therefore, they were recognized as the most unstable genotypes.

Discussion
In the present study, we observed highly significant differences among the investigated barley genotypes that indicate genetic differences between test materials (Table 2). In addition to the significant effect of genotypes, remarkable differences among environments and GEI effects reveal climate change from year to year and location to location clearly affects barley production. Furthermore, the results of AMMI showed that the magnitude of environment (E) and GEI in explaining the total variation of grain yield were higher than genotype (G), showing remarkable differences in the genetic background of barley genotypes across environments (Vaezi et al. 2018). Likewise, Khalili and Pour-Aboughadareh (2016), Vaezi et al. (2017Vaezi et al. ( , 2018, Ghazvini et al. (2018), Kumar et al. (2018), Vaezi et al. (2019), Ahakpaz et al. (2021) and Ghazvini et al. (2021) reported a significant level of GEI for grain yield in barley genotypes. The environments AHZ1 and DAB2 showed the lowest and highest productivity, respectively, and the difference between them was 3.11 tons ha −1 that is more than the grain yield produced in the AHZ1 environment (2.88 tons ha −1 ) ( Fig. 2A). This result could be appreciated for breeders and agronomist for optimization of barley productivity under low-yielding environments in the tropical climate in Iran. Investigation of barley genotypes in the METs reveals that there is a high level of variation in the case of grain yield productivity, yield stability and climate adaptations (Table 2).
Hence, further evaluation of tested genotypes could be used for either developing the high-yielding genotypes or improving grain stability performance and productivity due to suitable breeding strategies. Variability in response to various environments identifies genotypes with different ranks of yield potential and its stability. In addition, to achieve high performance across various environments, stability of genotypes in terms of grain yield is an important component to increase adaptation for cultivation (George and Lundy 2019). Various stability statistics and models have been proposed for evaluating stability of the test genotypes. In the current research, after revealing the existence of significant GEI for grain yield, several stability statistics were served to identify the most stable barley genotypes. In the AMMI analysis, the first two IPCAs accounted for 49.10% of the total variance. Thus, IPCA1 and IPCA2 can represent important details regarding the distribution of genotypes and environments. As shown in Fig. 2B, the distribution of genotypes based on mean grain yield and IPCA1 scores was a useful approach to evaluation of yield performance and stability, simultaneously. The results showed that genotype G12, G13 and G16 with the lowest IPC1 scores and the high mean yield relative to overall average yield can represent to have high yield stability. Besides, the wide dispersion of environments relative to investigated barley genotypes shows high environmental variability compared to genotypes (Ahakpaz et al. 2021). Stability methods commonly are classified into two main concepts: (1) static stability, where a stable genotype maintains its constant yield across environments, in other words, lower stability corresponds to greater sensitivity to diverse environments (Dyke et al. 1995), and (2) agronomic concept or dynamic stability, which points that the GEI is minimum and the mean yield of a stable genotype in each test environment is like the mean response of the tested genotypes. In other words, this concept of stability is dependent on the subset of specific environments (Lin et al. 1986). In this regard, for further dissection of the relationships among stability statistics, these associations were graphically assessed using a heatmap pattern based on Spearman's rank correlation matrix (Fig. 3). In the present study, our results showed a significant correlation between mean grain yield (GY) with Pi, NP (2) , NP (3) , NP (4) , S (3) , S (6) , KR, bi, Za, ASV and WAASB stability statistics. Hence, these statistics are related to the dynamic concept of stability and selection based on each of them is acceptable (Becker and Leon, 1988). This result agrees with findings reported by Scapim et al. (2012), Ahmadi et al. (2015), Khalili and Pour-Aboughadareh (2016), and Vaezi et al. (2018). The results showed that all AMMI-based stability statistics (ASV, SIPC, EV and Za) had a positive and significant correlation with each other and also with most parametric and nonparametric statistics, especially S (1) , S 2) , NP (1) , bi, W 2 , θ i ', S di 2 and R 2 . As an interesting result, WAASB statistic showed a positive and significant correlation with the most stability statistics. Since the WAASB is calculated based on the BLUP approach and uses all IPCAs, this result is not unexpected (Olivoto et al. 2019). However, a close correlation between WASSB and other stability statistics is considerable. Another result obtained from the map pattern is clustering the calculated stability statistics.
The WAASB along with all AMMI-based statistics and several nonparametric and parametric measures were grouped in the same cluster. The GY along with CV, Pi, NP (2) , NP (3) , NP (4) , S (3) , S (6) and KR was grouped together in the same cluster. These results indicate each of stability statistics given a similar ranking pattern for assessing genotypes.
Our results also indicated that various stability statistics gave similar ranking patterns in the selection of stable genotypes (Table 3), suggesting that any of them would be suitable for selecting desirable genotypes. As mentioned heretofore, one of the main objectives of the present study is to determine the efficiency of the WAASB index in identifying the ideal barley genotypes. According to the results presented in Table 3, the WAASB statistic indicates a ranking pattern, for at least one genotype, similar to all stability statistics except S (1) , S (2) , S (6) , NP (1) , NP (4) , Pi and bi. Therefore, the selection of genotypes based on this statistic can result in select high-yielding and stable genotypes. The biplot rendered using the WAASB and mean grain yield was used to better characterize superior genotypes in terms of yield potential and stability. The WAASB' biplot may be helpful over AMMI's biplot in selecting highly productive and broadly adapted genotypes. Indeed, in this way all IPCAs are used, which it turns allows GEI patterns that are not retained in IPCA1 to be considered in the genotypes' ranking (Olivoto et al. 2019). Among test environments, three environments DAB2, GON1 and MOG1 with higher productivity compared to the total mean showed a big role in GEI and therefore can be considered as the special environments for identifying high-yielding genotypes. The genotypes G1, G2, G5, G10, G11, G12, G13, and G16 with the lowest WAASB values and high grain yield were identified as the most stable genotypes (Fig. 2D).
To solve difficulty in selection of superior genotypes based on a single statistic, the average sum of ranks (ASR) of grain yield and stability statistics was calculated. In this method, the low ASR value indicates high level of stability. As per the results, we found that genotypes G1, G5, G10 and G12 had the lowest ASR values in terms of nonparametric statistics, while G10, G11, G12 and G15 had the lowest ASR values in terms of parametric stability statistics (Table 4). With respect to these results, both nonparametric and parametric stability statistics simultaneously identified G10 and G12 genotypes as the most stable genotypes compared to the control (G1 and G19) and other genotypes. These results support the fact that there is a strong correlation between the results of selection patterns achieved from both parametric and nonparametric stability statistics groups. The HCA is one of the powerful multivariate methods to classify genetic materials. In the present study, the HCA was performed based on the ranking matrix of genotypes. Based on obtained dendrogram, the investigated genotypes grouped into two main clusters which further divided into several sub-clusters. The control genotype G1 along with hiyielding and stable genotypes with relatively low ASRs was placed in the first cluster. Indeed, these genotypes may have specific adaptation to some of the environments (Vaezi et al. 2019). The second main group was further divided into three sub-groups, which each of them embraced some genotypes with a range of low to high performance (Fig. 4). However, we did not find any genotype with a low ASR value. Among the classified genotypes in this cluster, only two genotypes G5 and G19 have the potential to use in future barley breeding programs. Considering our results, genotype G10 and G12 showed high performance and acceptable ASR values. With respect to their performances, G10 and G12 produced 4.73 and 4.75 11 tons ha −1 grain yield, respectively, which showed a bit of performance less than both control genotypes. Considering the result of joint regression, the G10 showed bi ≈ 1, whereas G12 had bi = 1.22. We surmise genotypes G10 and G12 have