Climate envelope analyses suggests significant rearrangements in the distribution ranges of Central European tree species

Abstract

Key message

Climate envelope analysis of nine tree species shows that Fagus sylvatica L. and Picea abies H. Karst could lose 58% and 40% of their current distribution range. Quercus pubescens Willd and Quercus cerris L. may win areas equal with 47% and 43% of their current ranges. The ratio of poorly predictable areas increases by 105% in southern and south-eastern Europe.

Context

Climate change requires adaptive forest management implementations. To achieve climate neutrality, we have to maintain and expand forest areas. Impact assessments have great importance.

Aims

The study estimates the potential climate envelopes of nine European tree species for a past period (1961–1990) and for three future periods (2011–2040, 2041–2070, 2071–2100) under two emission scenarios (RCP4.5 and RCP8.5) based on the current species distribution.

Methods

Climate envelopes were estimated simultaneously using the random forest method. Multi-resolution segmentation was used to determine the climatic characteristics of each species and their combinations. Models were limited to the geographical area within which the climatic conditions correspond to the climatic range of the training areas.

Results

Results showed remarkable changes in the extent of geographic areas of all the investigated species’ climate envelopes. Many of the tree species of Central Europe could lose significant portions of their distribution range. Adhering to the shift in climate, these tree species shift further north as well as towards higher altitudes.

Conclusion

European forests face remarkable changes, and the results support climate envelope modelling as an important tool that provides guidelines for climate adaptation to identify threatened areas or to select source and destination areas for reproductive material.

1 Introduction

Climate change significantly alters the geographical distribution of forests worldwide. Droughts, which have become more frequent in Europe (Spinoni et al. 2015), not only negatively affect growth and production (Mátyás et al. 2018; Bodribb et al. 2020; Schuldt et al. 2020), but may also cause increases in mortality (Allen et al. 2015; Cailleret et al. 2017; Zscheischler and Seneviratne 2017; Buras et al. 2018; Choat et al. 2018; Zscheischler et al. 2018). Tree species have enough time to adapt to the conditions and retain their habitats as long as the climate is stable and there are no rapid trend changes or extreme fluctuations in weather patterns. In addition to climatic conditions, tree species have been able to adapt to other environmental factors such as edaphic conditions prevailing in their habitats over the centuries. Demonstrating sufficient resistance and competitive advantage against competitive, predatory, parasitic or pathogenic species present in the same areas has ensured the prolonged habitat presence of many tree species. However, climate change is threatening this status quo and putting selective pressure on populations through both abiotic and biotic factors. Consequently, tree species spatial distribution and forest composition are expected to change (Fekete et al. 2017; Scherrer et al. 2017). The distribution of economically important native European tree species may decrease or rearrange significantly due to climate change (Hanewinkel et al. 2013). The renewal and rotation period of European tree species that are important to timber production is usually slow and long (often more than 100 years). Compared to this, the climatic change is much faster for tree species to respond (Corlett and Westcott 2013). Most species will be unable to manage this relatively short and quick transition through local adaptation or natural migration. Consequently, forests need to be adapted to changed conditions by altering management strategies or by planting more adaptive or better pre-adapted and, thus, more resilient tree species (Chakraborty et al. 2015; Halofsky et al. 2018; Sousa-Silva et al. 2018).

Assessing the future vulnerability of dominant tree species in Europe is crucial due to their high economic and ecosystem service value. Unfortunately, the vitality of many tree species is likely to decline in the future (Walentowski et al. 2017), which has prompted researchers to actively investigate the extent to which other tree species or provenances with greater climate resistance could become suitable alternatives. Numerous studies have provided valuable insights into the potential response of selected tree species to climate change. Some of these studies have focused on tree-ring analysis (Scharnweber et al. 2011; Lévesque et al. 2013; Zang et al. 2014; Dulamsuren et al. 2017; Rehschuh et al. 2017; Buras et al. 2018), while others have investigated common garden experiments (Huang et al. 2017; Mátyás et al. 2021) or species distribution models (Walentowski et al. 2017; Thurm et al. 2018; Buras and Menzel 2019).

According to Booth (Booth 2017), the majority of studies that deal with species distribution modelling aim to capture the distributional changes based on the natural distribution of species—the realized niche. Studies that extend natural distribution data with outer locations of occurrence are rare. By outer locations we mean areas outside the natural range where the species was introduced for, e.g., production purposes. Using all occurrence data of the species would provide an opportunity to map the adaptation capabilities of species more realistically and, thus, provide a better estimation of fundamental niche. The rarity of such studies is due to the difficulty of data acquisition.

Most species distribution models use presence data as input parameters and very few refer to any confirmed absence data that are very difficult to obtain. Nevertheless, the majority of models (general linear models, classification tools, entropy models) require so-called pseudo-absence points that are inferred from available information about the presence locations of the species using different methods (Barbet-Massin et al. 2012; Dyderski et al. 2018; Thurm et al. 2018; Higgins et al. 2020). Remarkably, even fewer studies used exclusively presence distribution points (e.g., Buras and Menzel 2019). Recently, Buras and Menzel (Buras and Menzel 2019) have modelled the species composition changes on European forests using the method of climate analog areas combined with ensemble of downscaled climate projections and a European forest inventory (Strona et al. 2016; Mauri et al. 2017) albeit by using climatic data of low spatial resolution.

In our approach, we used the same idea to identify regions in space and time that have more or less similar climatological characteristics. These climatically similar regions likely provide suitable growth conditions for the tree species located in a particular region.

Thus, instead of (pseudo-) absence-presence based distribution models focusing on one tree species, the potential distribution of several tree species can be estimated simultaneously, and the absence data representing an uncertainty factor can be replaced by the effect of interactions between tree species. By analyzing the climatic niches of certain tree species for the past and the future we can determine and designate source- and destination regions in which the reproductive materials from the source regions can be utilized in the destination regions. Various international projects have increased their efforts in this area: (e.g., SUSTREE (https://www.interreg-central.eu/Content.Node/SUSTREE.html), REFOCUS (Sallmannshofer et al. 2021)).

In our research, we analyzed how the distribution of potentially suitable climatic areas of nine widely distributed and important European tree species and their combinations may change due to climate change over three periods (2011–2040, 2041–2070, and 2071–2100) and two emission scenarios (RCP 4.5 and RCP 8.5) using the method of climate envelopes. The aim of the present study is to outline the estimated climate envelope shifts based on the European distribution of the studied tree species, and to identify areas where the expected future climatic conditions are not covered by any existing European provenance. A novel approach (spatial multi-resolution segmentation) was used to determine climate analog areas by considering a total of 36 bioclimatic variables with high spatial resolution. Simulated tree species also include less frequently modelled tree species, such as the Hungarian oak (Quercus frainetto, Ten.), which had not been included in European distribution model studies before. Another novelty of our approach is that it examines not only the occurrence of individual tree species, but also the co-occurrence of tree species combinations.

The research seeks to answer the following questions: (1) how will the species-specific climate envelopes change in the future? (2) Are there currently any forest areas where none of the studied tree species is expected to be stand-forming in the future? (3) Could the selected tree species (area changes) be complementary to each other?

2 Materials and methods

In our modelling approach, we used only two basic datasets, which we describe in detail in the following sections. We used a dataset for species distribution data for Europe originating from a well-structured data harmonization action. We used another dataset for bioclimatic variables. This dataset also originates from a well-described climate database, and we used it as a base for defining climate envelopes and to serve as predictors for future scenarios.

2.1 Selected species and their distribution data

For the study, we selected nine tree species including six late-successional deciduous trees that play a dominant or co-dominant role in the closed associations of the large forest cycle in European temperate forests (Kuuluvainen 2016). These trees possess suitable characteristic traits (longevity, slow growth, large seed, and longer regeneration cycle) for such a role.

Besides the deciduous trees, we also selected three conifers that have a dominant role and high economic importance from forestry point of view. The selected conifers are widely distributed in Europe thanks to their ecological and economic characteristics. All the selected species are important for the European forestry sector and are significantly exposed to climate change impacts (Czúcz et al. 2011). Table 1 lists the selected species and their reported share of European forests.

Table 1 List of selected species

To generate current distribution maps of the tree species, we used the European tree occurrence dataset (Strona et al. 2016), which unifies species occurrence data for numerous European countries into a homogenised dataset.

2.2 Species occurrence data processing

From the European species distribution database, we selected training and test points to characterize the distribution of each species (Fig. 1A, B). The selection was completed by assigning a random number between 1 and 100 to each point. All occurrence points that had an assigned random number greater than 75 were included in the test database, while the others were among the training areas. We made rasters from the selected training points of each tree species with the same grid layout as of the bioclimatic variables—described later in this chapter. This was necessary because the European species distribution database is gridded to fit with the INSPIRE 1 km × 1 km grid, which does not overlap with the 1 km × 1 km grid of bioclimatic variables that we used. The rasters had unique values regarding those cells that contained the given tree species' occurrence points, while the cells without occurrences had the value of 0. All nine rasters of the tree species were combined into a final raster that contained the aggregated presence of individual tree species or their combinations in its cell values—in the cases of overlapping occurrences. (The detailed description of the handling of species combinations can be found in Appendix 1) Figure 1 demonstrates that the random selection kept the spatial characteristics of occurrence points, including density or spatial extent.

Fig. 1

Selected training (A) and test (B) samples of Quercus petraea, Liebl. from the European species distribution database

2.3 Bioclimatic variables

We used the Climate EU dataset (Marchi et al. 2020a; Marchi et al. 2020b) with a cell resolution of 1 × 1 km including 36 bioclimatic variables. The spatial extent of the database covers Europe between 34.26° and 71.24° degrees latitude and − 10.74° and 44.24° longitude. The bioclimatic rasters represent climate means for 1961–1990 and were used to calibrate species-specific climate envelope models: annual mean temperature (MAT), seasonal averages of precipitation, seasonal maximum, minimum, and mean temperature, as well as the average temperature of the warmest month (MWMT) and the average temperature of the coldest month (MCMT). In addition, the temperature difference (TD = MWMT-MCMT), the average annual precipitation (MAP), and the average amount of precipitation during the growing season (MSP-May to September) were included. The annual thermal humidity index (AHM = (MAT + 10)/(MAP/1000)) and the summer thermal humidity index (SHM = MWMT/(MSP/1000)) were taken into account. Degree metrics include the number of degree days of frosty days below 0 °C (DD < 0), the number of degree days above 5 °C or vegetation degree days (DD > 5), the number of degree days below 18 °C or heating degree days (DD < 18) and the number of days above 18 °C (DD > 18) or cooling degree days. They also include frost-free days (NFFD), the length of the frost-free period (FFP), the start and end time (bFFP, eFFP) of the frost-free period, and the amount of precipitation falling as snow between August of the previous year and July of that year (PAS). Finally, the parameters included the extreme minimum temperature (EMT) for the 30-year period, the Hargreave reference evaporation (Eref), and the Hargreave climate humidity deficit (CMD). We did not use further features from digital elevation models because they are already involved in downscaled bioclimatic variables.

For future projections, we used the ensemble estimates of 15 Global Circulation Models (AOGCM) with the best validation statistics based on RCP 4.5 and RCP 8.5 emission scenarios for the 30-year periods 2011–2040, 2041–2070, and 2071–2100 also included in the Climate EU dataset (Marchi et al. 2020b).

2.4 Data processing and statistical evaluation

2.4.1 Climate analog areas

Instead of directly assigning the values of the bioclimatic rasters to the species occurrence raster data points, we proceeded to create a coherent polygon network from the bioclimatic rasters by multi-resolution segmentation (Darwish et al. 2003) under Trimble’s eCognition v8. In this polygon network, each segment (polygon) encompasses a group of pixels that is more homogeneous but different from its environment, while forming a unified system in their superimposed topologies of different resolutions. The number of segments covering Europe in the finest resolution was 405,301. The segments designated by the training areas served as the basis for the definition of climate envelopes for tree species. (The detailed description of the segment implementation can be found in Appendix 2) The advantage of segmentation is that it organizes raster information into larger spatial objects that can be used for characterization, such as local, regional, or “global,” depending on the scale. In the present case, the building blocks of climate envelopes are the segments. An additional advantage is that the computational capacity requirement is reduced by grouping the pixels into segments.

2.4.2 Prevalence of current and future tree species-specific climate envelopes

For tree species-specific distribution of climate envelope models, we sought to estimate the occurrence of tree species distribution simultaneously using the random forest method by treating the climate envelope models as a classification issue (Breiman 2001). Random forest is a statistical tool consisting of a set of independent decision trees that allow for predictions on the same data domain. The forest prediction is the majority vote of all trees. Overall, we chose the random forest method for our model because, in our experience, it is well suited for modelling spatial patterns in most cases. The random forest method can handle continuous and discrete variables at once. It can also manage category variables and provide robust results. After defining the training areas, we performed the calibration of the random forest model with the bioclimatic data of the reference period 1961–1990. In the random forest algorithm, 301 decision trees were used, and for the decision at the nodes, we used randomly selected predictors equalling the square root of the number of variables used as predictors. The model was run iteratively three times with 10 replicates or until a consistent result was achieved, i.e., until there was a difference between the results of two consecutive runs of the random forest model. Using the model, we classified all segments covering Europe into the tree species combinations. This resulted in a raw distribution map of species (Fig. 2A).

Fig. 2

Example of the estimation phases of the climate envelope for Turkey oak in the case of the reference period (1961–1990): (A) raw estimate of random forest model, (B) the potentially predictable area based on the parameter ranges of training areas, (C) the extrapolation-free real climate envelope estimation after clipping the raw estimate with the predictable area

This map could not yet be considered as the final result for the reference period as it also included extrapolated areas. By extrapolation, we mean that the prediction also included a return value for segments whose bioclimatic parameter values for the given tree species combination fell outside the range defined by the corresponding training areas. In response to this, we had to limit the estimates of our model to the range of climate envelopes of each tree species taken for bioclimatic variables and covered by the training areas. This was achieved by creating a mask layer per tree species, which—taking into account all bioclimatic variables—selected the spatial range in which the values of all variables remained within the limits designated by the training areas (Fig. 2B). This delimited—for each tree species—the space in which the model estimate for that tree species could be considered as a valid estimation. Restricting the area to this spatial range assigned by the model to the given tree species—either stand-alone or in combination—outlined the final geographic projection of climate envelopes for the tree species without extrapolations (Fig. 2C).

By merging the extrapolation-free climate envelope areas of tree species and their combinations, the realistic pan-European tree species (combination) distribution map was compiled for the reference period from the available data. Climate envelope maps for the other three future periods and two emission scenarios were created using segmentation of the corresponding bioclimatic rasters. After segmentation, we ran the random forest algorithm on each of the segmented map variants corresponding to the six period and climate scenario variations. For each of the six maps, the extrapolation-free climate envelope areas were determined for each tree species. Finally, we created the realistic estimable optimal climate envelopes using the area masks. (Rasters of estimated periods are available in the following dataset: Illés and Móricz 2022.)

2.4.3 Grouping of tree species and their combinations by characteristics

Since the 9 tree species and their combinations formed 30 separate classes, we considered it worthwhile to merge the classes into larger composites for the sake of clarity. These larger composites are based on the water demand of the tree species and their presence in typical forest types. In our opinion, the macroclimatic water demand and the typical forest types that can be assigned to tree species can be well used to interpret the impact of climate change. When defining the composites, we regarded the special, characteristic tree species as the determining factor. As an example, one extreme is the case of fresh mixed conifer-broadleaved forest, which is characterized by the combined presence of beech and spruce among our tree species. The other extreme is the case of xerophytic broadleaved forest, which is a combination of downy oak and Turkey oak among our tree species. The composition of the composites and their constituent classes is shown in Table 2. We examined how the area of the composites defined in Table 2 changed among the different scenarios.

Table 2 Classification of tree species and groups of tree species into composites based on the basic species characteristics, e.g., of being mesophilic, xerophytic, thermophilic, or cold tolerant.

2.4.4 Statistics for evaluation of mapping results

We used the following statistics to evaluate mapping performance based on a test dataset according to Wunderlich et al. (Wunderlich et al. 2019). We compiled the confusion matrix of our model based on the area of correctly and incorrectly classified segments. The applied statistics were as follows:

• Overall accuracy: (true present area + true absent area)/total area;

• Sensitivity: true present area/(true present area + false present area);

• Specificity: true absent area/(true absent area + false absent area);

• True skill statistics: sensitivity + specificity–1;

where true present area represents the area of test segments, which were assigned correctly to the species whose test points fell into the segments; true absent area represents the area of test segments where the model does not predict the given species and they do not contain points for the given species; total area represents the area of test segments belonging to the given species; false present area represents the area of test segments, which were assigned to a given species incorrectly; true absent area represents the area of test segments, which are correctly classified as absent areas of the given species; false absent area represents the area of test segments, which were assigned to be absent areas of the given species incorrectly.

3 Results

The generated maps were subjected to an accuracy test for the reference period. Based on the occurrence data separated from the distribution points for validation purposes, we checked how accurate the spatial estimation of our model was for the reference period (Table 3).

Table 3 Accuracy of the model for the investigated tree species

The accuracy statistics of the climate envelope maps in Table 3 demonstrated that our model fitted well with the available factual data for most tree species on average, albeit with varying performance. Overall accuracy, sensitivity, and specificity have the same range between 0 and 1. Pinus sylvestris, L. presented the worst accuracy and specificity values (0.84). For all the other species, these values reached or exceeded 0.9. Sensitivity exhibited weaker results, especially in Quercus frainetto, Ten., which had the smallest and most localized distribution area. In general, the more widespread a species is, the higher the sensitivity value the model performed. The values of true skill statistics interpreted between − 1 and + 1 were above 0.5 in all nine cases, above 0.6 in five cases, and above 0.7 in two cases, ranging from 0.51 to 0.86. This indicated that the model effectively captured the distribution datasets.

According to the moderate scenario (RCP 4.5), the potentially suitable climate envelope occupied by beech, Scots pine, and spruce is expected to decrease significantly (by more than 40%) in Europe over the century (Table 4). A smaller but significant reduction of the potentially suitable areas is expected for sessile oak and Hungarian oak (approximately 30%). For pedunculate oak and black pine, a smaller decrease of less than 20% was predicted. However, the potentially available climate envelope can increase significantly (> 40%) for Turkey oak and downy oak. It is worrying that in both scenarios, the proportion of precarious areas is increasing significantly (by more than 100%), i.e., areas whose expected climatic conditions are outside the currently known European distribution areas of all nine tree species. This poses a significant adaptation risk. This is discussed in more detail in connection with Fig. 5.

Table 4 Modelled current potential climate envelope area of tree species and its change according to the two emission scenarios

During the present century, the potentially suitable climate for the Turkey oak, black pine and downy oak, which is dominant in the Mediterranean areas, mainly occupies the hitherto climatic areas of the currently dominant tree species in Central Europe. Among these tree species, the downy oak can also be highlighted, which can occupy significant areas mainly in Western and Central Europe. In general, the potential climate envelopes for almost all tree species are shifting further north and towards higher elevations. An exception to this is the Hungarian oak, where only a shift towards mountainous areas can be observed during the century (Figs. 3 and 4).

Fig. 3

Projected potential climate envelope shift of some selected tree species under the emission scenario RCP 4.5. Green areas indicate identical areas under current and future predictions, blue indicates areas that are potentially occupied, and red indicates areas that are potentially unsuitable

Fig. 4

Projected potential climate envelope shift of some selected tree species under the emission scenario RCP 8.5. Green areas indicate identical areas under current and future predictions, blue indicates areas that are potentially occupied, and red indicates areas that are potentially unsuitable

The prevalence of the currently dominant tree species in the Mediterranean and the Balkans, and thus tree species diversity, may also decline. These areas can be considered as potentially highly vulnerable areas.

Figure 5 shows the change of tree species composites according to emission scenarios and time periods. In most cases, the tree species composites shift to the north and west, while the area of the fresh forest types visibly decreases.

Fig. 5

Maps of tree species-composites for present and future periods according to RCP 4.5 and RCP 8.5 emission scenarios. The species composites are according to Table 2

The most significant lesson from Fig. 5 is that the proportion of precarious areas is increasing significantly in southern and south-eastern Europe—partly covering Central Europe—i.e., areas in which expected climatic conditions are unlikely to have pre-adapted origins (provenances) within Europe. At least this appears to be the case for the nine tree species examined in our study.

Table 5 also numerically demonstrates the change predicted by the maps shown in Fig. 5. The data confirm that, regardless of the scenario, the proportion of areas favorable to heat- and drought-tolerant tree species could increase significantly. At the same time, the proportion of areas with climatic characteristics typical of fresh forest types is expected to decrease significantly. There is also a significant increase in the proportion of areas that can be estimated uncertainly, which is particularly the case in the southern and central European regions (Fig. 5) most exposed to the effects of climate change.

Table 5 Estimated future area changes of tree species composites (%) compared to the period 1961–1990

4 Discussion

4.1 Changes in the geographic extension of species climate envelopes

At this point, it is worth noting that the distribution data, which is the starting point for modelling, fundamentally determine the applicability of the modelling results. We consider realistic mapping to be an important issue—especially in the border zones of distribution areas. It is also worth noting that the actual potential climate envelope may be larger than the modelling results, as the modelling is based only on the distribution data in the databases; thus, areas that may be appropriate for the tree species could still exist. Booth (Booth 2017) suggested dealing with global distribution data instead of localized data. Our study supports this initiative. The results indicate the potential presence of large areas with climatic conditions not covered by European provenances.

Our species distribution projections suggest that the face of the European forests will alter significantly during the 21^st century. The determinant factor for a species to gain or lose distribution area under climate change is strongly related to their current geographical distribution. The tree species with the center of gravity in higher latitudes (e.g., beech, Scots pine, Norway spruce) shift further northward and to higher altitudes with decreasing distribution area, particularly for RCP 8.5. On the other hand, the resulting gaps are filled partly with the more southern tree species such as Turkey oak, black pine or downy oak that may gain areas in Central Europe. The projected shifts and changes in distribution resulted in increased tree-species diversity in most regions of northern central Europe and a decreased tree-species diversity in the Mediterranean.

It is generally difficult to compare such study results due to the difference in the applied scenarios and explanatory variables. We found similar potentially highly vulnerable areas among the studies that can be compared to ours (Buras and Menzel 2019). It appears that for the Mediterranean and south-eastern Europe—partly covering Central Europe—the number of climatically suitable tree species is likely to decline. Consequently, these regions should be given particular attention in the context of adaptation to the changing climate and tree species should be selected carefully by involving several potentially suitable species.

The projected changes in the distribution of tree species are largely in line with the latest research findings. An example of the drought sensitivity of Scots pine is the mortality observed in Iberia and Central Europe immediately following drought events in 2003 (Rebetez and Dobbertin 2004), 2005 (Galiano et al. 2010), and 2015 (Buras et al. 2018). Dyderski et al. (Dyderski et al. 2018) classified Scots pine as a “losing” tree species because it is projected to suffer one of the largest changes (increasing endangered area, shifting north), which is also confirmed by our research results (Fig. 4). Norway spruce, also a significantly affected tree species due to climate change, may disappear from lowland areas in Central Europe according to several studies (Dyderski et al. 2018; Thurm et al. 2018). This forecast is also consistent with our results (Illés and Móricz 2022). The vulnerability of conifers is exacerbated by the sizes of the occupied areas, which decreases to the north, and by their drought tolerance, which is far lower than that of deciduous tree species (Dyderski et al. 2018).

Dendro-ecological (Zang et al. 2014; Walentowski et al. 2017) and experimental comparisons of the climate vulnerability of Norway spruce and beech tree species (Pretzsch et al. 2020) showed that the drought tolerance of beech is higher and that beech regenerates better than Norway spruce after drought events. Dyderski et al. (Dyderski et al. 2018) classified beech as a “winning” tree species because the potential distribution area in Central Europe does not decrease significantly in the period between 2061 and 2080 according to their modelling using the average emission scenario. In contrast, Thurm et al. (Thurm et al. 2018) reported a decrease in distribution area, which is similar to our results (Table 4) (− 56% RCP 4.5 and − 70% RCP 8.5) for beech for the period 2061–2080. Numerous regional and European-scale studies have found that oaks display good resilience to climate change (Dyderski et al. 2018; Walentowski et al. 2017; Perkins et al. 2018; Thurm et al. 2018). In our research, the potential area of sessile oak and Hungarian oak decreased the most among oaks (Table 4, Figs. 3 and 4), which is probably related to their current geographical distribution and associated climate adaptation in higher latitudes compared to the southern tree species such as Turkey oak or downy oak. The potential distribution area of sessile oak shifts further north and towards higher hills and mountains during the century, with more drastic changes expected in scenario RCP 8.5 in a similar spatial distribution as in Sáenz-Romero et al. (Sáenz-Romero et al. 2017). At the same time, several studies have shown that sessile oak is expected to provide stable production by the end of the century, even in areas with drier climates, although vitality loss is certainly likely (Hlásny et al. 2011; Mátyás et al. 2018). According to our analysis, the potentially suitable climatic space of pedunculate oak is reduced to a lesser extent than that of sessile oak. It should be noted that a favorable result for pedunculate oak might depend on its habitats, which are partly located in semi-arid areas with good moisture supply, excess water (shallow groundwater). Climatic conditions alone would probably be insufficient to sustain pedunculate oak in these areas. Although pedunculate oak is less drought-sensitive than beech (Scharnweber et al. 2011; Walentowski et al. 2017), dendro-climatological research and distribution modelling studies show that this tree species will also be endangered in the future due to climate change (Árvai et al. 2018; Sallmannshofer et al. 2021). Several studies demonstrated that Turkey oak may be one of the winners of climate change (Führer et al. 2011; Hlásny et al. 2014), due to its high tolerance against drought (Nardini et al. 1999; Móricz et al. 2021). Thurm et al. (Thurm et al. 2018) show that the potential distribution area of Turkey oak may more than double, which is supported to a lesser extent by our research (Table 4, Figs. 3 and 4). In our research, the center of gravity of black pine spread will shift from the Mediterranean to Central Europe during the century (Illés and Móricz 2022), while its potential climate envelope is likely to decrease. At the same time, its migration to the north supports its sensitivity to drought and high temperatures (Linares and Tiscar 2010; Móricz et al. 2018). On the other hand, according to several recent European distribution modelling results, its potential distribution area will increase during the century, mainly under unfavourable soil conditions (e.g., Thurm et al. 2018; Buras and Menzel 2019). We found that the potential distribution area of downy oak to be similar to that of Turkey oak. Downy oak currently spreads from south to north in Europe. Thurm et al. (Thurm et al. 2018) affirms this by predicting an area increase of 70 and 120% for the period between 2061 and 2080, respectively, applying the RCP 4.5 and RCP 8.5 scenarios. Dendro-ecological research also confirms that the drought tolerance of downy oak is even greater than that of Turkey oak (Tognetti et al. 2008).

4.2 Model and prediction limitations

Although several studies use pre-sorting and filtering between variables (Sallmannshofer et al. 2021), mostly to filter out autocorrelation, we chose not use any variable filtering. One reason for this choice is that the random forest method chosen to create the models is less sensitive to autocorrelation between predictors (Dormann et al. 2013). Another reason for not filtering predictors was that autocorrelation is really a problem when an extrapolation issue needs solving. The chosen random forest predictive method is unsuitable for such a task anyway (Hengl et al. 2018). Therefore, in our predictions, we proceeded to exclude the areas affected by extrapolation.

Although the inclusion of soil data in the models may improve their performance, we decided not to use any large-scale soil datasets like SoilGrids (Hengl et al. 2014), or the harmonized soil database of FAO (Fischer et al. 2008) for the following reasons. Gridded soil information is produced via modelling that generally includes features from elevation models, climate coverages, or even vegetation-based predictor variables (e.g., evapotranspiration). This would result in barely manageable, hidden relations between predictors and target values. Additionally, the use of gridded soil data would propagate uncertainty measures into our models according to spatially changing prediction quality of soil models in various ways. Furthermore, climate change is necessarily causing changes in the soil, as well. However, unlike climate change data, there is no available dataset about soil change, which would result in one of the predictors in the model remaining constant over time. Finally, despite gridded global soil maps having a nominal resolution of 1 km or even 250 m, in reality extensive areas have the same soil group, or soil type value in these maps.

Although our model was trained with long-term climatic data (1961–1990), it inherently contains variability. The change of climatic means in the future does not imply constant climate variability. Due to climate change, climate variability is expected to increase, leading to an increase in extremes (Rajczak and Schär 2017). During the modelling, we considered changes in 30-year climatic averages; however, this does not consider the increasing variability (e.g., droughts). Thus, our forecast for the occurrence of tree species is probably optimistic. The studied tree species are likely more endangered than predicted, e.g., due to drought-induced mortality (Allen et al. 2015; Senf et al. 2020). Nevertheless, our models were based on climatic conditions corresponding to 1 km²; thus, in highly variable topographical areas it could not take into account the micro and mesoclimatic conditions due to the different slope and exposure, where tree species may remain as relics even if the macroclimate has shifted otherwise. It is also necessary to note that our climate envelopes may be smaller than they are in reality because the species are more widespread than described by our presence data. A narrower estimate may mitigate the effects of extremes, as they act more strongly at the edge of the distribution than in the core area of the species. Consequently, we did not consider a positive effect (higher real prevalence) and, in return, did not consider the negative effect of extremes.

Due to its nature, our model cannot take into account the possibility that the physiological characteristics of tree species may change in the changing environment. For example, we do not know how it would affect adaptation capabilities if water use efficiency increases for some tree species. In the same way, we cannot take into account in this model if the opposite happens, that is, if the deterioration of the efficiency of physiological functions is amplified by environmental changes.

Our model only reflects climatic changes when studying the distribution of tree species. However, it cannot take into account the impact of new biotic pests in a given area due to climate change, which can significantly threaten forest production, health status and ecosystem service capacity (Hicke et al. 2012; Klapwijk et al. 2013; Anderegg et al. 2015; Kern et al. 2021).

To avoid difficulties of ‘model-in-model’ approaches, we did not use outside Europe occurrence data or other datasets to maintain data quality or data homogeneity. During the modelling, we calculated the distribution area of tree species from European areas. Areas outside Europe or non-native tree species that may be of interest primarily to Mediterranean areas as potential alternatives in the future were not included. This is a particularly significant issue because the area comprising the Mediterranean region is excluded from the climatic range covered by the nine tree species we studied, even when examined under the RCP 4.5 scenario (Fig. 5). We will probably not find pre-adapted propagating material for the regeneration of forest areas for these nine tree species. The question of nature conservation may arise. If we cannot find native or at least European tree species to maintain forest areas, then to what extent is the use of non-native tree species an option? Furthermore, the migration of tree species is not only a function of climatic variables. Other factors such as unfavorable soil conditions and late frosts can also act as migration barriers for tree species.

Finally, the modelling was presumably negatively affected by the impact of human interventions (artificial regeneration, selection) on the current distribution of tree species, as there were very significant anthropogenic effects on the tree species composition of European forests (Kirby and Watkins 2015).

5 Conclusions and outlook

According to the modelling results, the potential distribution area of the currently dominant tree species is about to change significantly. The extent of the suitable climatic range of some of the major tree species, such as spruce, beech, and pine, may decrease significantly in the plains during the current century in Europe. Oaks, mainly downy and Turkey oak, but to a lesser extent pedunculate oak, can become potential alternative tree species whose potential range may largely remain or even increase in Europe, depending on which scenario becomes a reality.

Based on the distribution models, the actual distribution data, and the modelled climatic data an attempt can be made to examine which areas in Europe will have the climatic conditions that meet the climatic conditions of a species’ past distribution areas and vice versa. In the first case, the provenances of the given area can be examined in order to determine which areas the gene pool of these populations could be transferred. In the second case, we can solve the replacement of declining populations with more adapted propagating material of the same tree species in order to ensure the survival of forests in the given area. This probably also enables the preservation of the source populations' gene pool.

Identifying source and destination areas in the climate envelope is only the first step. The implementation also requires comparing these potential areas with actual or perceived distribution data for tree species. To determine the specific areas, it is also necessary to have the actually occurring and sufficiently old stands of the given tree species in the source and destination areas in order to transfer the propagating material. Therefore, the theoretically possible source and destination areas should also be compared with possible distribution data. As climate change progresses, modelled future distribution maps can provide clues as to which real occurrences are likely to be the most stable or long-term. These maps can also help determine the appropriate pace at which to move propagating material bases over time.

Appendices

Appendix 1

1.1 Details of species occurrence data processing and handling species combinations

While shaping the training dataset, we wanted to proceed in a way that separated the species combinations as much as possible while maintaining their characteristics to improve distinctiveness and prediction accuracy. We also sought to avoid overlap between the training areas of combinations. For the sake of manageability, we wanted to reduce the number of combinations to ensure there would be no more than approximately 30 separate combinations in the end. At the same time, we sought to maximize the number of training points in each combination. In order to reduce the number of combinations, our analysis included only the species combinations whose total extent reached or exceeded 100 km² in Europe (exceptions were made only in the case of the Hungarian oak (Quercus frainetto, Ten.) due to its fragmented occurrence). This constraint excluded 0.45% of the total study area and reduced the number of combinations by 63%. The number of species combinations meeting the above condition of spatial extent still reached 41. Out of these 41 combinations, eight contained only one tree species, while the other combinations contained two or more tree species. In the case of QF (Hungarian oak), we only had mixed classes, which were used further as a merged combination for QF. In order to avoid overlaps between the training areas of the different tree species combinations we placed the cells including a single species occurrence in separate combinations from the mixed cells, which contained their mixed occurrences with other species. Therefore, we treated the mixed species combinations separately from their individual cases. For instance, we considered that the pure occurrences of Quercus petraea, Liebl. (QP) and Quercus cerris, L. (QC), differ from their combination (QP-QC).

According to the above, where two tree species were included in the same grid cell, the cell was characterized by indicating both tree species present at that cell. This accounted for 28 combinations. We found 13 species combinations where more than two tree species were included in one grid cell. Combinations containing more than two tree species were drawn back to an existing combination of two tree species if it was possible; or they were transformed into new two-species combinations while keeping the character of the original combinations.

During the simplifications, we proceeded by giving preference to deciduous tree species as a first consideration. For timber production reasons, both pine and spruce had been planted a great deal in the past; this extended both species far beyond their natural range (Boratynski 1991; Zerbe 2002; Szabó et al. 2017; Barbati et al, 2007). Because, e.g., black pine was commonly used in afforestation target programs, the species was planted in areas where it would not otherwise have appeared (Csontos and Cseresnyés 2015). Such afforestation programs for industrial purposes did not take place for deciduous species; hence, we believe deciduous species display a more natural distribution pattern. Therefore, mainly deciduous tree species were retained within the deciduous-pine combinations if two of the three or more species were deciduous. If two of the three species were pine, we chose the ecologically more suitable pine species next to the deciduous species. In the case of three or more deciduous species, the ecologically closer species were selected in pairs. The ecological suitability of the pairs was assessed by expert judgement based on the co-occurrence of tree species, shade tolerance, and water demand. The simplification of combinations containing more than two species is described in Table 6.

Table 6 Simplification scheme of grid points containing combinations of at least three tree species

Appendix 2

1.1 Details of the implementation of multi-resolution segmentation

Figure 6 shows the Bükk Mountains and their narrower surroundings in Hungary. The mountainous segments are more diverse, while the lowland segments are larger. We examined the rasterized training occurrence points of tree species to determine in which individual segments they appeared. A segment containing training occurrence point(s) of only one species combination became a training area of the combination represented by that point(s). However, it often occurred that several tree species occurrence points fell into the same segment. If these points belonged to the same tree species combination, the situation could be easily handled by removing the extra points (since the entire segment remains part of the training areas, it does not matter which actual occurrence point is left in it). If the points belonged to different tree species combinations, the segment was assigned to the combination with the majority within the segment. That is, we aimed to minimize the overlapping training areas and to obtain segments that belong exclusively to one combination as the training area. If it was impossible to make a clear decision on the combination of a segment on a majority basis, then that segment was excluded from the training areas.

Fig. 6

Example of the segmentation map based on 36 bioclimatic variables for the reference period 1961–1990. Points within the segments indicate the rasterized species occurrences. The background is a digital elevation model (source: Institute of Geodesy Cartography and Remote Sensing (IGCRS)) with a resolution of 20 m. Labels are abbreviations for each tree species and their combinations (Table 1 and Table 6)

All segments assigned to one given combination provided the data for climate envelopes of the tree species based on the bioclimatic rasters of the reference period between 1961 and 1990. Finally, 43,062 segments from the 405,301 segments (10.6%) became training areas.

Using the pixel values included in the segments, statistical indicators of bioclimatic variables were calculated for both training area and non-training area segments. For continuous variables such as temperature variables, mean and standard deviation were used, while for discrete data types such as number of frost-free days, median and quantile values (25 and 75%, respectively) were considered. Segments devoid of any training points were later assigned to one of the tree species groups during a random forest classification.

The number of training areas, test areas and training segments by species are shown in Table 7. Figure 7 displays the training segments belonging to QP. A comparison with Fig. 1 demonstrates that the segmentation did not alter spatial characteristics of training points.

Table 7 The number of training points, test points, and training segments by species

Fig. 7

Map of training segments assigned to QP