Interpreting genotype‐by‐environment interaction for biomass production in hybrid poplars under short‐rotation coppice in Mediterranean environments

Understanding genotype × environment interaction (GEI) is crucial to optimize the deployment of clonal material to field conditions in short‐rotation coppice poplar plantations. Hybrid poplars are grown for biomass production under a wide range of climatic and edaphic conditions, but their adaptive performance in Mediterranean areas remains poorly characterized. In this work, site regression (SREG) and factorial regression mixed models are combined to gain insight into the nature and causes underlying GEI for biomass production of hybrid poplar clones. SREG addresses the issue of clonal recommendation in multi‐environment trials through a biplot representation that visually identifies superior genotypes. Factorial regression, alternatively, involves a description of clonal reaction to the environment in terms of physical variables that directly affect productivity. Initially, SREG aided in identifying cross‐over interactions that often involved hybrids of different taxonomic background. Factorial regression then selected latitude, mean temperature of the vegetative period (MTVP) and soil sand content as main site factors responsible for differential clonal adaptation. Genotypic responses depended strongly on taxonomic background: P. deltoides Bartr. ex Marsh. × P. nigra L. clones showed an overall positive sensitivity to increased MTVP and negative sensitivity to increased sand content, whereas the opposite occurred for P. trichocarpa Torr. & Gray × P. deltoides clones; the three‐cross hybrid [(P. deltoides × P. trichocarpa) × P. nigra] often displayed an intermediate performance. This information can contribute toward the identification and biological understanding of adaptive characteristics relevant for poplar breeding in Mediterranean conditions and facilitate clonal recommendation at eco‐regional level.


Introduction
The continued development of sustainable practices for improving short-rotation coppice (SRC) plantations must be grounded, among other factors, on a proper deployment of genetic material as related to planting site characteristics. This goal is often challenged by the presence of genotype-by-environment interaction (GEI) (Yan et al., 2000). GEI is defined as the difference in the response of genotypes to different environments (Bradshaw, 1965). It complicates the identification of superior genotypes, highlighting the need to grow different material in different areas within a broader target region. In tree breeding, however, the implications of GEI for selection purposes are difficult to be considered due to the effort associated with the multi-site evaluation of new breeds. Indeed, GEI evaluation in forest tree species requires time, manpower and economic investment, thus hampering its proper characterization.
Populus species are present across a broad range of climatic and edaphic conditions, bearing an important ecophysiological variability that often underlies interand intraspecific adaptation patterns (Dickmann, 2001). Populus nigra L. (European black poplar) and P. deltoides Bartram ex Marshall (eastern cottonwood), both from the Aigeiros section, and P. trichocarpa Torr. & Gray (balsam cottonwood), from the Tacamahaca section (Eckenwalder, 1984), are probably the most relevant species for cultivation in intensively managed plantations, their hybrids being widely used in the northern hemisphere. The balsam cottonwood grows from sea level to 2000 m and prefers alluvial, fertile soils, although it can thrive with an acceptable performance on loessic soils; conversely, its growth is considerably limited in acidic soils (pH <6) (DeBell, 1990). On the other hand, species belonging to the Aigeiros section are considered to be more plastic: P. nigra is capable of growing on stony, poor and relatively dry soils, while P. deltoides is associated with lowland flood plains and sand riparian corridors with optimum growth on silty or sandy loam soils (Dickmann, 2001). However, P. deltoides can also thrive on deep, infertile sands and on clay soils (Cooper, 1990).
The intraspecific adaptive patterns of these poplar species and their hybrids have been less explored as compared with interspecific characteristics. Interspecific hybrid poplars usually show intermediate morphological and phenological characteristics between those of their parents, but they can also display hybrid vigor and specific adaptation to particular conditions (Bisoffi & Gullberg, 1996). Differences in productivity have been reported according to the genetic background (P. deltoides 9 P. nigra vs. P. trichocarpa 9 P. deltoides) (Marron et al., 2007), as well as in functional traits such as water-use efficiency and stomatal conductance (Dillen et al., 2008;Di Matteo et al., 2015). The identification of the adaptive properties of hybrid poplars for SRC plantations becomes a priority for maximizing biomass production, particularly under the very diverse array of growing conditions typical of the Mediterranean region. For example, Interamerican hybrids (P. trichocarpa 9 P. deltoides) display more vigorous growth than Euramerican hybrids (P. deltoides 9 P. nigra) at higher latitudes (Tabbush & Beaton, 1998). Studies relating productivity to site characteristics are common in forestry, including fast-growing plantations (Corona et al., 1998;Jug et al., 1999;Alvarez et al., 2010;Park et al., 2012). However, little information is available with regard to the effects of site or management factors on potentially contrasting performances of hybrid poplars in Mediterranean areas.
Clones are expected to be more prone to GEI than other breeding material, owing to the absence of genetic variability at the intraclonal level and the potential relevance of nonadditive genetic effects (in addition to additive variance) causing differential clonal performance (Eriksson et al., 2006). Previous studies characterizing the stability of poplar genotypes in multi-environment trials (MET) in the Mediterranean area have revealed an important influence of GEI on productivity. GEI has been reported for both initial growth during the establishment year (Sixto et al., 2011) and total biomass (TB) production (and its different woody fractions) at the end of the first rotation period (Scarascia-Mugnozza et al., 1997;Sixto et al., 2014).
The identification of physical variables (e.g. related to soil, climate or geography) responsible for differential clonal performance across environments would contribute to a better understanding of GEI patterns affecting productivity in SRC poplar plantations. Additionally, the role of site management strategies derived from fertilization, irrigation or weed management practices should also be considered as potential causes for GEI. In this regard, the incorporation of external information describing environments in statistical models provides an effective approach to understanding GEI. In particular, factorial regression models constitute an informative methodology for interpreting GEI (Denis, 1988;Van Eeuwijk et al., 2005) that can be applied to forestry trials in a straightforward manner (Al ıa et al., 1997).
Factorial regression describes GEI as differential genotypic sensitivity to explicit environmental factors (i.e. differences among genotypes in their slope of response). This approach can be extremely instructive to improve the understanding of circumstances under which particular genotypes perform better. Such a strategy can be combined with exploratory GEI approaches that exploit phenotypic information solely for the target trait (e.g. biomass-based information in the case of energy crops). Among these approaches is the sites regression (SREG) model (Crossa & Cornelius, 1997), also known as genotype main effects and genotype 9 environment interaction effects (GGE) model (Yang et al., 2009). SREG addresses the issue of cultivar recommendation in MET trials in a simple manner through a graphical technique: the so-called GGE biplot (Yan et al., 2000) (but see Yang et al., 2009 for a description of its limitations and most common misuses). The combination of SREG and factorial regression models enables not only the identification of superior material, but also the interpretation of key environmental aspects related to differential plant functioning underlying the observed superiority of certain genotypes in particular environments (Voltas et al., 2005).
Most often, the multi-environment evaluation of genotypes calls for a mixed model setting that extends the traditional fixed-effects model, in which genotypes and environments (as well as their interaction) are regarded as fixed. In the last decade, cutting-edge studies in the statistical modelling of GEI have advocated the use of mixed models as the preferred analytic framework, with either genotypes or environments, or both, being treated as random effects (e.g. Smith et al., 2002;Yang, 2007;Burgueño et al., 2011). In this regard, a mixed-model analog of the SREG model has been proposed that uses the factor analytic (FA) model for approximating the variance-covariance GEI structure (Piepho, 1998;Smith et al., 2002). Furthermore, in factorial regression, environmental factors or covariates can be easily embedded into a mixed modelling framework Malosetti et al., 2004).
The objective of this study is to get an insight into the interpretable patterns of GEI of hybrid poplar genotypes differing in their taxonomic background in SRC in Mediterranean environments at the end of the first rotation period. To this end, we combined SREG and factorial regression mixed models to assist in the identification and characterization of superior hybrids in the context of evaluation trials earmarked for biomass production. This work was conceived as an extension of a previous study by Sixto et al. (2014) that was aimed at assessing differences in genotypic stability for productivity traits employing some widely used stability measures: Shukla's stability variance, Finlay-Wilkinson regression and Eberhart-Russell stability analysis. In Sixto et al. (2014), we underlined the feasibility of exploiting specific adaptation for optimizing productivity in Mediterranean conditions, in which case hybrid type would play a relevant role. This realization opened the door to the application of more elaborated models in which additional information can be used in the form of explicit environmental characterization to model GEI for a target trait of choice (Malosetti et al., 2013). Here, we pursue a more comprehensive depiction, description and biological understanding of GEI patterns of SRC poplar plantations grounded on mixed modelling principles, with special emphasis on the taxonomic background of hybrids as regards the adaptive performance of clonal material.

Materials and methods
Trials were conducted at four representative sites of the cropping conditions in which poplars are currently grown in northern Spain. Their characteristics are listed in Table 1. We used nine poplar clones belonging to different interspecific hybrids. Clones were 'I-214', 'MC', 'Guardi', '2000 verde' and 'AF2' (P. deltoides Bartram ex Marshall 9 P. nigra L. = P. canadensis Moench. = P. euramericana (Dode). Guinier), 'Unal' and 'UW 49-177' (P. trichocarpa Torr. & Gray 9 P. deltoides = P. generosa Henry = P. interamericana Brockh), and 'Monviso' and 'Pegaso' [(P. deltoides 9 P. trichocarpa) 9 P. nigra]. The plant material used in the experimental sites was produced in a stool bed established for the provision of stock plant for research purposes. Trials were planted out manually in spring using cuttings at a density of 33 333 cuttings ha À1 , spaced 1 9 0.30 m, and a pre-emergence herbicide for weed control (Oxyfluorfen, 4 l ha À1 ) was applied. Three different treatments were tested at each site: (i) Additional weed control with herbicide (W): glyphosate application using a protective cover, (ii) Supplementary fertilization (F): granulated fertilizer containing N:P:K (12:12:17) and trace elements (0.1% Fe, 0.02% B and 0.01% Zn) applied at a rate of 800 kg ha À1 , and (iii) Control, with no further treatment (C). The trials were irrigated up to field capacity using a drip system to optimize the watering application. Irrigation is a common practice in SRC plantations in many Mediterranean areas because rainfall is not enough to optimize tree development.
A hierarchical split-plot experimental design was employed with four replicates per clone and treatment at each location. Two fixed factors were considered: genotype (nine clones, main plot) and treatment (three types, subplot). The allocation of treatments to subplots resulted in three distinct subenvironments for each of the four sites, so a total of 12 evaluation environments were considered in which each genotype had a particular response. Each replicate subplot consisted of 25 ramets of the same clone, resulting in square experimental units in which the nine central ramets were monitored. Yield was evaluated by recording total above-ground dry biomass per plant (TB, g) at the end of the first rotation (3-year-old plants). Data were expressed as dry weight after estimating the humidity content from a subsample of each plot, which was oven-dried to constant weight at 100°C. Detailed information on the series of trials and the experimental design can be found in Sixto et al. (2014).

Site characterization
The environmental characterization of each site was performed by considering physical variables that were likely to have a direct influence on above ground biomass in SRC poplar plantations. These comprised (i)

Statistical analysis
Combined MET analysis. The four experiments were combined across locations and analyzed as multi-environment analysis of variance according to the following model: where Y iklrp is the observation of the ith genotype and lth agronomic treatment in the kth location and rth block for the pth tree, l is the general mean, L k is the effect of the kth location, (LR) kr is the random effect of the rth block nested to the kth location, G i is the effect of the ith genotype, (GL) ik is the effect of interaction between the ith genotype and the kth location, (GLR) ikr is the random effect of the interaction between the rth block nested to the kth location and the ith genotype, T l is the effect of the lth agronomic treatment, (LT) kl is the effect of interaction between the kth location and the lth agronomic treatment, (GT) il is the effect of interaction between the ith genotype and the lth agronomic treatment, (GLT) ikl is the effect of interaction between the ith genotype, the kth location and the lth agronomic treatment, e iklr is the random residual effect of the ith genotype and lth agronomic treatment in the rth block of the kth location, and p ijkl is the random tree effect of the lth individual in the iklrth plot. Prior to the combined analysis across locations, error variances were tested for homogeneity using Levene's test.
Mixed modelling of MET data. Mixed models accounting for interaction and heteroscedasticity in two-way genotype-byenvironment tables were fitted complementing a number of stability measures, as reported in Sixto et al. (2014): Shukla's stability variance, Finlay-Wilkinson regression and Eberhart-Russell stability analysis. Model fitting was performed on a genotype-by-environment table of means. This implies that a two-stage strategy for analyzing MET data was applied (for details, see e.g. Malosetti et al., 2013;M€ ohring & Piepho, 2009). First, trials were individually analyzed and least square estimates of genotype-treatment means and plot error variances were obtained at each trial for TB. From these individual trial analyses, adjusted means were carried forward to the second stage, where different models were fitted to genotype-by-environment means. At this stage, each combination of location and agronomic treatment was regarded as a different 'environment'. We did not consider the use of weights (e.g. reciprocals of plot error variances) for estimates of genotype-treatment means, as the low range of plot error variances pointed to approximately equal precision of estimates at the trial level (largest error variance = 445 200 g 2 , lowest error variance = 91 404 g 2 ; approximate acceptable limit for a pooled residual = 10-fold range in error variances, current range = 4.8) (Williams et al., 2002). Single-stage analyses have certain theoretical advantages over two-stage analyses, but two-stage analyses are logistically and computationally easier to handle in complex mixed models (Malosetti et al., 2013). The standard additive two-way mixed model fitted to genotype-by-environment means can be expressed as: where Y ij is the observation of the ith genotype in the jth environment, l is the general mean, G i is the effect of the ith genotype, E j is the effect of the jth environment, and (GE) ij is the effect of interaction between the ith genotype and the jth environment. For our purpose of testing cultivars for a target agroecological condition (i.e. the irrigated cropping system for poplar cultivation in the western Mediterranean), we regarded the different environments as random. Therefore, the random effects E j and (GE) ij are assumed to be independently distributed with zero mean and variances r 2 E and r 2 GE . Note that the main effects of location (L k ) and agronomic treatment (T l ) and their interaction (LT) kl are not explicitly fitted in the standard model (2), being subsumed as part of the Environment factor (i.e. E j = L k + T l + LT kl ). The additive mixed model provided a starting point for comparison with other more complex but informative models: the SREG mixed model and the family of factorial regression mixed models. Each alternative model is outlined here as the sum of three components: the fixed terms, the random terms and the residual term. The residual term comprises the interaction not yet accounted for by the existing fixed and random terms. The residual term is very flexible for modelling purposes, and classical GEI approaches can be handled using appropriate variance-covariance (VCOV) structures for its approximation (Yang et al., 2009;Sixto et al., 2014). Following Denis et al. (1997), the general forms of expectation (E) and VCOV structures (Var) read as follows: where g is a column vector of size I and Γ is any covariance matrix of size I. Three different types of Γ were used here: Simple diagonal (R 1 ): GE is a common variance of GE effects across genotypes.
Diagonal (R 2 ): GEi is a genotype-dependent variance of GE effects for the ith genotype.
Site regression mixed model. The SREG mixed model is a mixed version of the GGE model (also known as SREG model) that uses a FA structure with two latent factors (FA(2)) to approximate the residual VCOV structure of (3) (Yang et al., 2009). Here, the residual VCOV structure models the main effects of genotypes plus GEI, while the environmental variance component r 2 E is included as a random term in the general model (3), as follows: This model, a mixed-model analog of the GGE biplot analysis (Yan & Kang, 2003), has environmental loadings and genotypic scores of the two latent components of FA(2) that are interpreted in a similar way to those for the SREG 2 fixed-effect model (Yang et al., 2009). Details on the identifiability of FA(2) and the principal component rotation applied to this dataset are provided in Burgueño et al. (2008). The number of variance components is 2(J + 1). The most common use of GGE biplot analysis in evaluation trials is the identification of mega-environments and the visualization of 'which-won-where' patterns, thereby summarizing the information contained in complex datasets. 'Which-won-where' patterns are only really identifiable if the correlation between the genotypic scores of the first latent factor and the genotypic main effects is positive and nearly perfect (Yan et al., 2001).
Factorial regression mixed models. In factorial regression models, linear regression terms for the explanation of GEI are incorporated in the form of environmental covariates to the levels of the environmental factor. In this way, physical variables underlying GEI can be identified. Three variants were fitted for this purpose.
Model 1 (mixed factorial regression model). This is an extension of the fixed factorial regression model as reviewed in Van Eeuwijk (1995). In factorial regression, explicit environmental information is included in the levels of the environmental factor to describe the interaction term, GE: where g is a I 9 2 matrix of fixed parameters containing the ith genotype main effect and the ith genotypic sensitivity to the environmental factor z j for the jth environment, and R 1 is defined above. For the sake of simplicity, we used here a very restricted set of environmental covariates ( Table 1) that were tested independently, given the much reduced number of available locations. It is important to note that the term g i2 z j of the product matrix g(1, z j )' becomes fixed because z j values contain known, explicit environmental information rather than unobserved variables, as e.g. for the Finlay-Wilkinson model (Piepho, 1997). The number of variance components equals two.
Model 2 (Shukla's mixed factorial regression model). This is a mixture between a general heteroscedastic model (model 2 in Sixto et al. (2014)) and a mixed factorial regression model, and was first proposed by Shukla (1972) in its fixed guise: Here, the information about the genotypes is concentrated into 'triplets' of parameters: a general level of performance (g i1 ), a measure of sensitivity to the environmental covariate (g i2 ) and a stability variance (r 2 GEi ). The diagonal R 2 structure is defined as above, and the number of variance components equals I + 1 .
Model 3 (Shukla's structured mixed factorial regression model). This model is a simplified version of Shukla's mixed model that allows for the definition of G different groups of genotypes with a priori different variability (i.e. heteroscedasticity). Here, three different groups of genotypes (i.e. taxonomic groups) were considered: P. deltoides 9 P. nigra, P. trichocarpa 9 P. deltoides, and [(P. deltoides 9 P. trichocarpa) 9 P. nigra]. The expectation and VCOV structures are essentially the same as in (6), but the diagonal R 2 structure is defined so as to assign the same residual variance to all genotypes within a particular group. The number of variance components is G + 1.
The best fitting environmental covariate was independently selected for the set of geographic, bioclimatic and edaphic factors under model 1. Next, the adequacy of the VCOV models incorporating stability variances at the genotype (model 2) or genotype group (model 3) level was tested against model 1 by computing the restricted log-likelihood statistics for each model and deriving information criteria such as Akaike's information criterion (AIC) and Bayesian information criterion (BIC). Both involve a penalty for the number of parameters in the VCOV structure, which favors parsimonious models, but BIC generally penalizes a large number of parameters more strongly than does AIC. Both statistics are in the smaller-is-better form. It should be noted that the environmental covariates have only four levels characterizing each of the four different locations, while the number of total environments (or location-agronomic treatment combinations) is 12. By using the complete genotypeenvironment matrix (instead of a simplified genotype-location dataset), we aimed at improving the estimation of stability variances for genotypes (model 2) or genotype groups (model 3); on the other hand, estimates of genotypic sensitivities are equivalent regardless of the chosen approach. The analyses were performed with the MIXED procedure of SAS/STAT v.9.3 (2004), using restricted maximum likelihood for estimation of variance components.

Multi-environment trial analysis
The MET analysis of variance for TB revealed significant differences among locations (L) and among agronomic treatments (T), and also a significant L 9 T interaction (Table 2). Overall, the highest biomass was achieved by the fertilization treatment (1237 g per plant), followed by the control and herbicide treatments (1085 and 1068 g per plant respectively). There were significant differences among clones (G) and also a significant G 9 L interaction. It is important to stress that neither the interaction between agronomic treatment and clone nor the second-order G 9 L 9 T interaction was relevant for TB (Table 2). Therefore, GEI patterns were mainly due to a differential reaction of clones to the testing sites.

Site regression mixed model
The rotated to principal component biplot of the FA(2) model displays the response patterns of nine hybrid poplar clones evaluated in 12 environments for TB (Fig. 1). The amount of GGE variation explained by the two components of the FA(2) model was 75% (first component = 55%; second component = 20%), suggesting that a rank-two approximation is adequate for defining mega-environments and evaluating genotypes. In particular, the identification of mega-environments and depiction of 'which-won-where' patterns in the biplot seem accomplished for this dataset, as the correlation between genotypic scores of the first latent factor and the genotype main effects (derived from Eqn 2) was highly positive (r = 0.92).
The environmental loadings of the first latent factor were mostly positive or close to zero with the exception of S4H and S4T, which exhibited negative loadings. The twelve environments fell into three apparent groups (i.e. mega-environments): most environments from sites 1 and 3 formed one group (hereafter EG1), most environments from site 2 formed a second group (EG2) and, finally, environments from site 4 formed a third group (EG3). As for genotypes, two groups could be visually identified by the scores of the first latent factor: 'UW 49/177', 'Pegaso' and 'Unal' (with negative scores) and 'AF2', '2000 verde' and 'MC' (with positive scores). On the other hand, clones 'Guardi', 'Monviso' and 'I-214' displayed close-to-zero scores for the first latent factor. 'Monviso' and '2000 verde' showed the most extreme scores (positive and negative respectively) for the second latent factor.
In a vector representation, positions of environments on the biplot plane determine lines starting at the origin (0,0), with their corresponding vector norms associated with GEI variability. Thus, it was feasible to visualize the approximate ranking of genotypes in a particular environment on the basis of their projections ordered on the environment vector. For example, the biplot suggested '2000 verde' as the best yielding clone at S3H, while the least productive clones were 'Pegaso' and 'UW49/177'. By using the inner-product property of the biplot, the SREG mixed model was useful for visually identifying cross-over interactions: clones 'AF2' and '2000 verde' had higher TB in EG1, whereas 'Monviso' had higher TB in EG2 and 'Pegaso' displayed better performance in EG3. In the latter case, however, the environments included in EG3 had close-to-zero loadings, †For fixed effects. indicating that the superiority of 'Pegaso' should be interpreted cautiously.

Site characteristics underlying GEI
Single-covariate factorial regression mixed models (model 1; Eqn 5) were first fitted to the two-way table of genotype-environment means for TB (Table 3). They provided an initial screening of relevant physical variables underlying genotypic responses to changing growing conditions. Significant geographic covariates were LAT, VP and ALT (in this order of significance), whereas all climatic factors were relevant for the explanation of differential genotypic responses (with MTVP as the most important). Also, all edaphic factors were found to underlie GEI patterns, with SAND and pH (in this order) being the most relevant ones ( Table 3). The best fitting variable of each class was used thereafter for testing two variants of factorial regression: Shukla's mixed factorial regression model (model 2; Eqn 6) and Shukla's structured mixed factorial regression model (model 3). Models 2 and 3 account for heteroscedasticity in the residual effects so different stability variances can be added for each clone (model 2) or taxonomic group (model 3). Model 2 had the most satisfactory fit to the residual variance structure in all cases according to AIC and BIC statistics (Table 4). Following (6), the expectation of any genotype-environment combination in model 2 takes into account both the genotypic mean and the sensitivity of the genotype to a selected physical variable (either LAT, MTVP or SAND) given by the factorial regression slope (g i2 ). In our case, clones 'Pegaso' and 'Unal' were sensitive to changes in all three variables (LAT, MTVP and SAND), whereas 'Guardi', 'I-214', 'MC' and 'Monviso' did not respond significantly to any particular variable. Other clones showed specific responses to LAT ('2000 verde') or to both MTVP and SAND ('AF2' and 'UW49/177') ( Fig. 2).
Shukla's mixed factorial regression models also explained genotypic stability in terms of independent clonal variances of residual GEI effects, but a distinctive (i.e. structured) taxonomic pattern was not detected here, irrespective of the external variable entering the model (Table 4). The relationship between genotypic stabilities estimated for the different Shukla's models (i.e. involving LAT, MTVP or SAND) was high (Spearman's rank correlation between stability estimates ≥0.82; P < 0.01). Clone '2000 verde' consistently showed the largest residual stability variance and 'Pegaso' and 'UW49/177' had the lowest values for this parameter (results not shown).

Discussion
The biomass produced by hybrid poplars in this MET depended simultaneously on factors associated with the plantation site and the genotype, but also on specific genotypic responses to particular site conditions, hence pointing to the existence of GEI (Sixto et al., 2014). Biplots are effective tools for visually identifying GEI patterns, allowing environments with similar characteristics and genotypes with comparable performance to be identified (Yan et al., 2000), but they are useful for recommendation purposes only if the target region is sufficiently sampled (Yan et al., 2001). Additionally, they should be considered an initial, exploratory step toward a comprehensive understanding of GEI (Yang et al., 2009). To gain predictability in GEI, factorial regression can produce more parsimonious models that help determine the underlying physical factors for In bold, most relevant environmental variables (P < 0.05) underlying GE interaction according to a chi-squared distribution (asymptotic approximation). interaction (e.g. Voltas et al., 2005). Here, we have illustrated how SREG and factorial regression mixed models can be combined using a common methodological framework to gain an insight into GEI for biomass production from hybrid poplar clones. This strategy is more in line with the recognized need to allow for random effects in the analysis of MET trials (e.g. Yang, 2007). The biplot analysis suggested a stable performance of certain clones, namely 'I-214', 'Guardi' and, to a lesser extent, 'MC', the latter displaying the highest production of this group. Selecting the most appropriate genotype is one of the main concerns of plantation managers when planning short-rotation forest crops; thus, the availability of such generalist genotypes can be particularly useful in the absence of precise site information. The remaining clones exhibited a clearer preference for certain environments in which their production was maximized, revealing a more site-specific behavior. For example, the biplot identified 'AF2' as a highly productive clone over the MET trial (with the exception of S4 environments), but its performance seemed to improve in those environments where additional fertilizer was applied (cf. Fig. 1). Variation in the response of poplar genotypes to treatments involving different degrees of fertilization has been previously reported (Kara ci c & Weih, 2006;Zalesny et al., 2007). Such variation may be attributed, among other factors, to differing behavior as regards nitrogen economy. However, perhaps with the exception of 'AF2', a differential clonal response to fertilization could not be observed in this MET trial. Similarly, only 'Monviso' seemed to perform relatively better in environments where more rigorous weed control was applied (S2H and S1H), or where weeds were already scarce prior to the application of additional herbicide, such as in S2C (visual observation). But, despite such approximate interpretation of GEI patterns, we did not detect a significant treatment by genotype interaction for TB in our dataset (Table 2). Makeschin (1999) stated that aspects related to management and, particularly, to weed control have an impact on the production of biomass in short-rotation plantations in central and northern Europe. Other studies have addressed weed-poplar competition dynamics (Otto et al., 2010), although no specific mention has been made, to the best of our knowledge, as regards differences in the productive performance of genotypes associated with competition from weeds.
In relation to the evaluation sites, S4 (with environments S4H, S4C and S4F), which presented the most extreme edaphoclimatic conditions for growth, was the least discriminant site for clonal differences as regards productive potential. This fact had already been noted for the establishment year when growth variables for the same set of genotypes were evaluated (Sixto et al., 2011). Conversely, S3 proved to be a good location, allowing for greater discrimination among genotypes as average yields were high (first latent score) and interaction was low (second latent score), regardless of the agronomic treatment applied.
Among the site characteristics related to geographic, climatic or soil factors, LAT, MTVP and SAND were determined to be most relevant for the explanation of clonal differences in biomass production. Soil texture has been identified as one of the most important site factors for production in Salicaceae (Labrecque & Teodorescu, 2003;Pinno & B elanger, 2009;Bergante et al., 2010). Mean temperature during the growing season has also been identified as a major climatic factor for Salix species (Labrecque & Teodorescu, 2003). Conversely, site altitude, which had previously been identified as highly explanatory of variation in clonal growth (Tabbush & Beaton, 1998), was not a very relevant (albeit significant) factor in this study despite strong differences in this variable among sites. In fact, the effect of altitude is partly accounted for by the observed differences in growing season length among sites, which turned to be a more important factor than altitude for GEI explanation. The combined effect of the length of the growing period (VP) and altitude is actually reflected in MTVP, the most relevant environmental variable underlying GEI.
Those clones identified as stable according to the biplot ('I-214', 'Guardi' and 'MC'), together with the 'Monviso' clone, showed no sensitivity to the aforementioned environmental covariates. However, the remaining clones exhibited a significant sensitivity to at least one covariate. The low productive clones 'Pegaso' and 'Unal' were sensitive to all variables identified as relevant in the analysis, benefiting from cool vegetative periods, sandy soils and high-latitude sites. On the other hand, the Euramerican clone '2000 verde' behaved better at low latitudes. 'Pegaso' and 'Unal', hybrid clones from [(P. deltoides 9 P. trichocarpa) 9 P. nigra] and (P. deltoides 9 P. trichocarpa), respectively, are better adapted to higher latitudes and, therefore, to shorter vegetative periods than the Euramerican hybrids. For example, the Euramerican clone 'AF2' benefited from warm vegetative periods and also from not overly sandy soils. This performance has been highlighted in a recent study (Di Matteo et al., 2015). An opposite performance was detected for the Interamerican clone 'UW49-177', which preferred cooler vegetative periods and sandy soils. The overall site productivity appears also to be influenced by the covariates identified as significant for GEI. In particular, the least productive site (S4) had the lowest MTVP and the largest percentage of sand. Conversely, the most productive and interactive site (S3) showed the highest values for these variables and also had the most southerly latitude. The remaining sites (S1 and S2) presented intermediate values. We acknowledge that four sites may not suffice to capture conveniently the diversity of agroecological conditions for SRC poplar plantations encountered in northern Spain, and that other potentially relevant characteristics for GEI such as soil fertility (e.g. N, P and K content) have not been considered in this study in the absence of a proper soil characterization. Moreover, we have not considered features related to soil water availability which, according to Bergante et al. (2010), is the most important factor explaining variation in biomass production under Mediterranean conditions. This factor, however, was considered to be of little relevance, given that plantations were irrigated during the entire vegetative period and the soil kept close to field capacity, which is a usual management practice in southern Mediterranean areas.
In summary, genotypic responses to environmental covariates depended strongly on taxonomic background. P. deltoides 9 P. nigra clones showed an overall positive sensitivity to increased MTVP and negative sensitivity to increased SAND, whereas the opposite occurred for P. trichocarpa 9 P. deltoides clones. The three-cross hybrid [(P. deltoides 9 P. trichocarpa) 9 P. nigra] often displayed an intermediate performance, which may be due to its intermediate genetic background in relation to the other hybrid groups. However, in the case of LAT, it was not possible to identify such a clear response pattern based on taxonomic affiliation, even though it was expected that the Interamerican hybrids would exhibit an overall positive sensitivity to increased latitude of the testing site. In this regard, the response of genotypes to latitude being dependent on the geographic origin of the Fig. 3 Genotypic sensitivities to environmental factors as related to Finlay-Wilkinson empirical sensitivities to improved environmental conditions. Empty circles refer to P. deltoides 9 P. nigra, half-filled circles refer to (P. deltoides 9 P. trichocarpa) 9 P. nigra, and filled circles refer to P. trichocarpa 9 P. deltoides. parental material has been accepted for decades (Pryor & Willing, 1965). Here, the relatively low number of testing sites and clones may have precluded the detection of clearer responses to latitude changes. In any case, the environmental characteristics linked to GEI point to the existence of different adaptive patterns depending on the taxonomic background, probably influenced by the particular geographic origin of the parental material.
Although only a limited number of site characteristics were analyzed, clonal performance was found to vary according to geographic (latitude) and edaphoclimatic gradients which explained, at least partially, GEI patterns in clones previously characterized as exhibiting specific adaptation (Sixto et al., 2014). This is confirmed by the close associations between the Finlay-Wilkinson stability parameter for each clone (Sixto et al., 2014) and the genotypic sensitivities to environmental factors obtained from Shukla's mixed factorial regression models (Fig. 3). The overall response patterns detected at taxonomic group level suggest differences in adaptive performance for each genetic background of hybrid poplars. These patterns could help facilitate clonal recommendation in the Mediterranean area as well as highlight adaptive characteristics to be considered in breeding programs. Although further research and more extensive testing is needed in relation to the factors that influence clonal response and contribute to GEI, the information provided in this study may contribute toward improving GEI predictability in hybrid poplars at Mediterranean eco-regional level.