Genetics of phenotypic plasticity and biomass traits in hybrid willows across contrasting environments and years

Abstract Background and Aims Phenotypic plasticity can affect the geographical distribution of taxa and greatly impact the productivity of crops across contrasting and variable environments. The main objectives of this study were to identify genotype–phenotype associations in key biomass and phenology traits and the strength of phenotypic plasticity of these traits in a short-rotation coppice willow population across multiple years and contrasting environments to facilitate marker-assisted selection for these traits. Methods A hybrid Salix viminalis × (S. viminalis × Salix schwerinii) population with 463 individuals was clonally propagated and planted in three common garden experiments comprising one climatic contrast between Sweden and Italy and one water availability contrast in Italy. Several key phenotypic traits were measured and phenotypic plasticity was estimated as the trait value difference between experiments. Quantitative trait locus (QTL) mapping analyses were conducted using a dense linkage map and phenotypic effects of S. schwerinii haplotypes derived from detected QTL were assessed. Key Results Across the climatic contrast, clone predictor correlations for biomass traits were low and few common biomass QTL were detected. This indicates that the genetic regulation of biomass traits was sensitive to environmental variation. Biomass QTL were, however, frequently shared across years and across the water availability contrast. Phenology QTL were generally shared between all experiments. Substantial phenotypic plasticity was found among the hybrid offspring, that to a large extent had a genetic origin. Individuals carrying influential S. schwerinii haplotypes generally performed well in Sweden but less well in Italy in terms of biomass production. Conclusions The results indicate that specific genetic elements of S. schwerinii are more suited to Swedish conditions than to those of Italy. Therefore, selection should preferably be conducted separately for such environments in order to maximize biomass production in admixed S. viminalis × S. schwerinii populations.


INTRODUCTION
In heterogeneous environments, including those exposing plants to frequent climate fluctuations, individual plants' fitness and productivity will depend on their ability to phenotypically adapt to rapidly shifting environments. In such environments, phenotypic plasticity, which is the ability of an organism to alter the expression of a phenotypic trait as a response to changes in environmental conditions (Bradshaw, 1965;Schlichting, 1986) can provide increased environmental tolerance. Therefore, populations composed of highly plastic individuals are expected to be able to inhabit a broad range of environmental conditions or be able to withstand large environmental fluctuations. As a result, phenotypic plasticity is expected to influence the geographical distribution of species (Agrawal, 2001;Sexton et al., 2009).
Phenotypic plasticity of an individual can be depicted by its norm of reaction, which shows the variation of phenotypic traits across different environments, thus displaying the environmental sensitivity of that individual (Schmalhausen, 1949). Individuals of a population can indeed vary in how sensitive they are to environmental changes, which is evident as significant genotype Â environment (G Â E) interactions, where the extent to which a trait value changes across environments is the phenotypic plasticity (Via and Lande, 1985;Wu, 1998).
As phenotypic plasticity of different traits can influence the performance of individuals across environments, it should be taken into account when conducting selection and breeding for variable climates and environments (Via and Lande, 1985). In order to do this, it is crucial to characterize the level of phenotypic plasticity and its variation among individuals used for breeding and to determine how much of the plasticity has a V C The Author 2017. Published by Oxford University Press on behalf of the Annals of Botany Company. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. genetic basis (de Jong and Bijma, 2002). Organisms that can be clonally propagated, e.g. poplars, are excellent study systems as phenotyping can be carried out in clonal replicates grown in contrasting environments (Wu, 1998;Wu and Hinckley, 2001;Fabbrini et al., 2012). Interestingly, clonally replicated poplar populations grown across climatic gradients and contrasting environments have revealed large variation in phenotypic plasticity between clones (Wu andStettler, 1997, 1998;Rae et al., 2008;Guet et al., 2015;Elferjani et al., 2016). Furthermore, a genetic basis of phenotypic plasticity has also been demonstrated for important biomass traits (Rae et al., 2009;Fabbrini et al., 2012). Willows in the Salix genus are closely related to the poplars and in temperate regions are grown as short-rotation coppice (SRC) primarily for bioenergy production (Kuzovkina et al., 2008). However, improving biomass yield and resistance to pathogens and pests is essential in order to be able to expand willow cultivations to facilitate commercialization. Breeding strategies, especially for perennial crops and trees such as willows, are highly dependent on the level of phenotypic plasticity since selection could be based on high performance across diverse environments aiming for one wide breeding zone or based on high performance in specific environments aiming at different breeding zones.
Salix viminalis and Salix schwerinii are two closely related willow species (Berlin et al., 2011;Fogelqvist et al., 2015) that are commonly used when breeding for fast-growing and resistant SRC willow cultivars. Hybrid crosses between the two species are easily done experimentally, and many commercial cultivars contain a mix of S. viminalis and S. schwerinii genetic backgrounds. This is mainly due to a major genetic factor originating from S. schwerinii that confers resistance to the leaf rust fungus Melampsora larici-epitea, a pathogen that leads to significant yield losses in willow cultivations (Samils et al., 2011). Thus, breeding programmes using both species have successfully developed hybrid cultivars with high biomass productivity and strong resistance to the leaf rust, and that are adapted to the conditions found in Sweden and northern Europe. It is, however, not yet known how hybrid populations are performing in contrasting environments, e.g. Sweden versus southern Europe, and/or in conditions with good water availability versus drought conditions. In particular, it is poorly known whether mixed populations are phenotypically plastic in this regard or whether they harbour any genetic variation for phenotypic plasticity. Addressing such issues becomes increasingly important when adapting breeding populations of long-lived perennial species to future possible climate changes (Alberto et al., 2013).
Keeping this in mind, it should be noted that the identification of genotype-phenotype associations may provide valuable information about genetic regulation of key breeding traits as well as phenotypic plasticity (Nicotra et al., 2010) and may furthermore directly assist breeding for these traits by the use of marker-assisted selection (MAS) (Lande and Thompson, 1990) and genomic selection (Goddard and Hayes, 2007). Genotypephenotype associations based on the segregation of alleles in bi-parental families can be identified by quantitative trait locus (QTL) mapping analyses using genetic linkage maps. This method can be particularly attractive when investigating largely unexplored traits in non-model organisms with limited genomic resources, as the number of markers required to cover the genome is modest and no a priori information about the underlying genes or their function is necessary. QTL mapping analyses using mixed populations have previously been done for biomass , Tsarouhas et al., 2002 and phenology traits (Tsarouhas et al., 2003;Ghelardini et al., 2014), as well as for rust resistance (Samils et al., 2011) and drought tolerance Pucholt et al., 2015).
The primary goal of the present study was therefore to study genotype-phenotype associations in key biomass and phenology traits as well as the plasticity of these traits in a hybrid willow pedigree population (called S 1 ) by QTL mapping analyses. The S 1 population was produced by crossing an S. viminalis female ('78183') with a S. schwerinii ('79069') Â S. viminalis ('Orm') hybrid male (cultivar 'Björn') (Berlin et al., 2010;Tsarouhas et al., 2002) and the phenotypic assessment of the S 1 progeny was done over multiple years in three common garden experiments located in contrasting environments: one in Sweden and two in Italy (one irrigated and one non-irrigated). The S. schwerinii accession 79069 is regarded as a key accession and has been frequently used as an S. schwerinii allele donor, usually via the hybrid cultivar 'Björn', in willow breeding programmes in northern Europe (Karp et al., 2011). Results from this study are thus essential for identifying markers that are associated with trait variation across a wide range of climates and conditions in order to implement MAS for highyielding SRC willow cultivars for diverse environments and a changing climate.

Plant material, common garden experiments and phenotypic assessment
In this study, phenotypic and genotypic data from 463 progenies from the S 1 pedigree population were analysed (Tsarouhas et al., 2002;Berlin et al., 2010). The population was produced by crossing the S. viminalis female '78183' and the S. viminalis Â S. schwerinii hybrid male 'Björn' and is maintained in a plant archive in Pustn€ as, South of Uppsala, Sweden (59 49 0 N, 17 40 0 E). For the study, three different common garden experiments were established (from now on called experiments): one located in Pustn€ as, Sweden, and two in Cavallermaggiore, Italy (44 43 0 N 7 41 0 E). The climatic conditions differed markedly between the two geographical regions. Cavallermaggiore has a higher average annual temperature (12Á8 C) in comparison with Uppsala (5Á7 C) and Cavallermaggiore also had a higher precipitation (746 mm) than Uppsala (551 mm) (www.climate-data.org).
In spring 2008 the experiment in Pustn€ as was established, employing a complete randomized block design with six blocks. To establish the experiment, each of the progeny in the S 1 population was clonally propagated by stem cuttings. The cuttings were planted in pots and after 5 weeks of growth in the greenhouse they were planted in the field. One clonal replicate per progeny was planted in every block, so that each progeny was replicated six times (from now on called a 'clone'). The clonal replicates were randomly positioned in every block. The experiment was manually weeded every year and fertilizer corresponding to 80 kg of N per ha was applied after each cutdown/harvest. The plants were cut down in 2009, 2010, 2011 and 2014 (Table 1, Supplementary Data Table S1). As a non-destructive measurement of biomass, the summed basal area (SumBA) was calculated for each plant based on the diameter of each shoot (Nordh and Verwijst, 2004). After the harvests in 2010 and 2014, the fresh weight (FW) of each plant was assessed. These fresh weights thus reflect 2 and 3 years of stem growth, respectively. In addition, two phenology traits were measured: leaf senescence (LS) in 2013, according to the scale in Ghelardini et al. (2014) and bud burst (BB) in 2014 with the scale used in Weih (2009). Data on LS from 2010 from Ghelardini et al. (2014) were added to the analyses for comparison.
In spring 2012, two experiments with contrasting water regimes were established in Cavallermaggiore: one irrigated (Cavallermaggiore-IRR) and one non-irrigated (Cavallermaggiore-NI). To establish the experiments, 12 stem cuttings of 20 cm each were taken from each clone in the Pustn€ as experiment in winter 2012. The same experimental layout as in Pustn€ as was employed, although the distance between rows was 140 cm (compared with 130 cm in the Pustn€ as experiment). One block in the Cavallermaggiore-NI experiment was omitted due to poor survival of the plants. The cuttings were planted on 2 May 2012 in a well-cultivated soil. During the establishment year mechanical weed control was done repeatedly between rows. No fertilizer was applied. Biomass and phenology traits were assessed on each plant in 2013 and 2014, and Nsh was counted and D measured in 2013 and in 2014. We assessed LS and BB in 2013 and 2014, respectively. A description of all measurements in every experiment is presented in Table S1.

Statistical analyses
To obtain trait predictors for each clone in every experiment, phenotypic data were subjected to variance decomposition by applying mixed linear models with the statistics software ASReml 3Á0 (Gilmour et al., 2009). The distributions of the SumBA and FW traits deviated considerably from normal and therefore additional analyses were performed where these traits were transformed by the square root. Nonetheless, because data transformation was observed to have a negligible impact on the results, only the untransformed data analysis results are henceforth treated. The mixed linear model applied was: where y ijkl is the phenotypic trait value in the ith block for the jth clone located at row k and position l. The overall mean is denoted as l, while b signifies the fixed block effect, c the random clone effect, r Â p the random interaction between rows and positions (spatial term) and e the random residual. Random effects were obtained using the best linear unbiased predictor (BLUP) approach (Robinson, 1991) and were assumed to be independent except for the spatial term (r Â p), which was restricted to follow a two-dimensional first-order autoregressive correlation structure (Cullis and Gleeson, 1991) across plant rows and positions. Thus, clone predictors (y pr ) adjusted for block, and spatial environmental effects were obtained from the effect estimates as y pr;j ¼ b l þ b c j and used in subsequent analyses.
Because the clone predictors are based on random terms, the broad-sense heritability of the clone predictors (H 2 pr ) can be regarded as equivalent to the square of the predictor accuracy (sometimes called r TI ; Henderson, 1984). Thus, following the approach of Wei et al. (2003), H 2 pr for each trait was calculated as: where (PE j ) is the prediction error of clone j supplied by ASReml for each of the n tested clones, while r 2 c is the clone variance estimated by eqn (1).
In order to assess the extent of G Â E interactions that truly affect clonal ranking and that are not due to trait scale  (Burdon, 1977). As the clone predictors largely reflect genotypic trait values, such correlations should closely resemble the corresponding genotypic correlations between experiments. Correlation estimates close to zero would thus indicate considerable crossover-type G Â E interactions, while estimates close to unity would indicate no or little G Â E interaction. Correlations and their 95 % confidence intervals were obtained by using the cor.test function in R (R Core Team, 2013). Plots of clone predictors between contrasting experiments were made with the program JMP V R 10 (SAS Institute Inc., 2010).
To get an estimate of the phenotypic plasticity for each clone, a plasticity trait variable was constructed based on standardized (mean ¼ 0, variance ¼ 1) BLUP values by taking the difference between BLUP values in contrasting experiments. Plasticity traits were constructed in this way for BB14, LS13, Nsh13, MeanD13 and SumBA13 between (1) Pustn€ as and Cavallermaggiore-IRR, (2) Pustn€ as and Cavallermaggiore-NI and (3) Cavallermaggiore-IRR and Cavallermaggiore-NI. A plasticity trait variable for SumBA14 was constructed between Cavallermaggiore-IRR and Cavallermaggiore-NI. In addition, to encompass between-year variation at Pustn€ as, plasticity trait variables were constructed for (1)

QTL mapping analyses
The linkage map developed for the S 1 pedigree population (Tsarouhas et al., 2002;Berlin et al., 2010) was used in QTL mapping analyses. For this study, 41 new markers were added to the map with the JoinMap software (Van Ooijen, 2009). Nineteen markers were positioned on linkage group (LG) XV, one on LG XII, 19 on LG XIX and two on LG XIII. The new markers were developed and genotyped as described in Berlin et al. (2010). In total, the linkage map contains 696 markers with a mean distance of 4.4 cM between the markers. QTL mapping was performed using the clone predictors for each trait in each environment and also on the plasticity traits. The program MapQTL 6.0 (Van Ooijen, 2009) with interval mapping and a regression model was used, where the genome was scanned at 1-cM intervals to determine putative QTL associated with the variation of each trait. In order to determine significant QTL, significance threshold values, estimated as the logarithm of the odds ratio (LOD), were determined with a permutation test of 1000 repetitions. A genome-wide threshold for a significant QTL was set at 0Á05. One and two LOD confidence intervals for each QTL were estimated using the LOD drop-off method (Lander and Botstein, 1989) based on the LOD value at the peak position of the QTL. The proportion of the clone predictor variation explained by each significant QTL was estimated (PVE %). The difference between maternal alleles was estimated as the absolute effect of: (3) the differences between paternal alleles were estimated as the absolute effect of: (4) and the paternal-maternal interaction effect was estimated as the absolute effect of: where l ac , l bc , l ad and l bd are the estimated phenotypic means of the four genotypic classes ac, bc, ad and bd obtained from an ab Âcd cross (Van Ooijen, 2009). For comparison between traits, the effects were calculated as the percentage of the trait mean value. However, for the plasticity QTL the absolute raw effect estimates were retained since the plasticity of traits was based on standardized predictor values. Non-parametric Kruskal-Wallis analyses were performed to verify significant differences between marker genotypic classes close to the peak positions of the QTL. Furthermore, in order to determine phenotypic effects of S. schwerinii alleles, haplotypes were constructed for QTL clusters based on phased data from the linkage map and genotypic data from grandparents. Linkage map and QTL positions were depicted using MapChart (Voorrips, 2002).

Phenotypic and genetic variation
Across the experiments, plant survival was high; in 2014 survival was 99Á9 % in Pustn€ as, 96Á4 % in Cavallermaggiore-IRR and 93Á0 % in Cavallermaggiore-NI. In 2013 the plants grown in Pustn€ as (2-year-old shoots, 6-year-old roots) had generally more shoots with lower shoot diameters compared with the plants grown in Cavallermaggiore (2-year-old shoots and roots) ( Table 2). This also resulted in a higher shoot biomass in Pustn€ as based on the summed basal areas (SumBA13). Independently of treatment, the biomass growth (SumBA13) in Cavallermaggiore was 52 % of the biomass in Pustn€ as. Phenotypic variation was observed for all traits and moderate to high broad-sense predictor heritabilities (H 2 pr ranging between 0Á51 and 0Á86) were detected, indicating favourable conditions for detecting QTL (Table 3). Due to the uneven and patchily distributed drought effects, the Cavallermaggiore-NI experiment exhibited a higher percentage of spatial variation (up to 50 %) than the other two experiments (maximum 28 %), resulting in somewhat lower predictor heritability estimates for the NI experiment (H 2 pr ranged between 0Á51 and 0Á82) compared with the other experiments (H 2 pr ranged between 0Á55 and 0Á86) (Table 3).
Clone predictor correlations across consecutive years were consistently high for all traits (0Á72-0Á94). In Pustn€ as, correlations across multiple years (3-5), including a harvest operation, were also moderate to high (0Á53-0Á85). The clone predictor correlations between Cavallermaggiore-IRR and Cavallermaggiore-NI ranged from 0Á56 (LS) to 0Á84 (BB), indicating moderate to low G Â E interactions (Table 4). For the BB trait, the correlation between Pustn€ as and Cavallermaggiore-IRR was 0Á57 and that between Pustn€ as and Cavallermaggiore-NI was 0Á55. For all other traits, correlations between the Pustn€ as and Cavallermaggiore experiments were considerably lower (0Á18-0Á39) than any correlation estimate across years and irrespective of water availability in Cavallermaggiore (Table 4). These correlation estimates thus indicate considerable crossover-type G Â E interactions between Pustn€ as and the Cavallermaggiore experiments that cannot be entirely attributed to differences in shoot or root age.
The S 1 progeny displayed large variation in phenotypic plasticity for biomass production (based on SumBA13) across the experiments (Fig. 1). Some clones, e.g. 481 and 553, were rather high-producing in all experiments; however, the best-producing clones in Pustn€ as were not high-producing in Cavallermaggiore and vice versa.

QTL are grouped in 18 clusters
We performed QTL mapping analyses in order to identify genomic regions associated with the regulation of phenology and biomass traits (trait QTL) as well as the plasticity of traits (plasticity QTL). Across all traits and experiments, a total of 133 QTL were identified, of which 126 were clustered into 18 groups with overlapping 2-LOD regions (Table 5, Fig. 2). Three of the clusters, located on LG II (cluster C), X (cluster J) and XVII (cluster P), contained QTL that explained most of the phenotypic variation. In these three clusters, 21 QTL out of 53 explained >10 % of the phenotypic variation (Table 5), thus demonstrating the importance of these regions in the regulation of the investigated traits. Some clusters contained QTL for the same trait measured in different experiments and during different years (clusters C, D, E, H, I, M and Q), indicating that the genetic regulation of these traits was fairly insensitive to yearto-year fluctuations and environmental variation. To analyse the stability of QTL across experiments in greater detail, LOD scores for QTL that were significant in at least one experiment were checked for in the other experiments. This resulted in 14 additional QTL with LOD score >3Á0, further supporting the results of common QTL in different experiments (Supplementary Data Table S2). However, some QTL were only found in either Pustn€ as or in Cavallermaggiore (Fig. 2,  clusters A, B, F, G, J, K, L, O, P and R), which suggests a more complex genetic background with different regulation of the corresponding traits in the different environments.
QTL for BB and LS were found in six different clusters (Fig. 2) and clusters C (LG II), D (LG III) and I (LG VIII) harboured QTL for BB and LS assessed in all three experiments, suggesting that these phenology traits are regulated similarly independently of the environment. Furthermore, in cluster C, five of the seven QTL explained >10 % of the phenotypic variation, which demonstrates that cluster C defines a genomic region of particular importance for the phenology traits (Table 5).  Most of the clusters harboured QTL for different biomass traits but only cluster H on LG VI-1 and M on LG XIV contained QTL for MeanD13 and SumBA13 measured in all experiments (Fig. 2), suggesting that some consistency of the genetic regulation of biomass traits across experiments exists even though each QTL explained <10 % of the phenotypic variation. Except for these two clusters, QTL identified for biomass traits in Pustn€ as and in Cavallermaggiore mapped in separate clusters. Especially, in clusters J and P several QTL explained between 10 % and 19 % of the phenotypic variation, indicating that these regions were important in regulating biomass production, particularly in Cavallermaggiore (Table 5). With very few exceptions, QTL identified for the two contrasting water treatments in Cavallermaggiore co-located, demonstrating that water availability alone did not impact the genetic regulation of the traits to any great extent (Fig. 2). The overall results rather suggest that the genetic regulation of biomass traits was much more affected by the climatic differences between Pustn€ as and Cavallermaggiore than the difference in water availability. Moreover, given the relative consistency of the phenology QTL across experiments, it appears that the genetic regulation of phenology was, in general, less affected by environmental differences than were the biomass traits.

Plasticity has a genetic origin
A total of 18 plasticity QTL were found between the contrasting experiments and 11 of them originated from the contrast between Pustn€ as and Cavallermaggiore-IRR (Fig. 2). Plasticity QTL for the biomass traits (Nsh13, MeanD13 and SumBA13) generally co-located with the corresponding trait QTL in Cavallermaggiore (clusters J and P), whereas the plasticity QTL for phenology traits (BB14 and LS13) most often did not. In particular, MeanD13 and SumBA13 plasticity QTL in cluster P for the contrast between the Pustn€ as and Cavallermaggiore-NI explained a considerable part of the phenotypic variation: 16 and 20 %, respectively (Table 5). On the other hand, biomass plasticity QTL very seldom co-located with the corresponding trait QTL detected in Pustn€ as (only two cases in cluster C exist). In general this shows that, for the biomass traits, plasticity and trait variation among the progeny grown in Cavallermaggiore were correlated, whereas this was not the case for plasticity and trait variation in Pustn€ as.
There was only one significant plasticity QTL for the contrast between Cavallermaggiore-IRR and Cavallermaggiore-NI (Table 5, Fig. 2). This MeanD13 plasticity QTL was located in cluster P close to the corresponding plasticity QTL for the other two contrasts (Pustn€ as versus Cavallermaggiore-IRR or Cavallermaggiore-NI), indicating a similarity in the genetic regulation of plasticity independently of contrast.
For the between-year contrasts, eight plasticity QTL were found, of which the majority were for biomass traits in Pustn€ as and for longer time intervals, but none of them explained >6 % of the variation (Supplementary Data Table S3).

Large influence of S. schwerinii alleles on phenotypic expression
In order to trace phenotypic effects of S. schwerinii alleles and haplotypes, we combined results from the interval mapping analysis, the Kruskal-Wallis analysis and phase information from the linkage map, and, if available, we also used the marker raw data from the genotyping of grandparents (Supplementary  Data Table S4). We observed differences in parental allelic effects of significant trait QTL of up to 38 % of the mean value (Table 5, Supplementary Data Table S4) and paternally segregating alleles were associated with differential effects of considerable size more frequently than were maternally segregating alleles. In fact, paternally segregating alleles had differential absolute effects >10 % of the mean values in 54 of the 115 trait QTL, whereas the corresponding number was 22 for maternally segregating alleles. The hybrid origin of the male parent with alleles originating both from S. viminalis and from S. schwerinii is probably the cause of the frequent occurrence of QTL whose underlying alleles show a high differential effect (Table 6, Table S4). In cluster C, which harboured QTL for both phenology and biomass traits, the effect of S. schwerinii alleles was traced by constructing an S. schwerinii haplotype for this region. The S. schwerinii haplotype gave both an earlier start and earlier leaf senescence independent of experiment. The haplotype also had a positive effect on the number of shoots in all experiments together with a decrease in MeanD in the Cavallermaggiore experiments, ultimately leading to an increase in SumBA13 and FW in Pustn€ as only. Paternally segregating alleles associated with large effects were found in cluster P, where specifically SumBA13 and SumBA14 in the Cavallermaggiore experiments showed effects ranging from 26 to 38 % (Table 5, Table S4). These effects are connected to an S. schwerinii allele that reduces the biomass growth of the heterozygote genotype (Table 6, Table S4). Similar haplotype effects on SumBA in Cavallermaggiore, although weaker in TABLE 4. Clone predictor correlations across pairs of contrasting environments (columns) for the assessed traits (rows) with 95 % confidence intervals in brackets. Pu 2-3 , Pustn€ as at root ages 2-3 years (assessed in 2009-2010); Pu 6-7 , Pustn€ as at root ages 6-7 years (assessed in 2013-2014). Cavallermaggiore was assessed at root ages 2-3 years only (2013)(2014). IRR 2-3 , irrigated Cavallermaggiore experiment; NI 2-3 , non-irrigated Cavallermaggiore experiment Trait Pu 2-3 -IRR 2-3 Pu 6-7 -IRR 2-3 Pu 2-3 -NI 2-3 Pu 6-7 -NI 2-3 NI, no information.
magnitude (5-29 %), were associated with QTL in cluster J. In cluster R on LG XIX, QTL for different biomass traits measured at Pustn€ as entailed effect magnitudes of >20 % in comparison with the mean associated with the replacement of one paternally segregating allele with another (Table 5). Also for plasticity traits, the largest effect magnitudes in general originated from paternally segregating alleles (15 cases out of 18; Table 5). The S. schwerinii haplotypes were, for most biomass plasticity QTL (clusters B, J and P), connected to a higher biomass production in Pustn€ as relative to production in Cavallermaggiore (Table 6).

DISCUSSION
The analyses of variation in key phenotypic traits in SRC willows across multiple environments and years led to the detection of extensive variation in biomass traits and substantial phenotypic plasticity among the offspring in the S 1 hybrid population. Clone predictor correlations within traits but across years were moderate to high for all traits, indicating a fairly consistent genetic regulation over time and that shoot and root ages only have a limited impact. The clone predictor correlations were likewise substantial between the two water contrasts in Cavallermaggiore and, in particular, higher compared with the correlations between the climatic contrasts between Pustn€ as and Cavallermaggiore, indicating greater G Â E interactions for the climatic contrast compared with the water contrast. The clone predictor correlations thus suggest that the climatic difference influenced the plants to a greater extent than differences in shoot and root ages or the difference in water availability. We measured phenotypic plasticity for each clone as the difference in standardized trait means across contrasting experiments and between years, and the greater the difference the greater was the phenotypic plasticity. We found different plasticity patterns, as illustrated for the summed basal area in 2013 (SumBA13), where some clones had high values in Cavallermaggiore and low values in Pustn€ as while other clones showed the opposite pattern. The occurrence of plasticity QTL for many traits demonstrates that phenotypic plasticity of phenology and biomass traits across the contrasting environments is genetically determined to a substantial degree. As expected from the low clone predictor correlations, most plasticity QTL were found for the two climatic contrasts compared with the water contrast. Interestingly, plasticity biomass QTL often colocated with trait QTL, demonstrating that they are potentially controlled by the same genetic loci. Furthermore, plasticity biomass QTL mostly co-located with the corresponding biomass trait QTL from the Cavallermaggiore experiments and not from Pustn€ as, indicating that the variation in plasticity in the population overlaps with the trait variation in Cavallermaggiore, but not in Pustn€ as. In contrast, phenology plasticity QTL overlapped only once with the corresponding trait QTL, thus suggesting that the genetic regulation of that plasticity is independent of that of the trait per se.
Willow biomass QTL displayed substantial climatic sensitivity as most QTL were detected in either Pustn€ as or Cavallermaggiore. The difference between the experiments was also evident by the greater biomass production in Pustn€ as   compared with Cavallermaggiore. The greater biomass in the Pustn€ as experiment in 2013 was influenced by the fact that the plant roots in the Pustn€ as experiments were older, and thus better established, than those in Cavallermaggiore when assessed in a given calendar year. Biomass production in Cavallermaggiore was considerably higher in the second year compared with the first year, supporting the idea that establishment as well as root age influence production, which has previously been demonstrated in willows (Nordh, 2005). However, given the considerable G Â E interactions observed between Cavallermaggiore and Pustn€ as and the environmental specificity of the detected QTL, it cannot be excluded that plants of the S 1 population were better adapted to the local climate in Pustn€ as compared with Cavallermaggiore and as a result contributed to a greater accumulated biomass in Pustn€ as. The environmental specificity of the genetic background of biomass traits shows how complex these traits are, which means that a large number of genetic loci and multiple pathways are involved in controlling biomass production, and these might vary depending on the environment. Furthermore, given the rather obvious climatic differences between Cavallermaggiore and Pustn€ as it is not surprising that different genetic pathways are involved in regulating biomass production. On the other hand, stable QTL were found across years within the same environment for biomass traits, indicating concordant genetic regulation over root and shoot ages, which is furthermore supported by the high clone predictor correlations across years. It should, however, be noted that a few plasticity QTL across years were identified, but they explain very little of the total variation. In taxonomically related poplars, stable biomass QTL were also identified across multiple years (Rae et al., 2009;Muchero et al., 2013). In contrast to the biomass QTL, most of the phenology QTL were stable across all experiments and years. Phenology QTL in Pustn€ as co-located with previously identified QTL for the same traits , which furthermore supports the stability of QTL across years. Photoperiod (i.e. daylength) alone or in combination with temperature is the main environmental cue regulating both autumn phenology traits (e.g. leaf senescence) and spring phenology traits (e.g. bud burst) (Howe et al., 1996;Keskitalo et al., 2005;Lagercrantz, 2009;Rohde et al., 2011a). Previous studies have found that these responses to photoperiod and temperature are under strong genetic control (Bradshaw and Stettler, 1995;Keller et al., 2011;Rohde et al., 2011b) and that populations are adapted to the local conditions at the site of growth (Morgenstern, 1996). Local adaptation has been described in several species, and it can be seen as clines in phenology traits (Morgenstern, 1996). Nonetheless, despite phenology QTL consistency, the low and moderate clone predictor correlations between Pustn€ as and Cavallermaggiore for leaf senescence and bud burst, respectively, still indicate an influence originating in climatic differences between the two geographical regions. In agreement with this, Fabbrini et al. (2012) found considerable G Â E interactions for bud set in Populus nigra across contrasting environments and concluded that temperature was one of the major causes.
The complex genetic background, particularly for biomass traits, which display extensive environment specificity, leads to major challenges for breeders. A rough trade-off between performance stability across environments and high production (SumBA13) was observed even though a few clones with relatively high biomass production in both experiments in the climatic contrast exist. Thus, it seems that the maximization of biomass production requires separate selection procedures for each environment. Given the frequent use of a specific S. schwerinii genetic donor cultivar ('79069') in the breeding of willows for biomass, we were interested to know what phenotypic effects haplotypes originating from that S. schwerinii parent actually had on the measured traits. Also, here we found a complex pattern, as phenotypic effects of S. schwerinii haplotypes differed among QTL clusters and experiments. However, S. schwerinii haplotypes derived from QTL always conferred an early bud burst and early leaf senescence independent of experiment. Moreover, S. schwerinii haplotypes were often associated with increased biomass production in Pustn€ as, but mostly with decreased biomass production in Cavallermaggiore. In terms of biomass production, this implies that several genetic key elements of S. schwerinii confer poorer plant adaptation to the climate in northern Italy in comparison with the corresponding S. viminalis elements, but such a pattern was not observed in Sweden. This is further supported by the observation that the S. schwerinii haplotypes for plasticity of biomass traits always conferred a lower biomass production in Cavallermaggiore in relation to that of Pustn€ as. Moreover, the biomass plasticity QTL co-located frequently with the regular trait QTL whose S. schwerinii haplotypes per se decreased biomass production in Cavallermaggiore. Since S. schwerinii has a north-eastern distribution in Asia (Skvortsov, 1999) and does not occur naturally in Sweden or Italy, the climate is probably more similar to Swedish conditions compared with northern Italian conditions. Our results thus suggest that individuals comprising a heavy genetic S. schwerinii component could readily be used when breeding for a more northerly climate but would likely not be the best suited for a southern one. The rationale behind comparing environments with contrasting climate is the desire to develop clonal cultivars that contain the valuable rust resistance factors originating from S. schwerinii but that nonetheless are suitable also for southern European conditions. The stability of some of the QTL indicates the possibility of using these QTL in MAS even though for biomass traits the effect of the S. schwerinii allele appears to vary depending on environment. Given the patterns and characteristics of the QTL detected and the considerable G Â E interactions (low clone predictor correlations) for the climatic contrast, it seems reasonable to separate the breeding material into two different breeding populations aimed at the different environments (Namkoong et al., 1988). In MAS the same allele could be selected for or against depending on the breeding objective and, because we identified QTL that could be used specifically in separate environments for both biomass and phenology traits, we conclude that divergent breeding in divided populations could even be facilitated by MAS.
Further studies are needed to be able to select markers within QTL for use in MAS, but these are the first important findings to identify genomic regions for willow biomass traits across contrasting environments in Europe.

SUPPLEMENTARY DATA
Supplementary data are available online at https://academic. oup.com/aob and consist of the following. Table S1: management and measurements in the field experiments. Table S2: overview of QTL across environments and environmental contrasts (plasticity). Table S3: summary of QTL analyses of across-year plasticity traits for all experiments. Table S4: grandparental genotypes and Kruskal-Wallis test for markers in peak positions of selected clusters.