Genotype by Environment Interaction Analysis of Agronomic Spring Barley Traits in Iceland Using AMMI, Factorial Regression Model and Linear Mixed Model

: Spring barley ( Hordeum vulgare L. ) is the most important cereal in Iceland and its national breeding program aims to select barley genotypes adapted to its environment. A critical step to understand the adaptation of Nordic barley material to a cool maritime climate is to assess the genotype by environment interaction (GxE). In this study, we evaluated the yield and thousand-kernel weight (TKW) of 32 spring barley genotypes in seven Icelandic environments. We applied three methods to analyze GxE: the additive main effects and multiplicative interaction model, a factorial model, and a linear mixed model. For yield, GxE was mainly caused by a better response of six-row genotypes compared to two-row genotypes in high fertility soils. For TKW, GxE showed a pattern along a gradient of daily mean temperatures. This pattern translated into a divergent TKW response between the 2-row and 6-row genotypes, with substantial crossovers along the temperature gradient. This GxE pattern was disentangled using all three methods, illustrating the value of cross-analysis. As yield is the main trait of interest for barley cultivation in Iceland, and few crossovers of genotype performance have been observed between environments, the deﬁnition of one mega-environment was recommended for Icelandic cultivation and breeding. We identiﬁed promising genetic material for both traits and highlighted the superiority of six-row genotypes for yield.


Introduction
Barley (Hordeum vulgare L.) is the fourth most cultivated cereal crop in the world after maize, wheat, and rice [1]. It is well adapted to a broad range of climatic and day-length conditions due to its large diversity [2]. Barley is by far the most important cereal crop in Iceland, despite its location on the northern margin of barley cultivation (63 • -65 • N latitude). Iceland is a subarctic agricultural region characterized by a cool maritime growing season, long photoperiod, risk of frost in spring and autumn, and strong winds. Barley cultivation in Iceland consists of both 2-row (2r) and 6-row (6r) spring varieties, and most of the crop is used as animal feed, while less than 1% is for human consumption [3]. In the last lustrum, the average barley yield in Iceland was 3.142 t ha −1 with an average production of 5560 tons per year [1], indicating that there is still room for increasing grain yield and production when compared with other northern countries such as Norway or Finland [4,5].
The economic viability of barley cultivation in Iceland depends mostly on yield, although grain quality traits are also important [6]. Thousand-kernel weight (TKW) is a highly heritable trait and is correlated with high starch content [7], which is desirable for malting and feed. Bragason [8] proposed that the most important breeding traits for Iceland included earliness and lodging resistance to improve stability. The response of selection for those traits has been positive. For example, in the last 30 years of breeding barley in Iceland, where many barley genotypes were found to be unstable and seemed to perform better in one soil type than the other [36].
The evaluation of GxE has not been carried out in Iceland for spring barley. Until now, the Icelandic barley breeding program has not aimed to select specific adaptations to the various growing regions in Iceland and no formal mega-environment analysis has been conducted. In this study, we focused on the analysis of GxE for the yield and TKW of spring barley. Therefore, the aim of this study was (i) to demonstrate the value of combining several methods for analyzing and interpreting GxE, including AMMI, factorial regression model (FM), and LMM (ii) to characterize the GxE direction and magnitude caused by the environmental factors underlying GxE for both traits of 2r and 6r barley and (iii) to evaluate the existence of mega-environments in the context of Icelandic barley cultivation.

Plant Material
The data we used was from recommendation trials that captured the diversity between Icelandic locations (E) for barley cultivation using 32 representative genotypes (G) to study the interaction of genotypes and environments (GxE) within the same year. In 2017, sixteen spring barley (Hordeum vulgare L.) cultivars (5 Finnish, 4 Swedish, 5 Norwegian, and 2 Icelandic) and 16 breeding lines (12 Icelandic, 3 Swedish, and 1 Norwegian) (Table S1) representing the germplasm adapted to the Nordic climate [37] were grown in seven locations across Iceland ( Figure 1). reports have been published about crossover-type GxE (i.e., different genotype ranks depending on the environment) for soil salinity [34] and acidity [35]. Soil type has been shown to cause variability in the yield and kernel weight of spring barley in Iceland, where many barley genotypes were found to be unstable and seemed to perform better in one soil type than the other [36]. The evaluation of GxE has not been carried out in Iceland for spring barley. Until now, the Icelandic barley breeding program has not aimed to select specific adaptations to the various growing regions in Iceland and no formal mega-environment analysis has been conducted. In this study, we focused on the analysis of GxE for the yield and TKW of spring barley. Therefore, the aim of this study was (i) to demonstrate the value of combining several methods for analyzing and interpreting GxE, including AMMI, factorial regression model (FM), and LMM (ii) to characterize the GxE direction and magnitude caused by the environmental factors underlying GxE for both traits of 2r and 6r barley and iii) to evaluate the existence of mega-environments in the context of Icelandic barley cultivation.

Plant Material
The data we used was from recommendation trials that captured the diversity between Icelandic locations (E) for barley cultivation using 32 representative genotypes (G) to study the interaction of genotypes and environments (GxE) within the same year. In 2017, sixteen spring barley (Hordeum vulgare L.) cultivars (5 Finnish, 4 Swedish, 5 Norwegian, and 2 Icelandic) and 16 breeding lines (12 Icelandic, 3 Swedish, and 1 Norwegian) (Table S1) representing the germplasm adapted to the Nordic climate [37] were grown in seven locations across Iceland ( Figure 1).

Experimental Design, Management, and Growing Conditions
Each trial consisted of a randomized complete block design with three blocks and plots of 10 m 2 (10 seed rows 8 m long and 1.3 m wide). The trials were seeded between 2017-04-26 and 2017-05-07. The trials were situated inside farmers' fields and management was conducted as a standard practice. Fertilizer application rates were decided according to soil types (Table 1), where the lowest N applied was 45 kgN ha −1 and the highest was 120 kgN ha −1 and K and P were not limiting factors. In a cold climate, N is set as a limiting factor to hinder tillering at the end of the growing season, which can be a problem in soils with high fertility. Plot yield was evaluated as tons of seed dry matter per hectare (t ha −1 ) and TKW as grams (g) per 1000 kernels, after drying the samples at 50 • C for three days and cleaning.  Table 1 shows the environmental conditions and management of the trials across the seven environments. The environmental code for each location is the geographical location of each trial denoted by the cardinal and inter-cardinal directions and the level of soil fertility. The locations of the farms were spread around the coast of Iceland, from the 63 • 32 N to the 65 • 46 N. The days from sowing to harvesting were fewest at SE_LF: 136 days, and longest at W_HF: 155 days. The Icelandic meteorology office (www.vedur.is accessed on 30 January 2021) collected all climate data at weather stations that were closest to the trial locations. The data were cropped from the sowing date of each location to its harvest date (all climate data collected are listed in Table S2). The soil types at each location showed great variability in texture. Here, the soil types were categorized into three levels of fertility: low fertility (LF), semi-fertility (SF) and high fertility (HF). The soil texture that was considered least fertile, received the highest amounts of nitrogen (N) while the soil types considered most fertile received the lowest amount of N.

Statistical Methods
To study GxE for yield and TKW, a set of statistical models were applied to the dataset, all derived from a base model as follows: where Y ger is the measured phenotype of genotype g in the environment e and block r, µ is the intercept, α g is the main effect of genotype g, β e is the main effect of the environment e, (αβ) ge is the GxE effect of genotype g in environment e, ω er is the effect of block r in environment e, and E ger is the error associated with genotype g in environment e and block r with E ger ∼ N 0, σ 2 E independent and identically distributed (IID), σ 2 E being the error variance. This base model allows the analysis of the phenotypic variance by splitting it into parts associated with the main effect of the genotype, the main effect of the environment, their interaction (i.e., GxE), and other environmental effects including those accounted for by the experimental design. In this study, we focused on three models to study GxE: the AMMI model, the FM, and the LMM.
The AMMI model, as proposed by Gauch [15], consists in applying the singular value decomposition to the fitted GxE matrix obtained from the base model, which can be decomposed as a sum of multiplicative terms between eigenvectors and singular values. It allows most of the GxE variation to be summarized into a subset of interaction principal components (IPCs), provided that GxE interactions show patterns of (dis)similarities between groups of individuals and groups of environments. The AMMI method consists of the following decomposition of the GxE term of the base model in Equation (1): where λ n is the singular value associated with IPC n and correspondingly λ 2 n is its eigenvalue, N is the number of IPC components, γ gn is the eigenvector of genotype g for IPC n, δ en is the eigenvector of environment e for IPC n, with all eigenvectors scaled as unit vectors. The IPC scores were scaled as λ 0.5 n γ gn and λ 0.5 n δ en , so that their products directly estimate interactions, without the need for another multiplication by λ n . The family of AMMI models was further referred to as AMMIn, where "n" corresponds to the number of IPC components included in the model (e.g., AMMI3 corresponds to an AMMI model with components 1 to 3).
Using the results from the AMMI model, the D statistics of Zhang et al. [38] were calculated for each genotype and each environment as the square root of the sum of the squared eigenvector over all IPCs. This statistic was used as an indicator of the stability of genotypes and environments regarding GxE.
An FM was applied to evaluate the importance of the environmental factor driving the most GxE for each trait: soil fertility for yield and daily mean temperature for TKW. For yield, the main effect of the environment and the GxE effect of the base model in Equation (1) were decomposed by accounting directly for soil fertility as follows: where θ f is the main effect associated with level f of soil fertility with f ∈ {LF, SF, HF}, ε f e is the residual main effect of environment e with the level of fertility f, (αθ) g f is the GxE effect of genotype g for a soil with a level of fertility f, (αε) g f e is the residual GxE effect of genotype g for environment e with a level of fertility f. For TKW, the main effect of the environment and the GxE effect of the base model in Equation (1) were decomposed by directly accounting for daily mean temperature as follows: β e = ηx e + ε e and (αβ) ge = (αη) g x e + (αε) ge (4) where η is the main regression effect of TKW on daily mean temperature, x e is the daily mean temperature measured in environment e (as presented in Table 1), ε e is the residual main effect of environment e, (αη) g is the GxE regression effect of TKW on the daily mean temperature for genotype g, (αε) ge is the residual GxE effect of genotype g for environment e.
The FM for TKW consists of a regression on daily mean temperature while the fertility factor used for yield includes three levels to better account for the GxE variability, which explains the differences in degrees of freedom between the two traits in the analysis of variances (AOV) tables (see below). We applied an LMM approach to evaluate how different environments translate into specific genetic variances and how genetic values correlate between environments. The main genetic effect and the GxE effect were merged into a random variable with the following distribution: where G ge is the random effect of genotype g in environment e, σ 2 Ge is the genetic variance in environment e, ρ ee is the genetic correlation between the environments e and e , and independence is assumed between G ge and E ger for all g, e, and r.
All analyses were performed using R (R Core Team). The significance of the factors of all models was evaluated using an F test. The AMMI model was implemented using the R package agricolae [39]. The GxE sums of squares (SS) associated with "noise" (GxE N SS) were calculated by multiplying the mean squares of the error with the degrees of freedom associated with GxE, as proposed by Gauch and Zobel [12]. The GxE SS associated with "signal" (GxE S SS) was then calculated as the difference between the total GxE SS and GxE N . The relevant variation for identifying mega-environments was calculated as (GxE S SS + G SS)/Tr SS, where G SS corresponds to the SS associated with the main effect of the genotype and Tr SS corresponds to the SS associated with "treatment" (i.e., the joint effect of genotype, environment, and GxE). For FMs, the figure of the expected response according to soil fertility or temperature was obtained by fitting a model including only the main effect of the genotype, the main effect of the tested environmental factor, and the effect of their interaction. Other residual GxE effects were discarded to limit the impact of standardization on such residuals. The LMM was implemented using the R package lme4 [40].

Base Model Analysis of Variance
For the base model, the AOV for yield revealed that the main effect of the genotype (p < 0.001) contributed the most to the variation, 34.2% of the total SS ( Table 2). The main effect of the environment (p < 0.001) contributed 31.2%, and the GxE effect (p < 0.001) 15.2%, of the total SS. Table 2. Analysis of variance for yield using the base model the additive main effect and multiplicative interaction (AMMI) model and the factorial model. With degrees of freedom (Df), sums of squares (SS), mean square (MS), the F-value, p-value and the proportion of SS associated with each factor over the total SS (%SS). For TKW, the AOV for the base model revealed that the main effect of the genotype (p < 0.001) contributed the most to the variation, 42.2% of the total SS, while the main effect of the environment (p < 0.001) contributed 15.9%, and the GxE effect (p < 0.001) 24%, of the total SS (Table 3). Table 3. Analysis of variance for the thousand kernel weight using the base model, the AMMI model and the factorial model. With degrees of freedom (Df), sums of squares (SS), mean square (MS), the F-value, p-value and the proportion of SS associated with each factor over the total SS (%SS). For both traits, the large SS associated with the main genotype effect indicated large differences between genotype means, causing most of the variations for the trait, especially for TKW. The large proportion of SS associated with the main environmental effect for both traits also indicated large differences between environmental means. The magnitude of GxE SS was less than half of the SS associated with the main genotype effect for yield, indicating that genotype-specific responses by environment contributed less to the overall variation than the main effects of the genotype and environment. For TKW, the GxE SS was 57% of the SS associated with genotype suggesting a substantial contribution of GxE to the overall genetic variation.

Additive Main Effects and Multiplicative Interaction for Yield
The AMMI analysis for yield showed that the first IPC accounted for 6.9% (p < 0.001) of the total SS and 45.4% of the GxE SS ( Table 2). The following main IPC components, IPC2 and IPC3, explained further 2.8% (p = 0.002) and 2.2% (p = 0.027) of the total SS, respectively. For yield, the GxE S was 53% of the total GxE SS, indicating that half of the GxE variation could be used for interpretation and prediction purposes. As IPC2-6 captured small proportions of the GxE SS and mostly noise, the AMMI1 model was the most relevant. The relevant variation for identifying mega-environments, defined as the ratio of the sum of genotype and GxE S to the total treatment SS, was 51%. The AMMI1 biplot for yield is shown in Figure 2A and summarizes the variation associated with both the mean response and the IPC1 scores of genotypes and environments. On this biplot, 89.7% of the treatment SS was captured ( Table 2). Other biplots with subsequent IPCs are presented in Figure S1. small proportions of the GxE SS and mostly noise, the AMMI1 model was the most relevant. The relevant variation for identifying mega-environments, defined as the ratio of the sum of genotype and GxES to the total treatment SS, was 51%. The AMMI1 biplot for yield is shown in Figure 2A and summarizes the variation associated with both the mean response and the IPC1 scores of genotypes and environments. On this biplot, 89.7% of the treatment SS was captured ( Table 2). Other biplots with subsequent IPCs are presented in Figure S1. The overall mean yield was 4.4 t ha −1 . The highest mean yields for specific environments were observed at N_HF (5.4 t ha −1 ) and W_HF (5.3 t ha −1 ). These were also the only two environments with high fertility (HF) soils. The lowest mean yield was observed at SE_LF (3.1 t ha −1 ). The highest yielding genotype was a Finnish cultivar: G32 (6r, 5.6 t ha −1 ), and the lowest yielding genotype was G23 (6r, 2.4 t ha −1 ). On average, 6r genotypes The overall mean yield was 4.4 t ha −1 . The highest mean yields for specific environments were observed at N_HF (5.4 t ha −1 ) and W_HF (5.3 t ha −1 ). These were also the only two environments with high fertility (HF) soils. The lowest mean yield was observed at SE_LF (3.1 t ha −1 ). The highest yielding genotype was a Finnish cultivar: G32 (6r, 5.6 t ha −1 ), and the lowest yielding genotype was G23 (6r, 2.4 t ha −1 ). On average, 6r genotypes achieved a higher yield (4.7 t ha −1 ) compared to 2r genotypes (3.7 t ha −1 ). The highest yielding 2r genotype was the Icelandic cultivar G19 (4.2 t ha −1 ) and the lowest yielding 2r genotype was G30 (2.7 t ha −1 ).
The results from AMMI1 model on the biplot in Figure 2A showed that the IPC1 scores for environments, scaled to be in phenotypic units, ranged from −0.9 (S_SF) to 1.7 (W_HF) and the environments N_HF and W_HF were the only ones that had positive IPC1 values, with IPC1 scores of 0.9 and 1.7, respectively. Environments with low fertility, N_LF, SE_LF, and S_LF, had IPC1 values close to zero, −0.1, −0.3, and −0.4, respectively. The semi-fertility environments NE_SF and S_SF had the most negative IPC1 scores of −0.8 and −0.9, respectively. The IPC1 score for genotypes ranged from −1 (6r, G23) to 0.7 (6r, G31). The 2r genotype that showed the lowest IPC1 score was G22 (−0.6) and the 2r genotype with the highest IPC1 score was G8 (−0.1). The 6r and 2r genotypes were mostly separated on the AMMI1 biplot (Figure 2A). On average, 6r genotypes had both a higher mean yield and IPC1 score compared to 2r genotypes. These results indicated that 6r genotypes yielded particularly better than 2r genotypes in the high fertility environments N_HF and W_HF that had both the highest mean yield and IPC1 scores. The gap between 6r and 2r genotypes was more limited in environments with low IPC1 scores such as S_SF or NE_SF.
The ranking of genotypes based on the AMMI1 model (Table S3), declared three winners (G32, G5, and G16) and suggested the existence of three mega-environments, as proposed by Gauch (2013). The genotype G32 showed the highest yield in four environments (SE_LF, W_HF, N_HF, and N_LF). The genotype G16 showed the highest mean yield for environment S_LF. The genotype G5 had the highest yield for environments S_SF and NE_SF. However, these three genotypes were always ranking high in all environments according to AMMI1 (Table S3), and only a few crossovers could be observed between all genotypes along IPC1, as shown in Figure S2A.

Additive Main Effects and Multiplicative Interaction for Thousand-Kernel Weight
The AMMI analysis for TKW showed that IPC1 accounted for 13.7% (p < 0.001) of the total SS and 57% of the GxE SS (Table 3). The following main IPC components, IPC2 and IPC3, explained further 3.7% (p < 0.001) and 2.5% (p = 0.002) of the total SS, respectively. For TKW, the GxE S was 69% of the total GxE SS, indicating that more than two-thirds of the GxE variation was considered reliable for interpretation and prediction purposes. Since IPC1 captured most of the GxE S , AMMI1 was considered as the most relevant model from the AMMI model family. For TKW, the relevant variation for identifying megaenvironments was higher than for yield (72%). The AMMI1 model for TKW is illustrated in a biplot shown in Figure 2B and includes the mean response for TKW and the IPC1 scores, capturing 87.5% of the treatment SS (Table 3). Other biplots with subsequent IPCs are presented in Figure S3.
The overall mean TKW was 34.4g. The highest mean was observed at S_LF (36.3 g) and the lowest at N_HF (32.5 g). Environments that were classified with low fertility (S_LF, SE_LF, and N_LF) had the highest TKW mean estimates. The genotype with the highest mean TKW was a Swedish breeding line, G25 (2r, 40.4 g). The genotype with the lowest mean TKW was G23 (6r, 28.7 g). On average, 6r genotypes had a lower TKW (33.5 g) compared to 2r genotypes (36.3 g). The 6r genotype with the highest TKW observed was the Icelandic breeding line G31 (35.9 g) and the lowest TKW mean observed for a 2r genotype was G28 (32.6g).
The AMMI1 biplot showed that all three environments located in the south of Iceland (S_SF, S_LF, and SE_LF) were the only environments with positive IPC1 scores (2.6, 2.1, and 0.5) for TKW ( Figure 2B). When looking at the IPC1 score for genotypes, the results showed that the 6r and 2r genotypes were mostly separated groups ( Figure 2B). The Icelandic breeding line G28 was the only 2r genotype that showed a negative IPC1 value (−0.2). The 6r genotype with the highest IPC1 score was the Icelandic breeding line G26 that scored 1.4. As for yield, the 6r and 2r genotypes were mostly separated on the AMMI1 biplot ( Figure 2B). On average, 2r genotypes had both a higher mean TKW and IPC1 score compared to 6r genotypes. Our results showed that 2r genotypes had a positive GxE with environments in the south of Iceland while 6r genotypes showed a negative GxE with these same environments.
The ranking of genotypes based on the AMMI1 model (Table S9) declared two winners (G25 and G14) and the existence of two mega-environments. The genotype G25 showed the highest TKW for five environments (N_HF, N_LF, NE_SF, SE_LF, and W_HF), while G14 showed the highest mean TKW for environment S_LF and S_SF. These two genotypes were always ranking high (top three) in all environments according to AMMI1. However, we observed several crossovers along IPC1 between the 2r and 6r genotypes ( Figure S4A).

Stability Analysis
Results from the stability analysis of genotypes and environments for both traits based on the D statistics are listed in Table S15. The range of genotype D values for yield ranged from 0.3 to 1.4, with low and high values indicating stable and unstable genotypes, respectively. The genotype G19 had the lowest D value and G23 had the highest D value. Environment D values ranged from 1.4 to 1.9. The highest value was observed for W_HF and the lowest for S_LF. The three environments with the lowest D values were located in the south of Iceland. For TKW, the range of genotype D values ranged from 0.4 to 2.6. The genotype G27 had the lowest D value and G22 had the highest D value. Environment D values ranged between 2.6 and 3.2. The environment with the lowest D value was NE_SF and the environment with the highest D value was S_SF.

Environmental Factors Driving GxE
For yield, soil fertility was included in a FM with three levels: LF (low fertility), SF (semi-fertility), and HF (high fertility). It explained most of the variability associated with the main effect of the environment: 19.3% of the total SS and 61.9% of the environment SS (Table 2). Soil fertility also accounted for a large part of the variability associated with the GxE effect: 7.5% of the total SS and 49.3% of the GxE SS. On average, environments characterized with high fertility soils yielded more compared to those characterized as less fertile. The response of each genotype according to the level of fertility is presented in Figure 3A. We observed a larger variability in high fertility soils compared to low and semi-fertility. This higher variability results mostly from a differential response between 2r and 6r genotypes; 6r genotypes responding better to high fertility soils compared to 2r genotypes (i.e., a larger increase in yield response is observed for 6r genotypes between SF and HF soils compared to the increase observed for 2r genotypes).  Figure 4 shows the genetic correlations between environments and genetic variance estimates from the LMM analysis. For yield, little or no pattern was observed based on genetic correlation between environments, with estimates ranging from 0.73 to 1. The highest genetic variances estimated were found in HF environments (1.37 for N_HF and 1.90 for W_HF), while the genetic variance estimates for SF and LF environments were all below 1.0. For TKW, a clear pattern was observed according to the geography based on genetic correlation between environments. Two groups of environments can be distinguished on the heatmap: (i) environments located in the south and southeast of Iceland (S_LF, S_SF, and SE_LF) and (ii) all other environments (NE_SF, N_HF, N_LF, and W_HF). Within each of these clusters, high genetic correlations were estimated between environments, while low correlations were estimated between the environments of different clusters (down to 0.23 between S_LF and W_HF). Genetic variance estimates in the south and southeast cluster showed higher values (from 6.34 to 13.27) compared to the cluster of west, north, and northeast (from 4.47 to 7.75). For TKW, the daily mean temperature was included in a FM. While the daily mean temperature explained a large part of the variability associated with the GxE effect (7.5% of the total SS and 31.2% of the GxE SS), this factor accounted for almost none of the variability associated with the main environmental effect: 0.5% of the total SS and 3.1% of the GxE SS (Table 3). On average, environments with higher daily mean temperatures showed lower TKW. However, the 2r and 6r groups responded in a divergence pattern along the temperature gradient. The higher the mean temperature, the lower the TKW for 6r genotypes, while 2r genotypes responded with higher TKW at higher mean temperatures. Figure 3B shows the response of each genotype according to the daily mean temperature gradient. The genotype with the highest mean TKW (G25) was stable across the temperature gradient. However, the genotype with the second highest mean TKW (G14) had a lower TKW in environments with lower daily mean temperatures, but exhibited the highest TKW in the environments with the highest daily mean temperatures ( Figure 3B). Figure 4 shows the genetic correlations between environments and genetic variance estimates from the LMM analysis. For yield, little or no pattern was observed based on genetic correlation between environments, with estimates ranging from 0.73 to 1. The highest genetic variances estimated were found in HF environments (1.37 for N_HF and 1.90 for W_HF), while the genetic variance estimates for SF and LF environments were all below 1.0. For TKW, a clear pattern was observed according to the geography based on genetic correlation between environments. Two groups of environments can be distinguished on the heatmap: (i) environments located in the south and southeast of Iceland (S_LF, S_SF, and SE_LF) and (ii) all other environments (NE_SF, N_HF, N_LF, and W_HF). Within each of these clusters, high genetic correlations were estimated between environments, while low correlations were estimated between the environments of different clusters (down to 0.23 between S_LF and W_HF). Genetic variance estimates in the south and southeast cluster showed higher values (from 6.34 to 13.27) compared to the cluster of west, north, and northeast (from 4.47 to 7.75).

Analyzing GxE Based on a Cross-Analysis
The strength of AMMI analysis is the possibility to simultaneously visualize genotypes and environments on biplots that facilitate the interpretation of main GxE features. On biplots, genotypes located close to environments with extreme IPC scores are likely to show a positive interaction in those environments (i.e., higher performance compared to their mean performance), and conversely for genotypes located far from those environments.
However, AMMI does not directly address the issue of identifying the environmental factors driving GxE but relies on correlation coefficients between IPC scores and environmental covariables to identify those factors [13]. To test this, we fitted a model with the environmental covariable that takes the most account of the GxE variability, which were soil fertility for yield and daily mean temperature for TKW. The response was almost identical to that observed along IPC1 using AMMI (see Figures 3B and S4A), which supports daily mean temperature as being the main driver of GxE for TKW. A similar observation was made for yield and soil fertility (Figures 3A and S3A).
Our results for yield showed that defining mega-environments in terms of winners as proposed by Gauch [13] could be misleading. The winner method did not categorize the seven environments into mega-environments according to clear patterns of geography, soil fertility, or climatic covariables ( Table 1). The differences between winners were most often negligible and did not reflect large differences between remaining individuals, as illustrated in Figure S2A. The LMM disentangled the part of the GxE variability that involves differences in genotypes ranks (crossovers) from the GxE variability that does not (convergence/divergence). While crossovers translate into low to moderate genetic correlations between environments, convergence/divergence translates into environmentspecific genetic variances along with high genetic correlations. The GxE variability observed for yield mostly resulted from different genetic variances between environments. Highly fertile environments tended to show higher genetic variances compared to low and semi-fertile environments. Genetic correlations for yield between environments were high ( Figure 4A) and suggested the good conservation of genotype ranks between envi-

Analyzing GxE Based on a Cross-Analysis
The strength of AMMI analysis is the possibility to simultaneously visualize genotypes and environments on biplots that facilitate the interpretation of main GxE features. On biplots, genotypes located close to environments with extreme IPC scores are likely to show a positive interaction in those environments (i.e., higher performance compared to their mean performance), and conversely for genotypes located far from those environments.
However, AMMI does not directly address the issue of identifying the environmental factors driving GxE but relies on correlation coefficients between IPC scores and environmental covariables to identify those factors [13]. To test this, we fitted a model with the environmental covariable that takes the most account of the GxE variability, which were soil fertility for yield and daily mean temperature for TKW. The response was almost identical to that observed along IPC1 using AMMI (see Figure 4B and Figure S4A), which supports daily mean temperature as being the main driver of GxE for TKW. A similar observation was made for yield and soil fertility ( Figure 4A and Figure S3A).
Our results for yield showed that defining mega-environments in terms of winners as proposed by Gauch [13] could be misleading. The winner method did not categorize the seven environments into mega-environments according to clear patterns of geography, soil fertility, or climatic covariables ( Table 1). The differences between winners were most often negligible and did not reflect large differences between remaining individuals, as illustrated in Figure S2A. The LMM disentangled the part of the GxE variability that involves differences in genotypes ranks (crossovers) from the GxE variability that does not (convergence/divergence). While crossovers translate into low to moderate genetic correlations between environments, convergence/divergence translates into environment-specific genetic variances along with high genetic correlations. The GxE variability observed for yield mostly resulted from different genetic variances between environments. Highly fertile environments tended to show higher genetic variances compared to low and semi-fertile environments. Genetic correlations for yield between environments were high ( Figure 4A) and suggested the good conservation of genotype ranks between environments. On the other hand, the GxE variability observed for TKW resulted both from higher genetic variances in the south of Iceland compared to other regions, and a clear clustering between these same regions according to genetic correlations ( Figure 4B). It suggested a ranking of genotypes substantially different between these clusters of regions for that trait.

Exploiting GxE Analysis for Breeding in a Marginal Climate
The yields observed in our experiment are comparable to the yields observed in similar trials across the globe [33] and the national averages in the Nordic region [1]. TKW, although defined as a yield component, showed less correlation with yield in the data presented here than in previous trials in Iceland [9].
Based on the FM (Figure 3) and the AMMI1 model ( Figure 2), a clear distinction between 2r and 6r genotypes can be observed in terms of yield and TKW, according to their mean performance and their GxE response. Contrasting yield responses of 6r and 2r between environments have already been described in Nordic trials in 1988-1989 by Nurminiemi et al. [30], where 2r yielded on average more than 6r genotypes in Iceland (e.g., the 2r Mari released in 1960). The results suggest that breeding strategies for Northern climates have since been more focused on the 6r rather than the 2r, or at least been more successful. The low TKW of Icelandic 2r breeding lines observed in our data (e.g., G19 and G28) may have resulted from larger breeding efforts on yield compared to grain quality.
Soil fertility showed to be the main driver of GxE for yield (Table 2) according to our FM, showing its usefulness in our investigation, as shown before in GxE analyses [23]. Higher yields of 6r genotypes compared to 2r has never been connected to soil characteristics, to the best of our knowledge. Peaty soil textures with high organic content in Iceland, like in environments W_HF and N_HF, have been reported to contain more available nitrogen [41] and a higher water capacity [42] compared to other soils. Barley has been shown to produce more ears under conditions of plentiful water [43]. Under such conditions, 6r shows a greater yield compared to 2r genotypes but led to a lower TKW in fertile soils as previously reported by Hilmarsson et al. [36]. The divergence pattern observed in Figure 3B, where the TKW of 6r genotypes decreases with increasing temperature unlike for 2r, can be explained by the higher number of ears at higher temperatures having a greater impact on TKW for 6r compared to 2r genotypes. This GxE pattern may result from a complex physiological divergence between 6r and 2r genotypes, and different germplasm [44].
It is important for breeding and cultivar recommendations to select genotypes that are stable across environments. The stability of genotypes is often assessed using AMMI biplots ( Figure 2) but when the first three IPCs in the AMMI analysis are significant (Tables 2 and 3), the D statistics is a better estimator of genotype stability [38,45]. The results presented here do not support the reports on the decline of yield stability for barley in the Nordic region [30]. According to AMMI1 model and D statistics, the most stable genotype for yield was G19, the Icelandic cultivar Iskria, which was also the highest yielding 2r genotype, and was stable for these two traits across soil types [36]. G24, an elite Icelandic breeding line was unstable based on D statistics but not AMMI1. This suggests that relying on a single IPC is not precise enough to estimate stability, as has been previously shown [24].
In plant breeding, trial environments should represent the cultivation regions so that GxE can be accounted for when selecting the highly performing genotypes. A breeding program does not necessarily need to include a large number of environments but rather include environments in which large variance can be observed [45,46]. The LMM showed that high fertility environments had a greater genetic variance for yield, graphically supported by the FM (Figure 3) and AMMI ( Figure S2A). Therefore, the fertile peaty soil texture at W_HF is ideal to compare genotypes for yield. The warm mean temperatures in the south of Iceland were associated with large genetic variances for TKW, making S_SF an ideal testing location for this trait.
The absence of major crossovers for yield did not support the definition of megaenvironments for cultivar recommendation or breeding. These results illustrate the issue of delimiting mega-environments using only the winners of a given AMMI model (three different winners for yield using AMMI1), as initially proposed by Gauch and Zobel [12]. We propose instead the use of all genotypes rather than just winners to determine megaenvironments. From a statistical point of view, we recommend the joint interpretation of AMMI results and the genetic correlation estimates obtained using the LMM. Moreover, it has been shown that LMM can outperform the AMMI model in terms of prediction accuracy (e.g., with oat yield [24]) suggesting adequate modeling of GxE using LMM. The data suggest considering Iceland as a single mega-environment for barley cultivation and breeding in Iceland. If TKW would become the main trait for breeding, there would be an opportunity to define mega-environments according to geography.
Future research should emphasize more on better characterizing GxE for these two traits by choosing relevant trial locations featuring a spectrum of soils and temperatures. Additionally, molecular markers can be used to improve the accuracy of GxE estimates and identify regions of the genome that are responsible for GxE by showing a differential effect on traits depending on the environments [11,47]. Such data are a great opportunity to understand further GxE in barley and the physiological pathways responsible for the interaction effect.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.