Common garden comparisons confirm inherited differences in sensitivity to climate change between forest tree species

Sources of data

In the framework of the European Union joint project FORGER (#289119), the response functions (reaction norms) of these four model species were investigated in European provenance tests by the Hungarian and French authors of this study. Part of the dataset from these multiple field provenance tests was used for this analysis. Fourteen populations of Scots pine (Pinus sylvestris) were selected from four trials performed by the VNIILM Institute, Pushkino, Russia (Shutyaev & Giertych, 1998; Shutyaev & Giertych, 2000); 34 populations were chosen from five trials of the international Norway spruce (Picea abies) trial series of IUFRO (Újvári-Jármay, Nagy & Mátyás, 2016); nine populations of European beech (Fagus sylvatica) were selected from six trials of an international IUFRO series (Horváth & Mátyás, 2016) and 14 populations of sessile oak (Quercus petraea) from 23 trials of the “Madsen collection” were used (Madsen, 1990; Sáenz-Romero et al., 2017).

Both the test sites and populations selected for the comparison were highly diverse. Layouts and test ages were not uniform, because the initial experiments had different individual goals and constraints concerning seed procurement and the planting of material due to seed crop fluctuations and differences in site accessibility. Randomized complete block designs were used for all four species, with blocks and replicated plots (for Picea abies, replicated single-tree plots) to account for within-site variation. The selection of populations for analysis was based on their presence in the trials (the populations included in trials were not the same everywhere) and data availability.

The field tests were managed independently by the participating partners, according to their own management guidelines and standard silvicultural practices. Details of the test sites and populations used are provided in the Supplementary Information. Geographic and basic climatic data for the chosen populations are provided, by species, in Table S1; data for the test sites are provided in Table S2. A map of the seed sources and test sites is provided in Fig. 1. Age at tree measurement varied between the available databases. Height was measured at the age of eight years (from planting) for European beech and sessile oak, 12–17 years for Scots pine, and 16–19 years for Norway spruce (Table S3).

The positions of the sampled populations and trial sites in the climatic space are presented in Fig. 2A, and the climate ranges covered by the provenances represented in the field tests and of the test sites relative to the species climate range are shown in Figs. 2B and 2C. The four species have very different geographic distributions, but the sampled ranges overlap at least partially and sufficiently for comparison. The test sites were in average slightly warmer than the sites of origin of the populations tested. At the time at which the tests were established, well before climate change became a major concern, it was usual to establish such experiments in favorable site conditions. The shift in climate conditions already underway at the time may have contributed to this. (Note that the reference climate periods for provenance locations and test sites are different—see below.) The locations of the provenance and field test sites for P. sylvestris cover a much smaller interval than the species as a whole, due to the relatively limited size of the range selected for investigation relative to the immense natural distribution of this species.

For the statistical analysis, we used means per population per block, except for the Picea abies trial network, which was established as a single-tree plot system with 25 replicate, for which means per site and per population were used. For technical reasons, our dataset was not balanced across species, as different numbers of populations and field sites per species were included in the analysis.

Hereafter, we refer to the assemblage of mother trees (forest stands, compartments, from which seeds were collected), represented by their seedlings in the experiments, as “populations”, whereas the origin of the population (the sampling site) is referred to as their “provenance” or “seed source”.

Climate data

Contemporary climate data for each provenance and test site were estimated from the data of Wang & Hamann (2012). Climate variables are described in Table 1. The time periods used as a reference are listed in Table S3. The climate estimates for the seed sources (provenance locations) and for the test sites were estimated for different time spans, as described by Horváth & Mátyás (2016). For the provenance climate, the earliest available climate data were used, to characterize the conditions to which populations had adapted in the past, through a long process of evolution. Climate data for 1961–1990 were used for the normalized climate period immediately before seed collection for Fagus sylvatica and Quercus petraea. These reference periods are the best estimates of the conditions responsible for the semblance of adaptive equilibrium before human-induced climate change began to accentuate (Sáenz-Romero et al., 2015). In the case of Pinus sylvestris and Picea abies, for which seeds were selected earlier than for the broadleaved species, we used the mean climate data estimates for the 1941–1970 period. The climate data estimates for the test sites were calculated as the mean values for the period from planting to measurement, as in Table S3 , to account for the effect of the climate ¹ in which the trees were actually growing.

Statistical analysis

The main challenge in the comparison of climate sensitivity between different species was identifying the climate variable giving the best fit to the response curve.

The climatic niches of different species are obviously not defined by the same climatic variables, as individual climatic variables can have different impacts on growth in different species. A compromise was therefore required to achieve an instructive comparison of response functions between the four forest tree species. For this purpose, we developed a common model, including the same climatic variables, for all four species. This may appear counterintuitive, but it should be borne in mind that the ecological tolerance range of a species is generally characterized with a single climate factor, even though the effect of this factor has different weights at the trailing (xeric) and front (thermal) limits of the distribution. We adopted a two-step approach to select the best climatic variables. We first explored the overall fit to a global model including the species and population components, by screening a large number of climatic variables. Once the most decisive climatic variables had been identified, we then performed a separate analysis for each species, to retrieve the species reaction norms, which we then compared.

Defining the response to transfer between climates: application of a mixed model

We checked the consistency of the results, using tree height as the dependent variable. We then used tree height and annual tree height increment as dependent variables, for modeling of the plant response to climate at species level, analyzing genetic differences between populations nested within species, with the approach described by Leites et al. (2012a) and Leites et al. (2012b). The mixed effects model (SAS/STAT computer software, release 9.1; SAS Institute Inc., Cary, NC, USA) includes both fixed and random effects. The fixed effects in this context correspond to genetic differences between populations due to climatic conditions, and can be split into three major components (Rehfeldt et al., 2003; Leites et al., 2012a; Leites et al., 2012b; Sáenz-Romero et al., 2017): (a) The effect of climate at the site of provenance, term “C” of the model, which accounts for genetic differentiation between populations due to climate at the provenance site, as a result of past climatic selection. (b) The climatic transfer distance, term “D” [D = (climate at test site) − (climate at provenance); linear or quadratic], which accounts for population differentiation due to differences in climate between the provenance and test sites (Mátyás, 1994); this term can be used to characterize the sensitivity or plasticity of populations as a response to climate transfer, based on the shape of the response function (reaction norm) curve. (c) The D∗C term, which accounts for the interaction between climatic transfer distance and the climate of the provenance site; in other words, the response of a population to transfer distance may depend on the climate of origin.

The random effects of the model account for the sources of variation not explained by the fixed effects: of test sites (environmental effects other than climatic effects, such as soil fertility, texture, drainage, silvicultural maintenance regimes, etc.), of populations (genetic effects other than those due to selection by climate, such as genetic drift, gene flow or inbreeding), and of blocks (Leites et al., 2012a; Leites et al., 2012b).

The full model is very similar to that used and described in detail by Sáenz-Romero et al. (2017), except that we have added a species term and nested population within the species term. Our model can be expressed as follows: (1)${\displaystyle Y_{ijkl}={\beta }_0+{\beta }_1{D}_{ijk}+{\beta }_2{D}_{ijk}^2+{\beta }_3{C}_{ik}+{\beta }_4\left(D_{ijk}\ast C_{ik}\right)+S_i+T_j+P_k\left(S_i\right)+e_{ijkl}}$ where Y_ijk is measured tree height (mean per block per population per site), of the kth population, nested in the ith species in the the lth block, nested in the jth test site; β₀ is the intercept; D_ijk is the climatic transfer distance for the kth population nested in the ith species at the jth test site; C_ik is the value of the climate variable at the provenance of the kth population nested in the ith species; D_ijk∗C_ik is the interaction between climatic transfer distance (for the kth population at the jth field test), and the climate variable at the provenance (of the kth population nested in the ith species); β₁–β₄ are regression coefficients; S_i is a random species effect of the ith species; T_j is a random site effect for the jth test. P_k(S_i) is a random population effect of the kth population nested in the ith species, due to factors other than climate; e_ijkl is the error term. D_ijk, C_ik and D_ijk∗C_ik are fixed effects. S_i, T_j and P_k are random effects (Leites et al., 2012a; Leites et al., 2012b).

This model was adopted on the basis of the ecological assumption that the response function (reaction norm) of populations to any change of climate follows a quadratic, symmetric “bell-shaped” function with a maximum at, or near a transfer distance of 0, i.e., close to the climate of the provenance, assuming local adaptation (Wang, O’Neill & Aitken, 2010; Kapeller et al., 2012; Leites et al., 2012a; Leites et al., 2012b; Yang et al., 2015; Sáenz-Romero et al., 2017). The response to climatic transfer includes a quadratic function (D² term), but the term for climate at seed source (C term) is based on a linear function. A large body of experimental results from common gardens of forest trees is available, indicating a linear, clinal variation of populations along climatic or geographic gradients (Mátyás & Yeatman, 1992; Alberto et al., 2013; for a review: Morgenstern, 1996).

A species-specific model was also applied. This model is similar to that described above, but without the species term: (2)${\displaystyle Y_{jkl}={\beta }_0+{\beta }_1{D}_{jk}+{\beta }_2{D}_{jk}^2+{\beta }_3{C}_k+{\beta }_4\left(D_{jk}\ast C_{jk}\right)+T_j+P_k+e_{jkl}.}$

We used a screening procedure to select climate variables for the linear and quadratic functions (Appendix S2). We found a single provenance climate variable appropriately expressing the clinal variation between populations. The inclusion of additional variables reduced the statistical quality of the models. In theory, multivariate climate transfer effects are intuitively appealing, but they increase statistical complexity and may lead to insurmountable constraints.

Fitting of an “overall” (full) model and species-specific models

We first fitted an “overall” transfer function model, with all four species, and populations nested within species (model 1), to identify the best mixed model with common climatic variables (one for provenance, one for climatic transfer distance). The best variable combination for the response functions (reaction norms) was determined simultaneously for the four species. We compared 49 models, and the one with the lowest AIC (Akaike Information Criterion, Akaike, 1973) value was selected as the best model. Model selection, including preselection of the seven best models for provenance and the seven best models for climatic transfer distance, yielding a matrix of 49 competing models, is described in Appendix S2.

We then refitted the model for each species separately, with the same selected variables, but using only observed data for the species concerned as input, to generate appropriate parameter values for the fixed-effect terms, and individual response functions for each species. These response functions were illustrated by plotting a curve for each species, using the model parameters estimated for the fixed terms for each species, and plotting the mean value of the selected climate variable for the species at the seed source against the best climatic transfer distance variable. The range of climatic transfer distance values (x axis) was extended beyond the measured values to make it possible to visualize the response function trend in full, as most of the field trials were not located at sites with extreme climates.

We checked that the differences in response curve between species persisted even if a different set of climate variables was used, by fitting response curves independently for each species. In this validation process, we selected the best climatic variable for each species separately with model 2. We thus obtained four pairs of climatic variables (for the D and C terms for each species). We then used each pair of variables obtained for a given species on the other three species, implementing model 2. This procedure yielded 16 models in total.

For further improvement of the comparison between species, we repeated the process for full model fitting with the selected climatic variables from models 1 and 2 above, this time using scaled growth data as an independent variable, to lessen differences in age and in species growth potential. We applied the following formula for each subset of data (for each age nested within each species): ${\displaystyle \mathrm{scaled height}=\left[\mathrm{height}-\left(\mathrm{minimum height}\right)\right]∕\left(\mathrm{maximum height}-\mathrm{minimum height}\right).}$

We then used the estimated fixed-term parameters to determine scaled-height response functions for each species, using the common climatic variables. In this analysis, “tree height” is the mean tree height per block.

We repeated the three-step model selection process for annual tree height increment (ATHI = tree height/age) and plotted the predicted values for each species against climatic transfer distance, to confirm the previous results and to eliminate age differences between species.

Exploring genetic clines

We explored the amount of genetic differentiation between populations, probably linked to the shape of the response function, by investigating the expression of clines along environmental gradients, by screening relevant provenance climate variables. We estimated Spearman’s correlation coefficients for the relationship between the mean ATHI per population per site and the provenance climate variables. We analyzed only “good sites” (sites at which the mean per population per site was above the overall mean), as differences in growth potential between populations are best expressed in mild environmental conditions. We ranked the absolute values of Spearman’s correlation coefficient by species. We then selected the three climate variables with the highest Spearman correlation coefficient values for each species. The site effect was eliminated by estimating the best linear unbiased parameter (BLUP) value for each site, and subtracting it from the provenance mean per site. BLUPs were estimated with a reduced mixed model, in which climate at provenance was a fixed effect and site was a random effect. The next step was to standardize the described corrected mean per population per site, by dividing it by the overall mean (across provenances and sites). In other words, we obtained standardized corrected means for each population and site, with a distribution centered on zero. Finally, we performed a linear regression analysis of the standardized corrected means per population per site against the provenance climate values of the three most relevant climatic variables per species identified in previous analyses, to compare the fit and slope of the regression lines. We plotted the regressions, to obtain a visual comparison of the different cline patterns between species.

Fitting a model for unilateral climatic transfer, towards warming

The global climate is changing in a single direction, with a marked shift towards warming. When modeling the future responses of species and populations to climate change, it is important to remember that only efforts to mimic such unilateral change, based on the transfer of populations to warmer and drier conditions, are of practical relevance.

We performed a regression analysis for unilateral climate transfer, with a shift towards warmer/drier climatic conditions, using the same climatic variables in the chosen model. In other words, we regressed the means per population per site against positive climatic transfer distance values only. In biological terms, a second-degree regression model would be plausible, but we nevertheless fitted a simple linear regression model, because the data were highly scattered and we wished to avoid artifactual overfitting. The linear model used for each species was: (3)${\displaystyle Y_{jkl}={\beta }_0+{\beta }_1{D}_{jk}+e_{jkl}}$ where D_jk is the climatic transfer distance for the kth population, of a given species, at the the jth test site.

Fitting of the “overall” (full) model: determining the best variables for the climate transfer function

After preselecting the climatic variables and identifying the overall (full) model best explaining the variability (for all four species), we applied a transfer distance based on annual dryness index (ADI) as the D term of the mixed model and used winter mean maximum temperature (°C) of the climate at provenance as the C term of the mixed model. The statistics of the model with the selected variables are presented in Table 2. Details of model selection and a comparison of fits with models based on other climatic variables can be found in Appendix S2 and Table S4.

An examination of the virtual response surface emerging from the results showed that linear and quadratic values of term D controlled the response of populations (i.e., height growth was determined principally by the effect of climate transfer), with term C, winter mean maximum temperature at provenance playing only a minor, secondary role.

Specific values for the two common variables selected (corresponding to terms C and D) were used for each species as well, to build a model with estimated parameters of the fixed-effect terms for each species. The quadratic response functions are shown for the species in Fig. 3. Curve height is related to the age of the populations at the time of measurement (see Table S3). The curves for the two conifers are therefore higher than those for the broadleaved trees because the trials were older. Thus, when comparing the curves, width and steepness should be considered, rather than maxima.

The model shows that growth decreases when populations are transferred either to warmer/drier sites (positive transfer distances), or to cooler/wetter sites (negative values). The shape of the response curves in Fig. 3 is of particular importance: flat and wide for Fagus sylvatica and Quercus petraea, but with more pronounced slopes for Pinus sylvestris and Picea abies. This means that the (adaptive) responses of populations to change are much more marked for the two conifer species. This results from the ADI² quadratic fixed term, which determines the shape of the response curves and highlights the differences between the modeled conifers Pinus sylvestris and Picea abies, and the modeled broadleaved species Fagus sylvatica and Quercus petraea.

Species transfer response

The trend shown in Fig. 3 was confirmed by the fitting of separate response curves for each species, with the three-step model selection procedure repeated for each species. The shapes of the curves obtained were generally similar to those in Fig. 3, except that the quadratic term yielded a positive value in some cases, resulting in an inverted, biologically unacceptable curve, which was therefore discarded. Curves of this type fitted for Fagus sylvatica (Fig. S1A) and for Pinus sylvestris (Fig. S1B) are presented in the supplementary information.

The response functions (reaction norms) did not have maxima close to a transfer distance of zero (i.e., corresponding to climates identical to the climate of provenance) for any species other than Fagus sylvatica, for which the maximum was very close to zero (on the mean curve for this species, the maximum tree height of 196.6 cm corresponds to an annual aridity index transfer distance value of 0.01; Fig. 3). For Picea abies and Quercus petraea, the maximum height was slightly displaced towards cooler and/or moister sites (for the mean curve of these two species, maximum tree height was predicted to occur at an annual aridity index transfer distance value of −0.02 and −0.04, respectively; Fig. 3). A similar phenomenon was observed for Pinus sylvestris (maximum tree height at an annual aridity index transfer distance value of −0.045; Fig. 3). These results are consistent with the species-specific results of analyses of field tests (e.g., Mátyás, Nagy & Újvári-Jármay, 2010; Rehfeldt et al., 2017).

Response functions may also be displayed as a three-dimensional surface, with the two climatic variables determining measured height (Fig. 4). This illustrates more clearly the difference between the responses of Q. petraea and P. abies. Q. petraea has a flat response, encompassing a larger climatic interval under the response surface (a range of 0.41 for annual aridity index transfer distance, from about −0.23 to 0.18), whereas Picea abies has a much steeper and narrower surface, encompassing a smaller climatic space interval (a range of 0.20 for annual aridity index transfer distance, from about −0.11 to 0.09; Fig. 4).

As in Fig. 3, the predictions shown in Fig. 4 are based on the climatic variables selected for simultaneous fitting to all four species, with subsequent refitting of the model for each species, for the estimation of regression coefficients for each species, with scaled heights (only two of the four species are displayed for the sake of clarity). The quadratic D term and the linear term C determine the shape of the response curve. C covers the climatic gradient along which population samples were collected, whereas D describes the response to climatic transfer distance.

We also modeled the annual tree height increment (ATHI) response (by three-step model selection, Appendix S2, as for tree height) against annual dryness index, which was identified as the best variable to represent climatic transfer distance and temperature difference (TD = MTWM − MTCM = continentality), which was selected as the best variable to represent seed source climate. The results displayed in the supplementary information (Fig. S2; Table S6) reveal a trend similar to that for tree height (Fig. 3), confirming the general pattern of flatter response curves for Q. petraea and F. sylvatica and more pronounced, narrow curves for the two conifer species (Fig. S2).

Genetic clines

The separate exploration of clines for annual height increment within in each species revealed striking differences between species. Linear genetic clines were much more pronounced in Pinus sylvestris, along gradients of spring temperature or precipitation as snow (Figs. 5A and 5B), accounting for 43% to 64% of the total variance (Table 3). Picea abies displayed less pronounced but nevertheless highly significant clines (see Figs. 5C and 5D; 11% to 13% of the variance explained, P <0.007; Table 3) along gradients of winter or annual precipitation. By contrast, the clines for Q. petraea were weaker (Figs. 5D and 5E), explaining less than 8% of the variance (Table 3), and those for Fagus sylvatica (Figs. 5F and 5G) were weaker still, accounting for less than 6% of the variance (Table 3).

Fitting a model for unilateral climatic transfer, with a shift towards warming

The linear regression of annual aridity index transfer distance (independent variable) against mean scaled tree height per population (response variable) clearly showed, for Q. petraea, that growth decreases following transfer to warmer/drier sites (Fig. 6). The regression was highly significant for Q. petraea (P <0.001), and non-significant for P. abies and Fagus sylvatica, with Pinus sylvestris having an intermediate P value (0.070), with growth tending to decrease with increasing climatic transfer distance towards warmer/drier sites. In this analysis, response functions were not extrapolated beyond the measured points, and data for negative transfer distance values were not considered. The high degree of scatter is largely due to the high within-species variance of individual populations.

Contrasting response curves and genetic clines

Response curves based on climatic variables identified in provenance tests have already been investigated separately for different species (e.g., Rehfeldt et al., 2002; Mátyás & Nagy, 2005 for P. sylvestris; Horváth & Mátyás, 2016 for Fagus sylvatica; Sáenz-Romero et al., 2017 for Q. petraea), but we used a novel comparative approach in this study. Interestingly, our results are consistent with those of the previous studies, in which drought-indicating variables, not unlike annual dryness index (e.g., Ellenberg drought index, Horváth & Mátyás, 2016; summer drought period length, Frank et al., 2017a) were found to be significant.

The flatter climatic response curves of Fagus sylvatica and Quercus petraea (resulting in modest growth over a wider range of aridity indices than for Pinus sylvestris and Picea abies) indicates that the populations of these species have higher levels of plasticity, enabling them to adapt to a broader range of environments. The response curves of these species were also flat towards the cooler and more humid environments. By contrast, Pinus sylvestris and Picea abies have steeper response curves with more pronounced slopes, encompassing a narrower range of environmental conditions to which they appear to be adapted. These response curves suggest that spruce and pine are more vulnerable to climate change, as the pronounced cline slopes indicate a higher degree of specialization. In a context of assisted migration, the conifer response curves suggest a large advantage of re-aligning genotypes to the environment to which they are adapted. Transferring seed sources to sites at which the future climate is predicted to be similar to their current “home” climate is being treated as a possible option in operational forestry.

This conclusion is supported by European forest mortality observations. Drought-triggered mortality is observed regionally in oak and beech (e.g., Lakatos & Molnár, 2009), but the destructive decline of Norway spruce due to climate shifts is continent-wide and cannot be explained solely by the known artificial extension of its distribution range by humans (Seidl, Schelhaas & Lexer, 2011; Marini et al., 2017). A Europe-wide investigation of genetic conservation units for the most important forest tree species (Schueler et al., 2014) also revealed that the two conifer species studied here are exposed to significantly stronger climatic threats than oak and beech. Several studies have reported high sensitivity to drought in P. abies (Bošela et al., 2014; Bouriaud & Popa, 2009; Desplanque, Rolland & Schweingruber, 1999; Huang et al., 2017; Lebourgeois, Rathgeber & Ulrich, 2010; Lévesque et al., 2014; Zang et al., 2014).

The lower adaptability of Pinus sylvestris and Picea abies to climatic transfer distance may result from their pronounced differentiation into genetic clines along climatic gradients, due to the divergent selection that has historically occurred along these gradients. In other words, steeper genetic clines probably constrain the response of pine and spruce populations to narrower reaction norms. By contrast, Fagus sylvatica and Quercus petraea have milder clines, which contribute to flatter reaction norms. We suspect that the widths of the climatic gradients to which populations were exposed were relatively similar for the four species, as they occupy similar geographic ranges, at least in Central Europe. Differences in cline steepness between species are therefore likely to reflect intrinsic genetic differences between species.

Mechanisms other than adaptation to climate-imposed selection may also play an important role in the greater phenotypic plasticity of Fagus sylvatica and Quercus petraea. For example, these species may have more efficient mechanisms for dealing with biotic interactions, such as competition, antagonism, mutualism, defenses against insect and pest attacks, and a greater capacity to recover from such attacks. This hypothesis appears plausible in light of the recent decimation of natural conifer populations by beetle infestations linked to ongoing climate change (Allen et al., 2010; Alfaro et al., 2014).

It is tempting to conclude that the more plastic responses of Fagus sylvatica and Quercus petraea might be due to the shedding of their leaves in the fall, or in periods of severe drought, enabling them to escape the worst climatic stresses, but this is probably an oversimplification. Evolutionary genetic and ecophysiological background, and, thus, adaptive strategy, appear to differ between species, making simple generalizations difficult. Even within closely related species, such as those of the genus Abies, the between-species differences in plasticity measured in common gardens are considerable (Frydl et al., 2018). The plasticity of adaptive specialists, with a high between-population variance for adaptive traits, may result in higher levels of sensitivity to climatic conditions, as in Norway spruce. By contrast, in beech, which had a flat response curve to climate changes in our analysis, high levels of within-population genetic variation were detected in comparison tests on seedlings (Frank et al., 2017a; Frank et al., 2017b).

Suboptimal presence of the four species

Our results indicate that the studied species (except for the adaptive generalist species beech, for which results were inconclusive) are preferentially present in environments that are drier or warmer than their optimum for growth. “Non-optimality” is a specific phenomenon of deviation from the expected symmetric shape of response curves. Asymmetric responses have been detected in common garden datasets (Mátyás & Yeatman, 1992; Rehfeldt et al., 2017). Early observations (Namkoong, 1969; Mátyás, 1990) indicated that response functions (reaction norms) constructed from common garden growth data displayed a systematic divergence of maxima from the mimicked position of climate similar to that of the provenance. This lag effect probably indicates a balance in the processes leading to local adaptation between directed selection and other random genetic effects, such as gene flow (see Savolainen, Pyhäjärvi & Knürr, 2007; Kremer et al., 2012). In conjunction with phenotypic plasticity, it may prevent “perfect” optimization, which would, in itself, lead to lower levels of survival in fluctuating environments (Mátyás, Nagy & Újvári-Jármay, 2010). Other factors (such as interspecific competition or other biotic interactions) may also have a strong effect on the contemporary distributions of natural populations.

We further extended “non optimality” towards drier and warmer sites, by analyzing responses to unilateral transfer. We thus focused on the most important practical question in forestry today: how the growth of populations is affected by a simulated shift in climate towards warmer/drier conditions. The simple linear regressions applied are free from extrapolation bias and confirmed the rapid loss of productivity under the anticipated changes in climatic conditions, but this trend has already been reported in species-specific analyses of common garden tests (Sáenz-Romero et al., 2017; Mátyás, Nagy & Újvári-Jármay, 2010; Újvári-Jármay, Nagy & Mátyás, 2016; Horváth & Mátyás, 2016). These results highlight the need to develop better models for the lower (xeric) limits, which has become particularly urgent in light of predicted climate change (Mátyás, 2010).

Limitations of this study

We acknowledge that our conclusions regarding the contrasting responses of species may be driven by constraints inherent to the comparative approach used. The selection of a single “universal” climate variable, such as ADI (annual dryness index) seems, intuitively, to hinder the combined analysis of four species with different ecological niches. This approach undoubtedly increased the proportion of the variance remaining unexplained and lowered the significance of the results. However, given the continent-wide shift in climatic conditions at all sites, all tree species are simultaneously confronted with a real threat to their vitality and a risk of mortality across the distribution area: extreme droughts (Allen et al., 2010; Allen, Breshears & McDowell, 2015). ADI transfer difference (ADI_d) expresses the balance between temperature (degree days above 5 °C) and water availability, and, indirectly, the probability, intensity and frequency of extreme events regardless of geographic location. Through the use of differential values (transfer distances), we were able to take the relative differentiation due to climatic change across the geographic range into account.

Another constraint inherent to the comparative approach relates to that the choice of model. We assumed a quadratic response for each species, based on the generally accepted ecological assumption that the response of populations to changes in environmental factors follows a symmetric “bell-shaped” curve (e.g., Griffiths et al., 2000). The quadratic response model has proved successful for explaining variation in boreal conifers, but appears to be less effective for other species, including broadleaved species in particular, which display greater phenotypic plasticity in their response to site attributes, including climate.

The calculated parameter statistics for the mixed model confirm that the species-specific growth responses of populations are confounded with tree age, and partially with local site conditions. The largest variance component was the species effect (63.0% of random effects, see Table 1), but this value explains above all differences in age and site potential. The contribution of the exclusively genetic between-population component appears to be modest (3.3%, Table 1). It should be stressed that this figure accounts for only a part of the genetic variation captured by the analysis. Three fixed components (D, D² and C; related to climate as a selective force), and the random component (resulting from effects other than climate) also contribute to genetic variation within species, and the transfer distance explains provenance-by-environment interactions.

The unbalanced, heterogeneous experimental setting may also have affected the outcome of this study. For researchers not familiar with field data from continent-wide common garden test series, the low level of model fitting across species may appear too uncertain for any practical conclusions to be drawn. A number of uncontrollable factors may have biased the response of individual species relative to others. In addition to management practices, important interacting biotic factors (e.g., differential pest pressures), and variations of local site conditions (fertility and water-holding capacity of the soil, hydrology) may be confounded with climatic effects. The consideration of site potential, as understood in forestry (an inseparable complex of soil, water and climate factors) might have accounted for much of the remaining variation. Insufficient detail in the description of experimental sites contributed to the high level of unexplained variation. In an international experimental series, in which data retrieval is dependent on different providers in different countries, it can be very challenging to obtain meaningful site data (exposure, soil type, rainfall, microclimate etc.) because of missing data for standard descriptors. Even the assessment of the most important abiotic factor in this analysis, climate, is dependent on the use of statistical concepts to interpolate digital surfaces from data points corresponding to meteorological observations. This situation, which is typical for large-scale common garden tests, highlights the need to describe site characteristics more accurately, to increase the predictive power of these valuable and often unique experiments.

Finally, the climatic range encompassed by the sampled populations and test sites for Pinus sylvestris and by the test sites for Picea abies (particularly for rainfall for Picea abies), being smaller than covered by the natural distributions of these species, may have limited our capacity to model the response function for the species as a whole with sufficient robustness. However, even with this limitation, we were able to capture the pattern of genetic differentiation between populations, with the detection of very pronounced clines for P. sylvestris, as shown in Fig. 5.

Despite these inherent limitations, our results were consistent across the various empirical validation steps applied, such as the use of different environmental variables, or raw vs. scaled data.