Analysis of genotype-by-environment interaction in a multisite progeny test with Scots pine for supporting selection decisions

In multisite forest genetic experiments, the presence of genotype × environment interaction (GEI) is common. GEI may negatively affect the estimates of genetic variance and hamper selection decisions in tree breeding programs. Several measures exist to evaluate the stability of tested genotypes’ performance across environments with a choice of the method likely affecting breeders’ decisions. In this study, we evaluated variation in diameter and height growth performance in the progeny test established at 4 sites with 80 open-pollinated half-sib families of Scots pine. We found significant variation among examined progeny at age 10, reaching up to 31% for diameter and 20% for height depending on site, and significant GEI in both traits. We estimated contribution of each family to GEI using various methods and tools of GEI analysis—AMMI, GGE biplots, heterogeneity of regression coefficients (bi’s), the deviation mean squares from regression (s2di) and Kang’s yield-stability index (YSi). Despite the presence of the cross-over interaction, family ranks did not vary much among sites. The selections based on the phenotype, YSi and restricted bi corresponded well to each other leading to the expected response to selection up to 7.8% on diameter and 4.4% on height, whereas those based on the AMMI stability variance were different and lead to a slight loss in both traits. We discuss our results in the context of the usefulness of those measures of genotype stability for tree breeding programs and propose the procedure to follow for making selection decisions in forest experiments with GEI.


Introduction
Progress in tree breeding programs is made by evaluation of genetic worth of selections through progeny testing (White et al. 2007). Progeny tests could also indicate possible deployment zones for families included in the tests, although they are not specifically designed for this purpose. However, estimation of genetic variance and genetic effects is commonly hampered in multisite experiments by the presence of genotype × environment interaction (GEI).
These interactions arise from a differential response of genotypes to environmental factors and biotic and abiotic stresses (Kang 2002;Li et al. 2017). Presence of GEI complicates breeding and deployment decisions, especially when important cross-overs occur in genotype rankings among sites. The GEI could be mitigated through homogenization of environmental conditions by grouping similar environments. However, if the objective is to evaluate the performance of genotypes over many environments, it requires appropriate ways of dealing with genotype × site interaction.
When phenotypic traits of genotypes are assessed across multiple environments it is imperative to investigate phenotypic stability. There are at least two concepts of stability. In a static or biological concept, a stable genotype's performance is unchanged regardless any variation in environmental conditions. In contrast, in a dynamic or agricultural concept of stability, the stable genotype has performance that closely follows the average response to environmental variability (Becker 1981;Becker and Léon 1988). It is a Communicated by Oliver Gailing.
1 3 dynamic stability that is typically sought in tree breeding programs.
Several methods exist that help breeders to analyse GEI, each with its advantages and limitations (Becker and Léon 1988;Crossa 1990;Reckling et al. 2021). The regression approach decomposes the GEI sum of squares from ANOVA into the component associated with heterogeneity of regression coefficients (b i 's) among genotypes and the remainder (Finlay and Wilkinson 1963). Eberhart and Russell (1966) proposed the deviation mean squares from regression (s 2 di ) as genotype's stability measure, and both the b i and s 2 di are used to describe genotype's contribution to GEI. However, both phenotypic stability and yield are important in assessment of genotype's performance. Therefore, Kang proposed the YS i statistic (Yield-Stability Index), integrating yield and stability into a single measure (Kang 1993). The YS i uses the Shukla's "stability variance" as a stability measure, which is derived as a linear combination of squares of residuals (Shukla 1972).
The multivariate methods of GEI analysis that use singular value decomposition (SVD) are very attractive, because they allow to recognize the patterns that exist in data with the help of bivariate graphs (biplots). The AMMI method (additive main effects and multiplicative interaction) integrates the analysis of variance and principal component analysis (PCA) (Crossa 1990). It applies ANOVA to the additive effects of genotypes and environments in the model, and the PCA to the interaction term (Gauch 1988;Gollob 1968). The other, somewhat similar method, is the GGE biplot analysis (Yan and Tinker 2006;Yan et al. 2007). The term "GGE" implies that genetic (G) and interaction (GE) effects should be considered simultaneously for the appropriate evaluation of tested genotypes and environments in GEI. Thus, although both the AMMI and the GGE biplot use the SVD, the GGE biplot method takes into consideration not only the interaction term, but also genetic effects, which allows for more effective visualization of those factors (Yan and Tinker 2006). The GGE biplots have been popular in analysis of GEI in agricultural research on cultivars tested in multiple environments (Yan et al. 2007), but has been rarely applied to forest trees (Correia et al. 2010;Jastrzębowski et al. 2018;Ukalski and Klisz 2016). However, it is a very useful tool for visualizing relationships among test environments and ranking genotypes based on their performance in relation to environments. Preferably, several stability measures need to be integrated with the assessment of genotype's performance for effective selection in breeding programs.
With various methods to analyse GEI and investigate patterns in the data, it is difficult to choose a single method that would ideally serve the breeder to rank genotypes and assess the stability of their performance across multiple environments. The primary objective of this study was to evaluate variation in diameter and height growth performance in the progeny test established at 4 sites with 80 open-pollinated half-sib families of Scots pine (Pinus sylvestris L.), an important European coniferous species. We aimed to estimate contribution of each family to GEI using various methods and tools of GEI analysis, and interpret the usefulness of those measures for tree breeding programs.

Planting material and test sites
Progeny test comprises of the open-pollinated progeny of 80 Scots pine trees. We refer to them as to families or genotypes for clarity. This nomenclature is common in forestry applications of GEI analysis, although those half-sib families represent the mixtures of genotypes, in contrast to pure lines used in agriculture. The four experimental sites were established in spring 2010 in Northwestern Poland (Table 1). At each site, the experiment was planted in a single-tree plots design. One-year-old containerized seedlings were planted in 1.5 × 1.5 m spacing (4444 trees ha −1 ). According to the original design, a total of 100 trees per family was supposed to be planted at each site. However, due to mortality of seedlings in the nursery, the tree number per family varied at the time of planting between 63 and 110 in Podanin, 69 and 112 in Tuczno, 81 and 103 in Płytnica and 88 and 105 in Sarbia. After the first year, the dead spots at planting sites were re-planted using the extra seedlings, preferably of the same Table 1 Information on planting sites of the analyzed Scots pine progeny tests. Survival is given without the re-planted trees a Currently in the Kalisz Pomorski forest district, b According to the Polish forest site typology, c Long term  MAT mean annual temperature, MAP mean annual precipitation from the WorldClim database (30-arc sec resolution, https:// www. world clim. org/) family. Those trees were excluded from further analyses. No thinnings were made to the age of the reported measurements, and thus, mortality was due to the natural causes and survival varied between 75 and 83% among sites (Table 1).

Measurements
In 2019, at tree age 10 years, the diameters (DBH at 1.3 m) of all trees were measured and heights on about the half of the trees' number. For the remaining trees, the heights were imputed with the Prodan 3-parameter model parameterized for family within site in the 'lmfor' package in R (R Core Team 2019).

Analysis
The first step in statistical analysis was to perform the ANOVA according to the model: where y ijk is the individual tree phenotype, µ is the overall mean, S i is the effect of site i, R j(i) is the effect of replication j within the site i, F k is the effect of family k, FS ik is the family × site interaction effect, and e ijk is the error. The same model was fitted with all random effects for estimation of variance components; then family heritability was calculated, according to the formula given in Jayaraman (1999). Since the number of trees per family was not uniform across sites and families, we used the harmonic mean of replications N = 76.65 for heritability calculation. Following the ANOVA, the AMMI analysis was performed. This procedure fits the additive main effects (for genotypes and environments) by the usual analysis of variance and applies the principal component analysis (PCA) to the interaction part of the model (Gauch 1988;Gollob 1968). The AMMI stability value (ASV) was calculated based on the first two AMMI PCs according to the formula given in Purchase et al. (2000). The lower ASV value indicates a greater stability of the genotype.
Kang's Yield-Stability index YS i is based on the Shukla's stability variance (Shukla 1972). The procedure first adjusts the genotype ranks based on the amount of variation among genotypes and combines it with the rank founded on the significance level for stability variance (Kang 1993). In other words, it penalizes the genotype's rank for excessive variation of performance across environments. The high YS i values are favorable in terms of selection decision.
The b i measure of GEI analysis is the regression coefficient from the regression of family performance on environmental index (Finlay and Wilkinson 1963). In the absence of independent environmental index describing site variation, the deviation of average for all families in a particular environment from the overall mean is used as an environmental index (Eberhart and Russell 1966). Genotypes with the b i values close to 1.0 show the average responsiveness to environment and are examples of dynamic stability (Becker 1981;Becker and Léon 1988). Genotypes with the b i values higher than 1.0 show a greater that average response to the improved environmental index, thus they are unstable in a good sense, but they may underperform in poor environments. On the other hand, genotypes with the b i values lower than 1.0 are not responsive to the improvement in environment; they represent a biological concept of stability. The deviation mean square from regression-s 2 di is related to the unexplainable remaining variation in genotype's performance (Eberhart and Russell 1966). The stable genotype should have those deviations close to zero.
The GGE biplot analysis, similarly to AMMI, uses the PC analysis. However, in contrast to AMMI, in this method, the contribution of each genotype to genetic (G) and interaction (GE) effects is considered simultaneously (Yan and Tinker 2006). With different methods of data centering, data scaling and singular value partitioning, the patterns of interaction could be neatly visualized in this graphic method to assess the contribution of environments and genotypes to the GEI (Yan et al. 2007).
Correlation coefficient between different measures of GEI was calculated. Further, the GEI measures were compared through the simulation of selection for 10 top families according to each method. The response to selection was calculated as selection differential × family heritability. Selection based on the b i was actually constrained to topperforming families having the b i values not significantly different from 1, and the s 2 di values not significantly different from zero. The stability variance was not used in that comparison, because it is incorporated into the YS i .
All the stability analyses were performed in the R environment (R Core Team 2019). The ANOVA was fit with the 'lmer' procedure in the 'lme4' package. The AMMI analysis and YS i calculations were performed in the 'agricolae' package. The b i and GGE analysis was done in the 'metan' package, and visualization of GGE biplots was done in the 'GGEBiplots' package.

Results
Analysis of variance revealed the existence of significant GEI for both diameter and height (Table 2). Sites were significantly different, with the same ranking of sites with respect to both traits. Greater values were observed at the Podanin site compared to the three other sites (Figs. 1, 2) with many cross-overs for families among sites (Fig. S1). Variation of tree diameters among families amounted to between 27% in Płytnica and Podanin and 31% in Tuczno  (Fig. 1), and for tree height it amounted to between 13% in Sarbia and 20% in Płytnica and Tuczno (Fig. 2).
Further analysis of GEI with AMMI revealed three significant principal components for both traits ( Table 2). The PC1 and PC2 jointly explained 72.5% and 71% of GEI for diameter and height, respectively. The AMMI stability values (ASV) based on the first two PCs are shown in Fig. 3, with lower ASV values indicating greater stability. The AMMI biplots are shown in Fig. S2. Although the diameters and heights were highly correlated at the family level (Pearson's r = 0.88), the correlation between the ASV values on both traits was only slight (r = 0.37, Fig. 3). All measures of stability from the GEI are shown in Table S1.
The analysis of GGE biplots did suggest the existence of separate mega-environments for diameter (Podanin vs. other sites) and height (Tuczno vs. other sites). The Tuczno site was the most representative for the average-environment for diameter, and the Sarbia site for height, although these sites were also the least discriminative among families (Fig. S3). The closest family ranks were observed between the sites in Podanin and Płytnica for height with the greatest differences between those two sites and the Tuczno site (Fig. S3). For diameter, the most disparate family rankings were between sites in Podanin and Sarbia (Fig. S3b). However, the family ranks were positively correlated among all sites (Spearman's rho between 0.52 and 0.7 for diameter, and between 0.55 and 0.71 for height, all p ≤ 0.0001), thus data were further treated as belonging to a single mega-environment ( Fig. 4 and Fig.  S3, Table S2). In the figure of mean performance vs. stability, families having the long distances between their point coordinates and the average-environment axis (the close-tohorizontal line running through the origin) contribute much to the GEI, and thus show poorer stability than families closer to that axis (Fig. 4). Among the families contributing the most to GEI were 1883, 1888, 8438, 6461, 3936, 7488 and 7493 for diameter, and for height these were 7488, 1881, 8436 and 3935 (Fig. 4, see also Figs. 1, 2).
Analysis of regression coefficient of family performance on environmental index showed that the heterogeneity of regression coefficients was not significant for both traits (Table 2). However, the b i values varied between 0.56 and 1.44 for diameter and between 0.80 and 1.19 for height (Fig. 5). For both traits five families had the b i values significantly greater than 1.0, indicating their greater than average responsiveness to improvement of environmental conditions, with two of those families in common (3718, 6461) for diameter and height (Fig. 5). Five and seven families had the b i values lower than 1.0 (lower than average responsiveness) for diameter and height, respectively, with four of them in common for both traits (Fig. 5). Twelve and 15 families had the s 2 di values significantly greater than zero (P ≤ 0.05) for diameter and height, respectively (Fig. 5). Seven of those families were also identified as contributing the most to GEI for diameter (as above), but only one (1881) for height (Figs. 4, 5, Table S2).
The Kang's Yield-Stability Index (YS i ) was highly correlated with phenotypic values. A weak positive, but significant correlation was found between the s 2 di and ASV and phenotype for diameter (Table 3). The s 2 di and ASV were also positively correlated with each other and with stability variance for both traits, but the strength of those correlations was greater for diameter than height. Correlations between the other measures of GEI were mostly weak and nonsignificant for both traits (Table 3).
The estimates of family heritability used for calculating selection differential were 0.894 (± 0.088 s.e.) for diameter and 0.891 (± 0.088 s.e.) for height. Family selections based on the phenotype and YS i corresponded very well and resulted in 5 mm (7.6%) response to selection on diameter and 21 cm (4.3%) on height (Table 4). There was a smaller agreement with family selection based on the b i (6.3% response for diameter and 3.3% for height), and much

Discussion
We have found significant variation among the open-pollinated families of Scots pine in the young progeny test with respect to height and diameter growth, with significant GEI present for both traits when tested at four sites. The majority of variation was associated with site effect. Because the sites differed little in climatic conditions, the majority of that variation was likely caused by variable soil fertility. This was reflected in the fact that site rankings in both traits closely followed the gradient of site type fertility according to the site typology used in forestry in Poland. Relative differences among sites were only slightly greater for height than for diameter. The examined stands were in the phase of initiation of strong inter-tree competition, thus the observed pattern indicates that height growth may gain priority over diameter growth in a closing stand, especially at the most fertile site in Podanin. Significant GEI was reflected in cross-overs of family ranks among sites. This finding ascertains that selection decisions should take this interaction into account. However, no single method of GEI analysis serves all the objectives that the breeder may have when making selection decisions.
The GGE biplots are very useful in disentangling the patterns of contribution of each genotype and environment into both genotype and genotype × environment interaction effects. The method proved its usefulness for the forestry applications (Correia et al. 2010;Klisz et al. 2017;Ukalski and Klisz 2016), which is also confirmed by our study. However, it requires balanced data which may limit its application. Yet, this constrain may be overcome with the use of BLUP values for examined genotypes in connection with GGE biplots (Ling et al. 2021). Conveniently, in our study we identified with the help of GGE biplots that even with the presence of the cross-over interaction, family ranks were not strikingly different among sites, thus there was no need for separate sets of selections for our sites. Despite its attractiveness, an important drawback of the GGE biplot method is the lack of uncertainty measures applied to it. The mean performance and stability/instability for genotypes could be easily identified on the biplot graph, however, any interaction patterns identified by this method should be tested with formal statistical procedures (Yan and Hunt 2002).
The AMMI procedure is also a multivariate method using the PC analysis. The advantage of AMMI over the GGE biplots is the availability of the ASV, which is a quantitative measure of genotype stability (Purchase et al. 2000;Sabaghnia et al. 2008). The ASV shows the distance from the origin on the biplot of the first two principal components weighted by the proportional difference between these two PCs, because the PC1 contributes the most to the G × E sum of squares (Purchase et al. 2000). The ASV values in our study were well correlated with the other measures of stability based on squared residuals (stability variance and s 2 di ) for diameter, but less so for height. Thus, it confirms its suitability as a metric of genotype stability, but it also requires further testing in broader set of environments for examination of factors affecting this measure. Despite of the availability of ASV in the AMMI method, the GGE biplots are easier to interpret than the AMMI biplots when examining genotype performance at particular sites and assessing its contribution to GEI (Yan and Tinker 2006;Yan et al. 2007).
The heterogeneity of regression coefficients among examined families was not significant, which limits the usefulness of regression approach for the interpretation of GEI in our dataset (Becker and Léon 1988). Relatively few families were either more or less responsive to the change in environmental conditions than the average. However, the regression coefficient should be interpreted together with deviation from regression, where the b i is regarded a response parameter and s 2 di a stability parameter (Becker and Léon 1988).
The method enabled identification of particularly interesting families in our study (see Figs. 1,2,5). Two families-7485 and 1899, had especially favorable growth on both diameter and height and average stability. Families 8444 and 8438 also had good growth, but with large s 2 di and with lower and greater than average stability, respectively. It indicates that the former one would underutilize the better site conditions, and the latter would underperform at poor sites. On the opposite end was the family 7484, which was poorperforming for both traits at all sites with below-average responsiveness to the environment. Our results indicate that, even when variation in regression coefficients for GEI is not large, the detailed examination of the results of regression approach is useful for disentangling the interaction pattern and making selection decisions.
The Kang's Yield-Stability Index (YS i ) was highly correlated with phenotypic values in our study. Perhaps the small level of variation in relative performance of families across sites has caused that, because despite incorporating the stability variance in itself (Shukla 1972), the YS i seems to put more weight on family performance than on its variation. Thus, the YS i was not informative in terms of the examination of the pattern of GEI present in the analyzed data. However, this method was mainly used for agricultural applications (e.g., Mohammadi et al. 2010). In situations similar to our study, the Kang's index should be used as an auxiliary method to other measures of GEI for supporting selection decisions.
Different measures of genotype stability used in our study were only weakly correlated, except for the correlation between ASV, s 2 di , and stability variance. The correlation between stability variance and s 2 di is not surprising as both parameters are based on squares of residuals (Becker and Léon 1988). However, stability of yield is not the only parameter important for selection, but rather an integrated analysis of stability and yield is needed for making informed selection decisions. Therefore, selection based on the ASV which focused only on stability, gave the most disparate results, compared to the selection based on YS i and b i , which gave results similar to selection based on phenotype only. However, if selection was confined to only stable genotypes (b i not different from 1.0 and s 2 di not different from 0) the two families 8444 and 8438 would be missed, although they may be especially suitable for poor and fertile sites, respectively. Selection based on these parameters, supported by careful analysis of patterns of family performance (see Figs. 1, 2) and stability provided by GGE biplots (see Fig. 4) and regression analysis results (see Fig. 5) is a promising approach to the use of GEI in forest tree breeding. The ASV was useful as stability measure in our study, but its suitability to support selection decisions should be tested.
A possible limitation of our study is the low number of, or rather the small range of variability among environments.
As pointed by Becker and Léon (1988), it is impossible to calculate useful genotype stability measures from a few environments only. However, how many is "a few" in this regard, was not stated. In forestry applications, commonly there are between 2 and 15 locations of experimental sites (e.g., Costa e Silva et al. 2001;Dutkowski et al. 2006;Fu et al. 1999;Schermann et al. 1997;Stoehr et al. 2010). Magnussen (1993) recommends to base the meaningful inferences of GEI on at least five sites. For that matter, the four sites used in our experiment might be at the lower end of the number of environments required for estimation of GEI and stability of tested genotypes. However, our results indicate that even this low number of environments would provide useful information about the extent of GEI for supporting selection decisions in tree breeding programs.

Conclusions
In conclusion, we found a significant GEI for both diameter and height in the progeny test with Scots pine. We used different methods of analysis of GEI and found their variable usefulness for supporting selection decisions. Successful selection decisions need an approach integrating information available from those methods. We propose the following sequence of steps after finding significant GEI in forest progeny tests. (1) Use GGE biplots with the G + GE data centering and different singular value partitioning to identify patterns of interaction. Examine the relationship among environments and the discriminative ability and representativeness of environments. Examine genotypes' performance and stability (contribution to GEI). (2) Perform regression analysis and identify genotype's response to environment (b i ) and stability (s 2 di ).
(3) Compare the patterns of stability based on mean square deviation and GGE biplot. (4) Finally, depending on selection objectives, make selections based on trait value (yield) and environmental responsiveness (b i ) taking stability into consideration.