We assessed population response of Norway spruce to climate change by a Europe-wide transplant experiment, the International Provenance Test 1964/68 with Norway spruce (IPTNS) initiated by the international Union of Forest Research Organizations (IUFRO). A total of 1,100 stands was sampled across Europe (Figure 1) and sown in spring 1964 (complete list of provenances may be found in Ujvári-Jármay et al., 2016). After uniform pre-treatment at one single nursery in Germany the material was distributed to be planted at 20 experiments established in 13 countries. Due to the large number of populations, the complete set was split into 11 blocks of 100 populations each, which were assembled, based on a stratified sampling to contain material from the entire sampled area, approximately representing a comparable climate space. These blocks can be regarded as a complete trial by themselves and were assumed to be equal in mean and within-block variance, thus easily comparable (Krutzsch, 1992). For seven experiments all 11 blocks were established at one site (e.g., EXP 7, 8 in Sweden), while others were established at two (e.g., EXP 5 in Norway and EXP 16 in Scotland) or even up to a maximum of nine sites for EXP 15 in Germany. Trials were designed in a similar way with single tree plots planted in a randomized complete block design with 20 or 25 replicates, at a spacing of 2 × 2 m. One exception is EXP 16 in Scotland, where six trees are planted per plot in four replications. All trials were managed independently by the various partner institutions, according to their own guidelines and silvicultural practices.

Measurements were taken at recurring intervals, but their timing varied among institutions for the same reason. Here, we concentrate at tree height as the most commonly assessed proxy for fitness, in which multiple ecophysiological traits affected by the environment are cumulated, e.g., taller trees compete more efficiently for light, water, and nutrients and are likely to have higher reproductive success at maturity (Aitken and Bemmels, 2016). We use juvenile tree height of 11 years available from 10 experiments, with different numbers of provenance blocks being distributed to a total of 32 site locations (97 blocks of 100 provenances) (Figure 1). Supplementary Table 1 provides further detail about site locations and reports the corresponding references to the used provenance mean data. Heights were either directly measured or interpolated to the common measurement age based on two concurrent measurements of a maximum age difference of 5 years (age 9–14 for EXP 20) using the following equation:

where HT = height, i and j = age of measurement, with i measured before j. The distribution of HT₁₁ used in the analysis is displayed by trial site location in Supplementary Figure 1.

Current occurrence of Norway spruce across Europe was obtained from EU-Forest by Mauri et al. (2017), the most comprehensive Dataset on European tree species occurrence up to date, which harmonizes national forest inventories from 21 European countries and aggregates them to an INSPIRE-compliant 1 × 1 km grid. All records of Norway spruce were treated as true presences (70,310 occurrences), while all other species records were summarized to a grid cell basis from which an equivalent number of pseudo-absences was selected (prevalence = 0.5). Similar to Iturbide et al. (2015), we first set an exclusion buffer of 10 km around the presence points, to then randomly select the pseudo-absences among the remaining target group absences. Barbet-Massin et al. (2012) found this to result in the highest predictive accuracy when using classification techniques.

Biologically relevant annual climate data was generated with the ClimateEU v4.63 software, publicly available at http://tinyurl.com/ClimateEU, based on methodology described by Marchi et al. (2020) (Table 1).

Climate parameters of provenance (suffix “p”) and test site (suffix “s”) for the calibration of the TraitSDM were extracted for different time intervals as described by Sáenz-Romero et al. (2019). The normal period of 1941–1970 was used for provenance climate characterizing the conditions to which the sourced provenances had adapted in the past. To characterize the climate of the test sites, mean values of climate estimates were calculated for the period 1968–1975. The OccurrenceSDM was calibrated with the widely used reference period 1961–1990 that largely predates a significant anthropogenic warming signal. For each of the three sets of predictor variables we performed a principal component analysis to inspect patterns of multicollinearity (Supplementary Figure 2).

Provenances were further characterized by geographic parameters to account for other factors affecting among population variation such as gene flow, genetic drift, or postglacial colonization in the TraitSDM that may result in spatial patterns independent of climate (Wang et al., 2010). These included latitude (LAT), longitude (LONG), and elevation (ELEV) and their transformations and interactions (i.e., LAT², LONG², ELEV², LAT × LONG, LAT × ELEV, LONG × ELEV, (LAT × LONG)², (LAT × ELEV)², (LONG × ELEV)², LAT² × LONG, LAT² × ELEV, LONG² × ELEV, LAT × LONG², LAT × ELEV², and LONG × ELEV²).

For spatial habitat modeling of both the OccurrenceSDM and TraitSDM, we used 1 × 1 km resolution climate grids in Albers Equal Area projection. To make predictions for the 1961–1990 normal period, and for three projected time slices, the 2020s (2011–2040), the 2050s (2041–2070), and the 2080s (2071–2100) ensemble projections provided by ClimateEU were used. These ensembles average across 15 Atmosphere-Ocean General circulation models corresponding to the IPCC Assessment Report 5 (CanESM2, ACCESS1.0, IPSL-CM5A-MR, MIROC5, MPI-ESM-LR, CCSM4, HadGEM2-ES, CNRM-CM5, CSIRO Mk 3.6, GFDL-CM3, INM-CM4, MRI-CGCM3, MIROC-ESM, CESM1-CAM5, and GISS-E2R), representing all clusters of similar AOGCMs, and having high validation statistics in their CMIP3 equivalents (Marchi et al., 2020). Each time period is predicted for a Representative Concentration Pathway (RCP) with a radiative forcing value of +8.5 W/m². This forcing value results in a projected global warming of approximately +2.0°C (±0.6) by the 2050s and +3.7°C (±0.9) by the 2080s (Marchi et al., 2020).

We followed the ODMAP protocol proposed by Zurell et al. (2020) to build the OccurrenceSDM that forecasts species occurrence probability based on observed species presence from the described national forest inventory data. First, those climate predictors of Table 1, which explain most of the variation in observed presence and absence were identified by assessing univariate variable importance in terms of Akaike information criterion (AIC). If these variables showed issues of multicollinearity, as visible by their grouping in the PCA (Supplementary Figure 2), the less important variable was removed from each pair of variables having a correlation |r| > 0.8 (Dormann et al., 2013).

The OccurrenceSDM was fitted using an ensemble approach accounting for different assumptions and extrapolations to new environments by different algorithms (Araújo et al., 2019). We implemented five popular algorithms, namely generalized linear models (GLM), classification tree analysis (CTA), generalized boosted regression models (GBM), artificial neural networks (ANN), and random forest for classification and regression (RF) using the multi-method platform biomod2 (Thuiller et al., 2020) implemented in R statistical environment (R Core Team, 2019). Default settings that are designed to balance model complexity and overfitting were used for all modeling algorithms, despite for GLM, where we specifically allowed pairwise interactions. Model output were occurrence probabilities, which were then ensembled to a consensus model by combining the median probability over the selected models. The median was chosen as it is known to be less sensitive to outliers than the mean. For each model as well as the final ensemble, true skill statistics (TSS), area under the receiver operating characteristic (ROC), model sensitivity (ability of the model to predict true presence), and model specificity (ability of the model to predict true absence) were recorded. To account for model uncertainty, the data was partitioned to 75% of the data to train the model while 25% was retained for model testing. This was done for three independent runs, each with a randomly selected set of training and testing data. For each model run model evaluation statistics were recorded and finally averaged for reporting.

The variable importance of the climatic predictors was assessed by an internal function of biomod2, which relies on the comparison of model predictions derived from the original dataset with predictions from permutated datasets (Thuiller et al., 2012, 2020). Permutated datasets are created by randomizing one environmental variable at the time, while keeping all others as original to then create a new prediction. Predicted occurrence probabilities resulting from the original dataset and the permutated dataset are then compared by calculating a Pearson correlation coefficient. A small correlation coefficient indicates a high importance, i.e., the predictor variable strongly affects the prediction, while a high coefficient indicates minor importance. This procedure was repeated three times for each model run and mean correlation values were derived across the three independent model runs. The finally reported variable importance corresponds to 1 minus the average correlation coefficient, so that high values reveal a high importance and low values a low importance. This procedure was repeated for each predictor variable and model algorithm.

To identify the climatic conditions of trial site location and provenance origin that are the main drivers of growth performance across the species range, we used a URF according to Wang et al. (2010). This function basically integrates response functions, being reaction norms of individual provenances to environmental variation, and individual transfer functions that express among population variation as genetic clines into a single model. In addition, other factors such as gene flow, demographic history and colonization events can have a considerable effect on among-population variation resulting in spatial patterns independent of climate. As stated by Wang et al. (2010), “these spatial genetic effects may be explained to some extent by geographical proximity of provenance, assuming population proximity is correlated with relatedness, and shared evolutionary history,” and were further included to the function. The mean tree height HT of population j at site i was fitted according to Equation 2:

where β represents the intercept and regression coefficients, S_i and P_j are one or more climate variables for site i and provenance j, respectively, S_i×P_j is the interaction between S_i and P_j, G_jis one or more geographic variables of provenance j, and e_ij is the residual error.

Multiple regression models were constructed using different variable combinations with a maximum of two climate variables for site and provenance and one additional geographic variable for provenance. We removed models showing high collinearity between climatic variables (Pearson coefficient r ≥ 0.70) and approved that the quadratic terms were of convex form (decreasing tree height towards both ends of climatic gradients), which is biological plausible and allows to estimate provenance climatic optima (β₂ < 0, β₄ < 0, or β₇ < 0). Significance levels for the predictors to stay in the model were set to α ≤ 0.001. Finally, we selected the best model based on the highest proportion of variance explained by partial R² and lowest AIC.

In order to validate the derived model, it was used to predict height growth performance for a set of independent data. These were 15 additional trial site locations (47 blocks of 100 provenances) of the IPTNS 1964/68 series, corresponding to seven experiments from six countries (Austria, Belgium, Czech Republic, France, Norway, and Sweden) that were not included in model calibration due to a different age of the available data (age 9, 10, 12, or 14 depending on trial site, Supplementary Table 1). Predicted heights were compared to observed tree heights using Spearman’s rank correlation as proposed by Chakraborty et al. (2016).

Spatially explicit predictions at the continental scale are derived according to the procedure described by Wang et al. (2010). Predictions are made once using “local” provenance, i.e., provenances originate from the climate conditions that are the same as at the planting site, as well as “optimal” provenance, where the, according to the provenance term in the model, best performing provenance is assumed to be planted. Growth of local provenance under current climate is predicted by setting values for site climate in the URF equal to provenance climate (TD_p = TD_s). To estimate growth under the three future scenarios (2020s, 2050s, and 2080s) current local climate is used for provenance and future climate for site. To identify the best performing provenances under reference and future climate conditions, the first order partial derivatives of the URF were drawn for climatic and geographic variables of provenance. These partial derivatives were set to zero and then solved to obtain the values for the chosen variables of the provenance that yields greatest growth at the site, i.e., at any given grid cell of the study area (Wang et al., 2010). As a log-transformation was applied on the response variable height for model fitting, predictions were back transformed prior to spatial mapping.

To translate occurrence probability from the OccurrenceSDM to binary presence/absence we identified the maximum sum of sensitivity and specificity as optimal threshold (Liu et al., 2005).

The TraitSDM predicts growth performance, which may not be directly compared to binary species occurrence and therefore requires conversion. Predictions were translated to binary presence/absence applying a threshold surface accounting for decreasing growth potential towards the margins (high elevations, northern, and southern edges) of the distribution. For example, average performance on trial sites in western Germany is over 2 m at age 11, while only about 0.5 m at the northern most site in Lappkojberget, which needs to be accounted for when translating predictions of height growth potential into species presence and absence. Therefore, we applied a productivity-oriented threshold according to the following procedure: First, all provenance means were standardized and a quadratic regression was fitted based on geographic latitude and elevation. Elevation was only used as linear term, as a negative linear cline with increasing elevation is commonly observed. The derived regression was then applied spatially, giving each pixel across the landscape a specific threshold value according to its latitude and elevation. Height growth predictions by TraitSDM were likewise standardized and compared to the threshold value to segregate presence (higher than threshold) and absence (lower than threshold).

Finally, the spatial extent of species presence predicted by different models (OccurrenceSDM and TraitSDM) and projected time slices (2020s, 2050s, and 2080s) was compared to a common baseline: the spatial extent of the binary presence for the reference period 1961–1990 as modeled by the OccurrenceSDM. The change in species range was quantified with the biomod2 built-in function BIOMOD_RangeSize().

All statistical analysis was conducted in the R programming environment for statistical computing and visualization (R Core Team, 2019).

After the inspection of multicollinearity in the predictor variables of the presence/absence data, four weakly correlated climate variables remained for fitting the OccurrenceSDM, namely climate moisture deficit (CMD), mean coldest month temperature (MCMT), mean driest month precipitation (MDMP), and the temperature difference between mean warmest and coldest month (TD) as indicator of continentality. MCMT and CMD contribute strongly to correctly predicting presence and absence in the forest inventory data, with GLM, CTA, and GBM depicting CMD as most relevant followed by MCMT, while it was the opposite for the algorithms ANN and RF (Table 2). Partial response curves are provided in Supplementary Figure 3 illustrating the sensibility of each algorithm to the climatic predictors. The predictive performance of the models was excellent with both training and testing data, with mean TSS ranging from 0.87 to 0.95 and mean ROC ranging from 0.97 to 0.99 (Table 3).

According to the OccurrenceSDM, the range of Norway spruce decreases for all predicted climate periods (2020s, 2050s, and 2080s) as compared to the reference period, first back towards the extents of the natural distribution and thereafter even further (Figure 2). Habitat suitability declines particularly in Southern, Central, and Eastern Europe, whereas it remains stable in the in the North (>60°N) and at higher elevations (e.g., in the Alps).

From the 14 climate variables tested for their performance in predicting height growth, two variables of the trial site, one climatic and one geographic variable for the provenance, as well as one interaction between site and provenance climate were found to be significant predictors (Table 4). The iterative process of comparative model selection is illustrated in Supplementary Table 2 (step 1–7). The final model explained 75.8% of the variance in the data and included the length of the frost-free period (FFP_s), which can be considered an indicator for the length of the growing season, as well as the summer-heat-moisture index (SHM_s). For provenance, it included the temperature difference between warmest and coldest month as indicator of continentality (TD_p) and the geographic latitude (LAT_p). As interaction term FFP_s × TD_p stayed in the model (the second interaction term tested SHM_s × TD_p was rejected as it did not fulfill α ≤ 0.001). Site climatic variables explain large amounts of the variation in tree growth with 82.7%, while 16.8% is explained by geography and climate of the provenance. Herein, provenance climate could only capture a limited proportion of 2.2%, whereas provenance latitude and its squared term shape a great proportion of the observed variation with 14.5% (Table 4). Despite being significant the interaction term between site and provenance climate accounts for less than 1%.

The model was evaluated with independent data from 15 locations (47 blocks of 100 provenances), well distributed across the study area (Supplementary Figure 4). Spearman’s rank correlations between predicted and observed provenance height had positive coefficients ranging between 0.06 and 0.70, with r > 0.5 for 29 blocks of 100 provenances. The full model (equal to Table 4) performs moderate to well for sites located within the southern center of the natural distribution, those close to the western border of spruce’s artificial distribution in northern France and Belgium, as well as for the southern edge of Norway. However, model performance is poor for Hjuleberg in Southern Sweden (r < 0.4 for six of 11 provenance blocks). The combination of geographic and climatic variables for provenance in the model improved the predictive power of the model, as illustrated by the comparison with two reduced models that exclude either provenance climate or provenance geography from the full model (Supplementary Figure 4).

Predictions of juvenile height based on the TraitSDM are illustrated in Figure 3. All areas with a positive height growth potential basically represent the fundamental niche of Norway spruce, i.e., the potential distribution of the species being restricted by climatic constraints. For the reference period 1961–1990 it covers the greatest extent, with large parts of this potentially inhabited area being located outside the natural distribution. Height growth performance is particularly high in the center, while it decreases towards high elevations and the margins. Areas above 60°N are characterized by a sharp growth decline (<70 cm), while countries along the east side of the Baltic sea (Estonia, Latvia, and Lithuania) and large parts of Central and Eastern Europe (Poland, Germany, Czech Republic, Slovakia, Hungary, and Belarus) are predicted to have superior growth (>200 cm). With changing climate conditions there is a substantial decline in height growth potential across the study area, particularly for areas of low elevation. Former disjunctions (e.g., the Polish disjunction) do again become a spruceless belt. The range with superior height growth >200 cm retreats towards the natural distribution and further within (from 0.87 × 10⁶ km² for 1961–1990 to 0.08 × 10⁶ km² by 2080s, −91%, Table 5). This also applies for the southeastern part of the Baltic part of the species distribution, while some small improvements in height are predicted for areas above 60°N. Higher elevations, e.g., in the Alps or Carpathians will profit from climate warming by an increased growth potential.

Choosing “optimal” provenances has a positive effect on height growth (Figure 4). Towards the second half of the century, potential height gains are particularly prevalent at high latitudes. Compared to the first half of the century, total area where height gains may be achieved will then decrease (Table 5), but the absolute gain will increase (Figure 4). By the 2080s, it will be possible to achieve superior height growth on an area of 0.39 × 10⁶ km², if provenances were chosen according to the derived model. However, these areas are mainly restricted to the higher elevations and the north.

For the ensemble OccurrenceSDM, the optimal threshold (based on true skill statistics) to translate occurrence probability to binary presence/absence corresponds to a probability of 0.412. Growth potential predicted by the TraitSDM is likewise translated using a productivity-oriented threshold (surface shown in Supplementary Figure 5), which assigns a specific threshold value to each grid cell according to its latitude and elevation.

The largest suitable habitat according to the OccurrenceSDM is predicted for the reference period 1961–1990, covering an area of 3.95 × 10⁶ km², that is used as common baseline scenario (Table 6). During this time period, the TraitSDM identifies an area of 5.24 × 10⁶ km² as suitable habitat for “local” provenances, which is 15.9% more than the baseline scenario. The 10.3% range loss compared to the OccurrenceSDM do represent mainly subalpine and alpine zones in the Alps that become unsuitable due to their low productivity and at the south-east facing range margin in Russia, which the TraitSDM predicts to run slightly shifted northwards. On the other hand, it indicates a range gain of 26.2%, with a potential spread to the Massif Central in France and further south to the Pyrenees, as well as some additional habitat along the Dinaric Alps and the Eastern Carpathians. With forthgoing climate change the border of suitable Norway spruce habitat climbs up in elevation, as growth potential increases.

Range size of both OccurrenceSDM and TraitSDM decrease over time, facing the same direction with a retreat towards the higher elevations and the northern latitudes. The former does already decrease by 18.3% in the 2020s (Table 6), while the latter has approximately the size of the baseline scenario to then pass the baseline extent shortly after. According to the OccurrenceSDM Norway spruce will lose 39.8% of its current range by 2050s and as much as 55.3% by the 2080s. This equals to an area that remains suitable of approximately 1.77 × 10⁶ km². The retreat of suitable area proceeds likewise according to the TraitSDM with a remaining suitable area of 2.56 × 10⁶ km² in the 2080s (−43.3%). With only 1% range gain for each time slice predicted by the OccurrenceSDM, the species will not be able to expand to new suitable habitats, as those are, for example, restricted by the northern edge of Scandinavia.

If there had been a model-based provenance transfer during stand establishment in the past using “optimal” provenance, the species range might have been even 24% larger than the baseline scenario 1961–1990. Choosing “optimal” provenance for future stand establishments may lead to substantial increase in height growth potential (Figure 3), however, this has only a marginal effect to slow down species range change. By the 2080s the area potentially occupied by Norway spruce will also be 39% smaller (approx. 2.75 × 10⁶ km² remain suitable).

Site climatic variables contribute 82.7% to partial R². This is a strong indicator for a particularly high phenotypic plasticity in height growth of young Norway spruce. The two most important drivers of the plastic response identified by the TraitSDM are the frost-free period (FFP_s) and the summer heat-moisture-index (SHM_s). The quadratic pattern with SHM shows great growth declines towards both ends of the quadratic response within the climate range of observed sites. Low values of SHM generally represent cold and moist conditions, while high values imply warm and dry conditions. In essence, increasing temperatures paired with decreasing water availability, particularly during the summer, do strongly restrict future growth potential of Norway spruce (Figure 3). The consequences of this effect were already noticed in the species’ reaction to drought events, where it shows low resistance and resilience (Vitali et al., 2017). According to projections of future climate, SHM will greatly increase in the second half of the century, with a negative effect on Norway spruce, as conditions will surpass the maximum of the quadratic relationship (SHM_max = ∼45) after which height growth will drop quickly. Sites with an SHM exceeding 85, e.g., large parts of the North European Plain, become completely unsuitable for Norway spruce. Kapeller et al. (2012) likewise identified a heat moisture index, i.e., the annual heat moisture index, as best variable to predict tree height applying response functions for aggregated clusters of Austrian provenances. They hardly found areas throughout Austria, where conditions will be already too warm and too dry for Norway spruce by the 2080s. This is corroborated by our observation that the alpine mountain ranges, including Austria, will remain suitable for Norway spruce into the 2080s.

For these areas, FFP_s has an additional positive effect. With an increasing length of the frost-free period, height growth potential will increase, which is particularly noticeable at the higher elevations and the northern latitudes. Biologically this may relate to bud phenology in spring, which is largely driven by heat sum accumulation in Norway spruce (Hannerz, 1999). If the frost period ends earlier and conditions get warmer, required heat sums are fulfilled earlier and buds break. As budset on the other hand is mainly driven by the photoperiod, which remains constant, the length of available growing season between budbreak and budset increases, resulting in a higher growth potential. Considering the climatic range of the test sites included in our analysis, Norway spruce has just reached the upper limit at which a further extending frost-free period may not provide further benefit.

The observed dominant effect of environmental variables on the variation in height growth has been reported previously for Norway spruce (Chakraborty et al., 2019a,b) as well as for other conifer tree species (Wang et al., 2010; Rehfeldt et al., 2014; Chakraborty et al., 2015). Despite this, these studies still identified a noticeable effect of local adaptation to provenance climate. Chakraborty et al. (2015) could attribute a contribution of 18% of the observed variation in height at age 24 to the mean annual temperature of provenance for Pseudotsuga menziesii, and Wang et al. (2010) attributed 21% to the same climatic variable for Pinus contorta. In our TraitSDM of Norway spruce, provenance climate was found to have only a minor influence (partial R² = 2.2%). Provenance latitude, however, had a much larger effect (partial R² = 14.5%) and may be interpreted as a sign of local adaptation as well.

The geographic latitude of the provenance closely resembles patterns of local adaptation in the northern part of the distribution where latitudes are >55°N (Supplementary Figure 6), but fails to picture the high observed variation in height growth for the southern part of the range. We interpret this as an effect of the post-glacial migration history of the species. Fossil pollen (Latałowa and van der Knaap, 2006) and genetic analysis [RAPD (Collignon et al., 2002); mtDNA (Tollefsrud et al., 2008)] suggest postglacial colonization from at least three geographically distinct refugia in the present-day area of Moscow, the Carpathians and the Dinaric Alps. From these starting points, the so-called Baltic, Hercyno-Carpathian and Alpine domain colonized the European continent. Tollefsrud et al. (2008) explain the high population differentiation along the pronounced latitudinal gradient in the Baltic domain as a result of post-glacial migration in form of a front-like spread at high population density, paired with repeated bottlenecks resulting in low within population diversity. Considerable differences in the annual amplitude of the photoperiod (up to 8 h between 55° and 65°N) strongly shaped natural selection along the way (Collignon et al., 2002). The southern range of the species is dominated by the Central European mountain ranges that compromise migration. The Western Alps are characterized by decreasing genetic diversity at relatively short distance, as high mountain ridges and dry valleys acted as physical barriers and competing species introduced from the north (Fagus) and south (Abies). In contrast to this effect stand genetic patterns in the West Carpathians and the Bohemian Massif with low differentiation among, but high within population variation, which are most likely the result of admixture with high levels of geneflow between different colonization routes (Tollefsrud et al., 2008).

Despite these strong effects of postglacial migration mirrored in the model by the geographic latitude, the continentality (TD_p) appeared to be the most important provenance climate variable in predicting juvenile tree height. It significantly interacts with the frost-free period, which basically represents a form of genotype × environment interaction. Populations originating from milder climate (smaller difference between mean warmest and mean coldest month temperature) grow better at sites with a longer frost-free period, while those of harsher climate are better suited for sites with a short frost-free period, e.g., in Central Sweden.

Generally speaking, the TraitSDM estimates the fundamental niche, the full range of environmental conditions that a species can occupy and use, without any other limiting factors such as competition, biotic threats, or dispersal barriers. The size of the predicted fundamental niche of Norway spruce appears quite large with 7.29 × 10⁶ km² (Table 5 and Figure 4), and largely exceeds the suitable habitat predicted by the OccurrenceSDM. However, the OccurrenceSDM basically predicts the realized niche, making this direct comparison inappropriate. This has been recognized previously by Chakraborty et al. (2019b), who used true skill statistics to identify a single value as optimal height threshold segregating presence from absence. While this approach worked well in their case for Pseudotsuga menziesii, an introduced species in Europe and mainly distributed at the low elevations, it was not expedient here for Norway spruce, a species naturally associated with the harsh, cold, and moist conditions of the higher elevations and latitudes. By the use of a regression surface as threshold (Supplementary Figure 5) to approximate the realized niche, we were able to account for the different site potential and, naming just one example, could successfully identify Northern Scandinavia as suitable habitat despite the species’ low growth potential in this area.

Comparing the approximated realized niche of the TraitSDM with that of the OccurrenceSDM, the latter does underestimate the range Norway spruce may actually inhabit thanks to its high phenotypic plasticity. This high plasticity, paired with the ease of propagation, made the species the most preferred option to reforest degraded areas over the last two centuries (Jansen et al., 2017), which has significantly enlarged the species range. Still, our TraitSDM identifies an area of +16%, where there might have been further potential for range expansion in the past (Table 6). Obviously, Norway spruce will not realize this further range expansion any more now. However, when it comes to the future, the range decline may be delayed. According to the OccurrenceSDM Norway spruce will lose as much as 55% of its range by the 2080s, while it will lose only 43% according to the TraitSDM. This difference translates to a remaining suitable area of 2.56 × 10⁶ km² instead of only 1.77 × 10⁶km².

Despite this difference, both models do agree on a substantial range decline, first back towards the extents of the natural distribution and then further within. The general direction and spatial pattern of this decline to the north, i.e., towards the Boreal forest, and the high elevations closely resembles the predictions from other SDMs, regardless their use of different sets of predictor variables, global climate models, or even species occurrence data (van der Maaten et al., 2017; Dyderski et al., 2018; Chakraborty et al., 2021).

Despite the fact that a broad application of TraitSDMs across species would provide great benefit to support forest management decisions, they are limited by the availability of provenance trials spanning large geographic gradients (Urban et al., 2016). It is crucial to use data of the same age, which greatly restricts their application. To date, very few comparable studies were able to assemble even aged datasets of age 20 (Wang et al., 2010) or older (Chakraborty et al., 2015), while most studies use observations available from younger age (Leites et al., 2012; Gárate-Escamilla et al., 2019; Hu et al., 2019). In our case the most abundant comparative data is height at an age of 11 years, which is relatively low compared to a full rotation of Norway spruce between 50 and 120 years (rotation length largely depends on site potential, planting material, and management recommendations; e.g., Felton et al., 2017). However, juvenile growth is considered to be closely connected to fitness, as intraspecific competition among trees is strongest among seedlings and juvenile trees (Kapeller et al., 2017). Further, Norway spruce shows strong age-age correlations, which makes it a suitable candidate for early selection (Isik et al., 2010). Bałut and Sabor (2002) found stable provenance reactions between 11 and 25 years at the polish site of the IPTNS 1964/68 series. And at an unthinned site in Germany, provenance ranks were found to be consistent for an even longer time period. Correlations between height of age 10 and DBH of age 51 – an age difference of 41 years – still had a coefficient of 0.75, while they ranged between 0.83 and 0.98 for shorter time intervals (Liepe et al., 2020).

Norway spruce has been largely translocated during past reforestation waves, far beyond the natural distribution, but likewise within (Jansen et al., 2017). Artificial populations were included in the presented trial series on purpose, as it was intended as an inventory to “test as many provenances as possible, regardless whether the seed sources were autochthonous or not” (Krutzsch, 1974). This likely introduces some bias to the model as artificially translocated populations may have only adapted for one generation to the climatic conditions at seed source location. While it is impossible to identify artificially planted populations within the natural distribution, those located outside may be spotted in the data and removed easily. However, when fitting a TraitSDM using a dataset reduced by those artificial populations, we identified the same climatic predictor variables (analysis not shown), and therefore do present here the analysis using the complete set of provenances.

The two statistical models used in the present study only consider climatic effects as potential drivers and delimiters of height growth potential and species occurrence. However, it is expected that other disturbance agents will strongly limit successful cultivation of Norway spruce. Drought stress causes sharp decreases in vitality (Marini et al., 2017), significantly lowering tree resistance against the colonization by bark beetles (mainly Ips typographus) (Schuldt et al., 2020). These do by themselves already profit from warmer temperatures with an earlier propagation in spring and subsequent exponential increase in population size. In addition, devastating wind throw events are expected to become more frequent (Gardiner et al., 2010), but their timing and strength are unpredictable.

Apart from climatic or climatically induced effects, population demography, structure and dynamics, dispersal, and competition with other species significantly shape future species distributions. Process-based models that simulate these aspects based on observations on individuals in natural populations or experiments have been proposed as biologically more realistic (Morin et al., 2007). However, these models are very “data-hungry,” limiting their broad application (Guisan and Thuiller, 2005). In case of our study species, Takolander et al. (2019) observed a high accordance between predictions of a process-based dynamic vegetation model and those of a bioclimatic envelope model. Further, their predictions closely resemble our predictions, allowing the assumption that the OccurrenceSDM reflects the processes responsible for niche realization close enough.

In the past, Norway spruce was planted far beyond its ecological niche limits (Figure 2, 1961–1990). In the future, the suitable range will decrease drastically towards the natural distribution and further within (Figures 2, 3 and Table 5). Norway spruce will be limited to a restricted area, remaining mainly at the high elevations and latitudes, which basically represent its core habitat in the boreal and subalpine conifer forests (Schmidt-Vogt, 1977; Caudullo et al., 2016). Here, form the practitioners point of view, Norway spruce may actually profit from the effects of climate change with better growth potential, as growth is promoted by the extending frost-free period, and there is still room for improvement in terms of an increasing summer-heat-moisture index (cf. Kapeller et al., 2012).

The presented TraitSDM was also used to identify “optimal” populations solving the partial derivatives for the provenance variables in the model. The height gain indicated in Figure 4 quantifies the potential merits of assisted migration, which has been proposed as key forest management strategy to mitigate negative impacts of climate change on the one hand, while making use of the potential opportunities for increased growth associated with warmer climate on the other (Pedlar et al., 2012; Williams and Dumroese, 2013).

Solving the model for latitude of the provenance, one particular area is identified to provide optimal growth across planting sites, i.e., the northern edge of the Hercyno-Carpathian domain. Provenances from this area do actually frequently stand out by their superior growth vigor, e.g., on sites of the presented IPTNS 1964/68 trial series (Liesebach et al., 2001; Ujvári-Jármay et al., 2016; Liepe et al., 2020), but likewise in numerous other provenance trials (e.g., Rau, 1983; Schmidt-Vogt, 1986; Giertych, 2007; Zeltinš et al., 2019). Further the interaction between provenance continentality and frost-free period affects the choice of the “optimal” provenance, in a way that provenances from milder, less continental climate than currently observed at the planting location are assigned as the frost-free period extends.

Our TraitSDM does suggest a movement of provenances northwards to enhance growth in the northern part of the distribution. Such a movement has already been successful in the past, for example, with provenances translocated from the Carpathians to southern Sweden achieving superior height and diameter in these new environments (Milesi et al., 2019). We do, however, caution this movement to span too large latitudinal gradients, despite the predicted high potential for improvement. The highly consistent shape of individual transfer functions across sites up to a latitude of 60°N, does not apply to the northern most site in Lappkojberget (63.4°N, Supplementary Figure 6). At this site, local provenances appear to be better adapted to the harsh climatic conditions. For the southern part of the range provenance choice may further enhance growth in those mountainous areas, where the increasing frost-free period already shows some positive effect. In terms of range decline, the choice of optimal provenances for reforestation may have some positive effect (total range loss of 39% compared to 43% using “local” provenance), but it may not significantly change the general pattern (Figures 3, 4).

Even though bark beetles have caused vast devastations, e.g., during the extraordinary dry summers 2018–2020, Norway spruce may still present a preferable option in certain areas, if some means of caution are applied. For example, in the German Harz mountains or the Sauerland, where the OccurrenceSDM already indicates great range loss by the 2050s, it may still be possible to make use of the densely emerging natural regeneration, as height growth potential is predicted to remain high in the 2050s (Figures 3, 4). As long as the conversion towards less susceptible mixed forests is initiated to mitigate the unforeseeable risks and increase stand resistance, one short generation of Norway spruce may be exhausted as productive element drawing upon her economic value (Paul et al., 2019). Heterogeneity (instead of largely common homogeneity in species and age in the past) appears to rule the future. This is recognized by the forest practice: management prescriptions urge for species mixtures, uneven-aged structures, small scale patches, or a mixture of the different approaches.