Continent-Wide Tree Species Distribution Models May Mislead Regional Management Decisions: A Case Study in the Transboundary Biosphere Reserve Mura-Drava-Danube

Abstract

The understanding of spatial distribution patterns of native riparian tree species in Europe lacks accurate species distribution models (SDMs), since riparian forest habitats have a limited spatial extent and are strongly related to the associated watercourses, which needs to be represented in the environmental predictors. However, SDMs are urgently needed for adapting forest management to climate change, as well as for conservation and restoration of riparian forest ecosystems. For such an operative use, standard large-scale bioclimatic models alone are too coarse and frequently exclude relevant predictors. In this study, we compare a bioclimatic continent-wide model and a regional model based on climate, soil, and river data for central to south-eastern Europe, targeting seven riparian foundation species—Alnus glutinosa, Fraxinus angustifolia, F. excelsior, Populus nigra, Quercus robur, Ulmus laevis, and U. minor. The results emphasize the high importance of precise occurrence data and environmental predictors. Soil predictors were more important than bioclimatic variables, and river variables were partly of the same importance. In both models, five of the seven species were found to decrease in terms of future occurrence probability within the study area, whereas the results for two species were ambiguous. Nevertheless, both models predicted a dangerous loss of occurrence probability for economically and ecologically important tree species, likely leading to significant effects on forest composition and structure, as well as on provided ecosystem services.

1. Introduction

Human-driven environmental impacts caused by climate change and globalization threaten forest ecosystems [1]. This may lead to reduced habitat suitability for local forest communities [2], increased risks of climate interactions between native pests and diseases [3,4], and intensified disturbance events such as storms [5], droughts [6,7,8], and floods [9,10,11]. Globalization also contributes to the accidental spread of non-native species [12,13,14]. In forests, this can result in an increased incidence of non-native pests and diseases [15], to which native species often lack resistance [16], and of non-native plants [17,18] that cause competitive pressure on native plant communities [19].

Highly dynamic forest communities with special ecological processes and versatile characteristics, such as riparian forests, are especially vulnerable to environmental changes and external impacts. Major man-made environmental changes such as large-scale river regulations [20] have serious ecological up- and downstream effects [21,22] and power plants interrupt the downstream supply of gravel. More than 1.2 million barriers are built in European rivers, which fragment river systems and hamper connectivity [23]. In combination with increased flow velocity and truncation of sediment transport, they cause a rapid river bed incision and lowering of water levels [24,25,26,27]. In addition, groundwater extraction affects the water availability [28]. These modifications to the water supply have crucial long-term effects on natural floodplain forests [29]. Historically, agricultural cultivation is estimated to have already destroyed 90% of Europe’s floodplain forests over the last century [28]. As a result, manifold ecosystem services [30] provided by riparian forests are threatened; this includes the role of riparian forests as natural buffers and for protection against floods [31,32], the maintaining of high levels of biodiversity [33,34], the provision of recreational and aesthetic values [35], and high forest productivity [36]. Consequently, climate change’ morphological and hydrological modifications [20,26]; and the spread of non-native plants [35], pests, and pathogens [37,38] can have massive ecological effects and fundamentally influence the provision of ecosystem services [21,39,40,41].

According to the current state of knowledge, the selection of tree species and provenances are among the key measures in adaptive forest management under climate change and can help to increase their resilience and resistance concerning biotic threats in order to maintain and restore ecosystem services [42,43,44]. In this regard, species distribution models (SDMs) have increasingly emerged as an important tool to support decision making in forest [45,46] and conservation management [47,48]. SDMs correlate species occurrence with environmental data to predict occurrence probabilities. In most cases, the occurrence data originate from national forest inventories [49,50,51,52]. In combination with other meteorological and soil data, the potential distribution can be modeled in time and space. Multiple future climate scenarios are commonly used to address the uncertainty of future climate developments [2,53]. The increasing availability of high-resolution climate data, wide-range distribution data and computing power make SDMs a popular tool for mapping species-specific suitability. SDMs have mostly been applied to exploring the future suitability of species and developing policies to increase the resilience and carbon storage capabilities of forests. SDMs are also suitable tools for estimating potential forest cover for landscape restoration projects, as one of the most effective strategies to mitigate climate change [54].

However, SDMs are predominantly based on bioclimatic variables [48,55] and rarely on other environmental predictors such as soil parameters [56,57,58,59], land cover [60,61], or geomorphological light availability [62,63]. Even though contemporary studies have found that including environmental variables was beneficial [57,58,59,64] and reviews likewise recommend their use, such data are often neglected [62,65]. Furthermore, SDMs have rarely been applied on a local scale as specific site conditions can strongly influence the predicted suitability of climate SDM models [62] and high-resolution data on a local scale are often a key constraint. There is thus a need to broaden the field of SDM applications towards the prediction of local variations in species occurrence that could be explained by additional predictors at a fine scale [66]. It is particularly difficult to select proper predictors and develop accurate SDMs for riparian tree species, as the respective forest communities occur in very limited areas and the occurrence of the species is strongly related to riparian water influence, which is rarely represented in environmental data sets. However, predictions of riparian tree species distributions are important for forest restoration of riparian habitats, which are among the most endangered and the most diverse terrestrial ecosystems worldwide, providing a multitude of ecosystem services [67,68].

The aim of our study was to understand the spatial patterns and the major drivers for the distribution of seven riparian tree species in Europe at the regional and continent-wide scale. We developed two types of species distribution models (SDMs) incorporating climate, soil, and river data to identify expected regional as well as large-scale changes to species distributions caused by environmental changes. Furthermore, we discussed the benefits of bioclimatic SDMs on a continent-wide scale and of SDMs enriched with soil and river data on a local scale in the context of forest management and restoration.

2. Materials and Methods

2.1. Study Area

Our study area encompasses riparian forest ecosystems in central to south-eastern Europe, namely, forests within the transboundary Biosphere Reserve Mura-Drava-Danube (TBR). Since this study area is too small to study large-scale climate effects on species distribution, models were established on two levels. Firstly, we demonstrated large-scale bioclimatic effects on riparian tree species distribution on the spatial scale of Europe as a whole (−10.62 to 44.24° E longitude and 34.56 to 71.19° N latitude). Secondly, the specific ecological niche including soil and river variables were investigated in regional models for an area of 8300 km² along over 500 km of major European rivers in the TBR (15.647 to 19.448° E longitude and 45.165 to 46.765° N latitude, Figure 1). The TBR, which is still under development, aims to become an international model area for sustainable regional development. It consists of four spatially connected approved Biosphere Reserves in the countries of Austria, Slovenia, Hungary, Croatia, and Serbia. New parts of the TBR were recently nominated to be approved as the first five-country TBR in the world by UNESCO in 2021. The entire core zone of this important ecological corridor—a belt of riparian forests along the three rivers Mura, Drava and Danube—has been designated as protected areas of various categories. A share of 27% (2250 km²) of the TBR is covered by forests, increasing to 61% within the core zone.

The TBR is characterized by highly fertile plains along the rivers with intense agricultural use for cereal and pasture farming on the one hand, and forestry on the other. The rivers are embedded in eutric Fluvisols (33%), surrounded by Luvisols (14%) and Cambisols (5%). Across the entire TBR, Phaeozem (35%) is the dominant soil type. The annual mean temperature ranges from 9.3 °C in the north-western part of the study area to 11.7 °C in the area between Đurđevac (Croatia) and Barcs (Hungary). The TBR exhibits a gradient from an Illyrian climate in Styria and Slovenia to Pannonian conditions in the east, resulting in a strong variation of average annual precipitation, ranging from sites with nearly 1000 mm in the west to less than 500 mm in the north-eastern Hungarian part of the BR.

2.2. Tree Species

We selected the seven most important riparian tree species occurring in south-eastern Europe, which cover an ecological gradient from soft- to hardwood riparian forests and are important for both nature conservation and timber production. The target species were black alder (Alnus glutinosa (L.) Gaertn), narrow leafed ash (Fraxinus angustifolia Vahl), European ash (Fraxinus excelsior L.), black poplar (Populus nigra L.), pedunculate oak (Quercus robur L.), field elm (Ulmus laevis L.), and European white elm (Ulmus minor L.).

2.3. Species Occurrence Data

Different occurrence data sets were used to construct the two types of models. The continent-wide models were based on the EU-Forest tree occurrence dataset, which was derived from national forest inventories and featured a spatial resolution of 1 km [69]. We employed 66,060 occurrence points from this dataset after the data preparation. For the regional models we collected 269,460 spatially precise occurrences at the tree and stand levels directly from national forest inventories and forest enterprises in the countries of Austria, Hungary, Serbia, and Slovenia. Stand-level data were transferred to the polygons’ centroids. To overcome a lack of information for Austria and Croatia, we employed approximately 2500 georeferences from field observations carried out across the TBR in 2019. As recommended by Titeux et al. [70] for restricted area models, the applied data cover the entire territory of the countries and thus originate from beyond the study area to include conditions not yet present but possible in the future. After data preparation, we were able to employ 236,604 occurrences in total for the regional models (

Table A1. Number of employed occurrences for developing the two types of models of the seven selected target species. The continent-wide models were trained with occurrences and climate data only, whereas the regional models were based on precise occurrences on a smaller spatial scale and a larger environmental dataset of climate, river, and soil variables.

Model Type	Alnus glutinosa	Fraxinus angustifolia	Fraxinus excelsior	Populus nigra	Quercus robur	Ulmus laevis	Ulmus minor
continent-wide	9488	1046	17,596	2115	30,700	171	1836
regional	39,918	4821	45,445	1929	51,539	816	1952

).

The statistical software R version 3.6.2 [71] was used for data formatting. Both models were run with occurrence data and pseudoabsence data, created by selecting background points at a certain geographic distance from the presence points. This strategy avoids the selection of absences too close to the observed presences and therefore such of the same niche and pseudo-replication. For the continent-wide model we applied a buffer of 2 degrees [72,73] and 0.007 degrees for the regional model. Following Barbet-Massin et al. [74], we randomly selected an equivalent number of pseudo-absences outside the buffers with presences from other species to achieve the highest possible model accuracy. We removed all points with one or more missing environmental variables.

2.4. Environmental Data

To explain the tree species occurrences, we took variables from three grid-based datasets; of these, only climatic data were used for the continent-wide model, whereas all three types were employed in the regional model. We used the geographic information system QGIS 3.4.14 [75] and R [71] for the preparation of geographic data.

(1) Biologically relevant bioclimatic variables were obtained from the ECLIPS 2.0 dataset [73]. These variables comprise annual trends, seasonality, and extremes of temperature and rainfall derived from monthly records of precipitation and temperature. This dataset was developed from dynamically downscaled and bias-corrected regional climate model results from the EURO-CORDEX projections with a resolution of 30 arcsec. The dataset comprises 80 annual, seasonal, and monthly climate variables for past and future periods for two representative concentration pathways (RCPs) adopted by the Intergovernmental Panel on Climate Change [76]. RCP 4.5 was selected to reflect moderate climate change effects as it assumes that emissions will peak around 2040 and then decline, causing a rise in global temperature of 2.4 °C by the year 2100. RCP 8.5, however, considers severe climate change effects, in which the global temperature increases up to 4.9 °C if no adaptive measures are taken into consideration. For this study, the more extreme RCP 8.5 and RCP 4.5 predictions can be found in the Electronic Supplementary Material S1.

(2) River bodies are fundamentally important for riparian species. In order to characterize the influence of the river systems, we calculated the variables ‘relative elevation difference’ and ‘horizontal distance’ of a certain presence/absence location to the horizontally closest river shore or any other type of inland waterbody connected to a major watercourse. We prepared two raster layers by calculating the attributes for the center of each raster cell. To simplify the computing process while maintaining the explanatory power regarding riparian influences on tree species occurrence, each of those rasters with 50 m final resolution was assembled using two different resolutions. Within a range of 1.5 km along the rivers, we employed a 50-m grid size, whereas a grid of 1 km was applied outside this buffer and resampled after the calculation to adjust the resolution. Elevation data were taken from the European Digital Elevation Model (DEM), version 1.1, available at 25-m spatial resolution, which was obtained from the EU Copernicus program [77]. Vector data of the major river networks were received from the EU-Hydro River Network Database, Version 1.2 (EU-Hydro) [78] for the Danube river basin. These data provide precise contour lines of river and water bodies derived from photo interpretation of image data (2011–2013) at a resolution of 2.5 m.

For all rivers in our study area, a trend of river deepening has been identified (Mura: [26]; Drava: [24,25]; Danube: [27]). Therefore, we performed an additional case study (Section 3.5) to analyze the effect of further river deepening on the future tree species occurrence probability. We assumed a deepening of the rivers by 0.5 m throughout the study area, which might result in spatial relocations of potential habitats within the TBR.

(3) A subset of the soil data in the European Soil Database v2-Raster Library (ESDBv2 RL) [79], with a spatial resolution of 1 km, was derived from the European Soil Database v2 [80]. We focused only on 20 key variables within the primary categories of chemical, mechanical, and hydrological properties. Even though these soil variables might have matched the spatial resolution of the occurrence and bioclimatic data, we decided against using them in the continent-wide model in order to be able to compare a standard bioclimatic approach to the regional one based on more environmental predictors.

All variables applied in the variable selection process for both model types can be found in

Table A2. List of variables used in the variable selection process for the continent-wide and the regional models.

Continent-wide and regional model	AHM	Annual heat:moisture index
	bFFP	Beginning of FFP
	DDabove18	Degree-days above 18 °C
	DDabove5	Degree-days above 5 °C
	DDbelow0	Degree-days below 0 °C
	DDbelow18	Degree-days below 18 °C
	eFFP	Ends of FFP
	EMT	Extreme minimum temperature
	FFP	Longest frost-free period
	MAP	Annual total precipitation
	MAT	Annual mean temperature
	MCMT	Mean coldest month temperature
	MWMT	Mean warmest month temperature
	NFFD	Number of frost-free days
	PPT_at	Mean autumn precipitation
	PPT_sm	Mean summer precipitation
	PPT_sp	Mean spring precipitation
	PPT_wt	Mean winter precipitation
	SHM	Summer heat:moisture index
	Tave_at	Mean autumn temperature
	Tave_sm	Mean summer temperature
	Tave_sp	Mean spring temperature
	Tave_wt	Mean winter temperature
	TD	Continentality
	Tmax_at	Maximum autumn temperature
	Tmax_sm	Maximum summer temperature
	Tmax_sp	Maximum spring temperature
	Tmax_wt	Maximum winter temperature
	Tmin_at	Minimum autumn temperature
	Tmin_sm	Minimum summer temperature
	Tmin_sp	Minimum spring temperature
	Tmin_wt	Minimum winter temperature
Regional model	min_direct_dist	Horizontal distance to the next water body (major rivers and inland water bodies)
	diff_dem	Elevation difference between site and next river. Calculated from a 25 × 25 m Digital Elevation Model
	awc_sub	Subsoil available water capacity
	awc_top	Topsoil available water capacity
	bs_sub	Base saturation of the subsoil
	bs_top	Base saturation of the topsoil
	cec_sub	Subsoil cation exchange capacity
	cec_top	Topsoil cation exchange capacity
	dgh	Depth to a gleyed horizon
	dimp	Depth to an impermeable layer
	dr	Depth to rock
	eawc_sub	Subsoil easily available water capacity
	eawc_top	Topsoil easily available water capacity
	hg	Hydrogeological class
	oc_top	Topsoil organic carbon content
	pd_sub	Subsoil packing density
	pd_top	Topsoil packing density
	peat	Peat
	pmh	Parent material hydrogeological type
	str_sub	Subsoil structure
	str_top	Topsoil structure
	text	Dominant surface textural class (completed from dominant STU)

.

2.5. Data Analysis

2.5.1. Variable Selection

The univariate Spearman correlation was first calculated between explanatory variables (Table A2) and occurrence data (Table A1) and significant variables (p < 0.05) were identified. From the list of significant variables, those with a high multicollinearity (0.7) were excluded for each species, following Dorman et al. (2012). When two variables were highly correlated, the one with the higher Spearman correlation coefficient was chosen. This procedure resulted in two types of models:—the continent-wide model, using three to five climatic variables per species from the total list of 32 variables; and the regional models, using 20 to 24 predictors from the total list of 54.

2.5.2. Species Distribution Models

Both types of models (continent-wide and regional) were built using the machine learning package “randomForest”, version 4.6-14 [81] available in R [71]. Based on the employed environmental data, the spatial resolution was 1 km for continent-wide models and 50 m for the regional models. We predicted all models with RCP 4.5 and RCP 8.5 scenarios for the near-term (2041–2060) and long-term (2081–2100) period. To identify the main drivers for tree species occurrence in both models, we estimated the percentage increase in mean square error (InMSE) when a variable was dropped from the model and the other predictors were left unaltered. The number of trees to grow was set to 100. Further accuracy metrics were calculated using the “evaluates function“ in R package “sdm” version 1.0-49 [82].

Random forest estimates are not suitable for solving extrapolation problems due to their inability to project trends outside the ranges of the training dataset, in contrast to the standard regression procedure [83]. Instead of extrapolating, the model will always yield the closest output value covered by the training dataset. This also implies that, in the case of climate change projections, the reliable validity range of future maps is determined by the geographic extent to which the range of predictor variables at a future time overlaps with the range of predictor variables in the training period (present). All other geographic areas need to be handled more carefully, and we highlighted these areas in all predictive maps.

3. Results

3.1. Model Training Data, Model Performance, and Accuracy

Species distribution models were trained to elaborate the ecological niche of the target species by combining explanatory variables and occurrence data. Occurrence data across Europe cover a much wider climatic range (Figure 2a) for all species than the climatic spectrum of the TBR. However, the climate range also varies among species, with P. nigra and U. leavis revealing the smallest climatic amplitudes. The relocation of the TBR’s climatic envelope indicates that most of the species currently present in the area of the TBR are likely to be forced out from their macro climatic niche due to climate change.

Occurrence data for the regional models demonstrate the relationship with riparian habitats for P. nigra and U. leavis occurring only in close proximity to rivers, not exceeding a 1–3 km horizontal distance and a 10–20 m altitudinal difference. The other species are less bound to the immediate river and occur also at larger distances from the river or at higher altitudes (Figure 2b).

Overall, the variation explained by the continent-wide models ranged from 91% to 99%, with highest values for F. excelsior (99.1%) and Q. robur (98.0%) and the lowest values for F. angustifolia (92.5%) and U. laevis (90.6%). The variation explained by the regional models ranged between 53% and 81%, with F. excelsior (81.4%) showing again the best performance and U. minor (53.5%) and Q. robur (53.2%) showing the weakest. In both datasets, an increasing number of occurrences resulted in a higher share of variation being explained (Figure A1).

The accuracies of the regional models were slightly lower compared to the continent-wide models. The specificity and sensitivity of the regional models ranged from 0.89 to 0.95, whereas the continent-wide models ranged from 0.98 to 1.00. The true skill statistic (TSS) of the regional models ranged from 0.79 to 0.91, whereas it varied from 0.99 to 1.00 for the continent-wide models. For the regional models, the accuracy for species such as A. glutinosa and the Q. robur was slightly lower compared to other species. In case of the continent-wide models, the accuracy between the species does not differ significantly (

Table 1. Model performance and accuracy statistics for the continent-wide and the regional models. To validate species distribution model (SDM) results in terms of accuracy, the true skill statistic (TSS) was considered a straightforward approach to compare the model types. The proportion of correctly predicted presences (sensitivity) and absences (specificity) were obtained by comparing the predictions to the original dataset. TSS values range from −1 to +1 and are derived as: sensitivity + specificity − 1. The higher the output value, the better the SDM performance [85]. Cohen’s kappa coefficient is considered to correct for the model fit expected by chance [86] and is therefore frequently applied to support model accuracy evaluation [50,53].

Model Type	Species	Sensitivity	Specificity	TSS	Kappa
Continent-wide	A. glutinosa	1.000	1.000	1.000	1.000
	F. angustifolia	0.998	0.998	0.996	0.996
	F. excelsior	1.000	1.000	1.000	1.000
	P. nigra	0.999	0.999	0.998	0.998
	Q. robur	1.000	1.000	1.000	1.000
	U. laevis	1.000	1.000	1.000	1.000
	U. minor	0.999	0.999	0.999	0.999
Regional	A. glutinosa	0.895	0.895	0.790	0.790
	F. angustifolia	0.938	0.938	0.875	0.875
	F. excelsior	0.958	0.958	0.916	0.915
	P. nigra	0.943	0.943	0.886	0.886
	Q. robur	0.899	0.899	0.797	0.796
	U. laevis	0.958	0.957	0.914	0.914
	U. minor	0.908	0.908	0.814	0.816

).

3.2. Drivers for Tree Species Distribution

Across the seven continental models, the most important variable (InMSE ± SD) was the maximum autumn temperature (98.5% ± 39.5%). It was selected for four species and explained up to 68% of the observed variation. In addition, the predictors continentality (25.3% ± 3.5%), annual precipitation (24.4% ± 3.3%), and mean summer precipitation (21.5% ± 13.5%) were significant predictors in four models. The summer heat-moisture index (52.8% ± 4.2%) and the maximum summer temperature (43.2% ± 1.7%) were suitable parameters for two species models. Other climatic predictor variables selected for only one species or with less importance (<21% on average) were found to be degree-days below 0 °C, mean coldest month temperature, longest frost-free period, mean spring precipitation, mean winter precipitation, annual heat-moisture index, and maximum spring temperature (

Table A3. Number of employed occurrences in the continent-wide models, the models’ variation explained, selected predictors, and the corresponding percentage increase in mean square error (InMSE) when the variable is dropped from the models for the seven selected target species.

Species	Occurrences	Variation Explained (%)	Explanatory Variables: InMSE
Alnus glutinosa	9488	97.89	Tmax_at: 76.9	SHM: 48.4	MAP: 27.1	TD: 24.1	-
Fraxinus angustifolia	1046	92.53	Tmax_at: 42.5	TD: 25.3	MAP: 23	PPT_wt: 15.2	-
Fraxinus excelsior	17,596	99.10	Tmax_at: 144.5	SHM: 56.8	PPT_sp: 9.8	-	-
Populus nigra	2115	95.21	Tmax_sm: 44.9	DDbelow0: 21	PPT_sp: 20.1	PPT_sm: 14.0	-
Quercus robur	30,700	98.04	Tmax_at: 122.1	PPT_sm: 44.6	TD: 30.6	MAP: 27.9	-
Ulmus laevis	171	90.57	MCMT: 17.7	PPT_wt: 12.7	AHM: 12.7	Tmax_sp: 12.3	PPT_sm: 10.9
Ulmus minor	1836	95.54	Tmax_sm: 41.5	TD: 20.6	MAP: 19.5	PPT_sm: 16.7	FFP: 16.0

).

Given the narrow climatic spectrum of the regional occurrence data set, climate variables explained only 21.8% ± 2.5% of the observed variation (34.5% ± 5.3% InMSE). In contrast, soil explained 28.7% ± 6.3% (44.4% ± 4.4% InMSE) and river proximity explained 13.6% ± 3.8% (21.0% ± 4.1% InMSE) of the variation.

The highest mean cumulative InMSE (±SD) was found for the soil variables (151.6% ± 55.4%), which included 18 potential predictors in total. The group of climate variables contained eight potential predictors, with a cumulated mean InMSE of 120.4% ± 62.1%, followed by the river variables (75.6% ± 42.7%).

The modeled species occurrence of F. excelsior and P. nigra were found to be mostly driven by the soil variables (41.3% and 32.2%), to which Q. robur occurrence was least sensitive (19.3%). The models for U. laevis and F. angustifolia were strongly affected by climate predictors. Compared to the other species’ models, P. nigra and F. excelsior were most strongly affected by the group of river predictors (20.3% and 18.2%), to which U. minor was least sensitive (9.0%,

Table A4. Number of employed occurrences in the regional models, the models’ variation explained, selected predictors, and the corresponding percentage increase in mean square error (InMSE) when the variable is dropped from the models for the seven selected target species.

Species	Occurrences	Variation Explained (%)	Explanatory Variables: InMSE
Alnus glutinosa	39,918	61.66	diff_dem: 78; AHM_199110_1km: 61.8; DDbelow0_199110_1km: 49.4; DDabove5_199110_1km: 42.2; dr: 18.2; hg: 17.7; pmh: 17.4; awc_sub: 15.5; cec_sub: 14; bs_top: 13.5; awc_top: 12.1; oc_top: 11.5; pd_top: 10.6; bs_sub: 10.5; str_top: 10.5; text: 10.1; cec_top: 9.3; pd_sub: 8.9; dgh: 6.9; peat: 6.2;
Fraxinus angustifolia	4821	64.56	DDbelow0_199110_1km: 39.5; AHM_199110_1km: 37.6; min_direct_dist: 25.8; DDabove5_199110_1km: 25.6; diff_dem: 23.4; hg: 15.9; awc_sub: 11.6; cec_sub: 10.5; cec_top: 9.8; awc_top: 9.7; text: 9.7; dr: 9.4; dgh: 8.7; str_top: 8.3; bs_top: 7.4; pd_sub: 7.2; pd_top: 7.1; dimp: 5.3; str_top: 4.2; peat: 4.1;
Fraxinus excelsior	45,445	81.38	AHM_199110_1km: 57.0; min_direct_dist: 50.1; DDbelow18_199110_1km: 42.8; diff_dem: 32.5; awc_sub: 16.9; cec_sub: 16.4; cec_top: 16; hg: 13.4; text: 12.2; pd_sub: 11.9; pmh: 11.4; dr: 11.2; bs_top: 11.2; awc_top: 10.9; bs_sub: 10.9; dgh: 10; oc_top: 8.9; str_top: 8.7; pd_top: 7; str_top: 6.7; peat: 3.7;
Populus nigra	1929	71.26	min_direct_dist: 50.5; diff_dem: 26.5; DDabove18_199110_1km: 25.2; DDbelow0_199110_1km: 24.7; AHM_199110_1km: 21.8; pmh: 10.9; hg: 10.4; cec_top: 10.4; text: 9.7; str_top: 8.9; pd_sub: 7.6; bs_sub: 7.1; awc_top: 7; dr: 6.8; bs_top: 6.8; pd_top: 6.7; cec_sub: 6.6; oc_top: 5.9; dgh: 5.7; dimp: 5.6; str_top: 3.4; peat: 3.1;
Quercus robur	51,539	53.21	DDabove5_199110_1km: 89.2; EMT_199110_1km: 88.9; min_direct_dist: 87.1; diff_dem: 82.4; AHM_199110_1km: 79.5; awc_sub: 20.5; dr: 19.4; str_top: 17.3; cec_sub: 17.2; text: 17; hg: 16.2; bs_top: 15.9; awc_top: 15.6; bs_sub: 15.5; pmh: 15.4; oc_top: 15.2; cec_top: 13.3; dgh: 12.3; pd_top: 11.3; str_top: 9; pd_sub: 6.8; peat: 4.4;
Ulmus laevis	816	63.57	min_direct_dist: 20.7; DDbelow0_199110_1km: 20.2; DDabove18_199110_1km: 19.4; AHM_199110_1km: 13; FFP_199110_1km: 10.2; diff_dem: 7.3; hg: 5.8; text: 5.4; dr: 5.2; awc_sub: 4.7; oc_top: 4.5; awc_top: 4.4; cec_sub: 3.8; pmh: 3.4; cec_top: 3.4; bs_top: 3.1; str_top: 3.1; dimp: 2.9; pd_top: 2.8; dgh: 2.8; pd_sub: 2.2; bs_sub: 2.1; str_top: 1.2;
Ulmus minor	1952	53.25	MAP_199110_1km: 38.5; DDabove5_199110_1km: 29.8; DDbelow0_199110_1km: 26.9; diff_dem: 22.6; min_direct_dist: 22.4; hg: 10.7; dr: 10.7; pmh: 10; awc_sub: 9.7; bs_top: 8.9; bs_sub: 8.7; str_top: 8.5; text: 8.2; oc_top: 7.9; pd_top: 7.3; cec_top: 6.9; cec_sub: 6.7; pd_sub: 6.4; awc_top: 5.5; dgh: 4.9; str_top: 4.2; peat: 0.9;

, Figure 3).

3.3. Comparing Model Type Predictions

Among the seven riparian forest tree species across Europe, six species are expected to expand their species range significantly under climate change. Increasing mean probabilities of occurrence show gains between 6% (A. glutinosa) and 27% (U. minor). The species P. nigra, F. angustifolia, and U. minor are expected to gain additional habitats in Western, Northern and Eastern Europe without larger losses in Southern Europe. F. angustifolia and Q. robur are expected to gain additional areas in north-eastern Europe, but at the same time, a loss of climatically suitable habitat in Southern and Western Europe was shown. In contrast, U. laevis is expected to lose significant parts of its current distribution in Central Europe and showed an overall decrease in its probability of occurrence of 9%. Its distribution range is expected to shift to north-eastern Europe and will be significantly smaller than the current one (Figure A2).

At the regional scale of the TBR, the continent-wide model generally predicted drastic future losses for all seven species. The regional models predicted a fine-scaled, more complex occurrence pattern. In contrast to the continent-wide models, the regional models were able to depict the regional core areas of tree distribution around major riparian landscapes by including a higher number of limiting environmental variables (Figure 4 and Electronic Supplementary Material S1). Future species distributions of both model types were more similar to each other than the current predictions. For six species, current and future probabilities of occurrence of the regional model were significantly lower than those of the continent-wide model. On the contrary, U. laevis displayed higher future probabilities in the regional compared to the continent-wide model (Figure 5). For U. leavis, the continent-wide model showed a strong decreasing trend across central and south-eastern Europe as well as within the Biosphere Reserve. The regional model confirmed this trend but still predicted future habitats along the rivers (Figure A3).

Patterns of the regional model for high future habitat suitability were in line with the continent-wide model but varied locally, as different predictors were used. Differences between the models were demonstrated by low or negative correlations for most species. Diffuse changes across the area resulted in low correlations between both models for many species (

Table A5. Correlation coefficients between the continent-wide models and the regional models for (a) current predictions, (b) future predictions, and (c) differences between current and future predictions calculated for 6000 random points within the transboundary Biosphere Reserve Mura-Drava-Danube.

Correlation	Alnus glutinosa	Fraxinus angustifolia	Fraxinus excelsior	Populus nigra	Quercus robur	Ulmus laevis	Ulmus minor
(a) Current	−0.17	0.71	−0.44	0.07	0.27	−0.45	−0.04
(b) Future	0.29	−0.33	−0.03	−0.12	−0.18	0.03	0.08
(c) Change	−0.02	0.64	0.54	0.10	0.14	0.39	−0.15

).

3.4. Predicted Co-Occurrence of Species

To investigate potential co-occurrence between the seven species, we estimated Pearson’s r correlation coefficient of occurrence probability for the current and future scenarios (

Table A6. Pearson’s r correlation coefficient for the predicted occurrence probabilities between the seven tree species at 6000 randomly distributed points across the transboundary Biosphere Reserve for the current and future scenarios of the regional model.

). The range was −0.64 to 0.80. Applying a minimum threshold of 0.2 for the correlations, three groups of spatially associated tree species were derived: (I) F. angustifolia, P. nigra, U. laevis, and U. minor; (II) A. glutinosa, and F. excelsior; and (III) F. angustifolia, Q. robur, and U. minor. Under expected future conditions, the correlations of group (III) decreased below the threshold and instead two new groups evolved: (IIIa) A. glutinosa, F. angustifolia, and Q. robur; and (IIIb) A. glutinosa, and Q. robur. Changes in habitat suitability will unbalance the competition and are likely to affect the future mixture and functioning of forest stands.

3.5. Consequences of River Deepening for Habitat Availability

To investigate the potential impact of river deepening on the future tree species occurrence, we tested the effects of river deepening by assuming a drop of 0.5 m throughout the study area. Despite the importance of the variable elevation difference to rivers for all species, only two tree species were significantly affected in occurrence probability distribution for the TBR in the models. The habitat of F. angustifolia was found to decrease by 5% occurrence probability on average and by 6% for the third quartile, which better covers the actual distribution. This may add to the decreasing occurrence probability induced by climate change. For the mean occurrence probability, the effect is similar to 102% of the magnitude of climate change until 2100 and for the third quartile to 36%. A. glutinosa showed the opposite trend. By applying the changed elevation differences, the average occurrence probability within the TBR increased by 3% and mitigated the predicted negative trend of climate change (Figure A4 and Figure A5). The predicted shifts are visible in the spatial pattern (Figure 6). These changes in occurrence probabilities are strongly site-dependent, but the observed trends of both species are also valid across the five countries.

4. Discussion

4.1. Consequences for Managing Riparian Forests of the TBR

Climate change is likely to force the currently present species out of their macro climatic niche. Our regional models showed a decrease in occurrence probability within the TBR for five of the seven species, whereas F. excelsior and U. minor could benefit. Paradoxically, these two species are already heavily threatened by severe diseases: F. excelsior by ash dieback, which is caused by the fungal species Hymenoscyphus fraxineus [87], and U. minor by Dutch elm disease, which is caused by the tubular fungus Ophiostoma novo-ulmi [88]. Overall, our results are in good agreement with studies analyzing the velocity of climate change, which showed that south-eastern European lowlands face among the highest velocity climate change across the continent [89,90] and a European-wide analysis of winner and loser species in conservation areas showed that our study region contains a high share of loser species in expected climate change [91].

The expected decreases in habitat suitability will not only impact the (potential) habitats of native plants, but can be expected to facilitate in particular a further spread of non-native plants [92,93,94]. This is because riparian forest development is strongly driven by regular flooding regimes [95], which drive small-scaled mortality with subsequent canopy gaps and irregular light regimes [96]. Additionally, river networks are essential vectors [97] and many non-native tree species are typical wetland species of riverine environments in their country of origin. Stimulated by these conditions, and in interaction with the prevailing propagule pressure, alien trees can be established [98]. Therefore, riparian forests have already developed as hotspots for this group [99,100]. With the (ongoing) dieback of trees or the disappearance of certain tree species, canopy gaps will become more frequent and environmental niches will become unoccupied, thereby further increasing susceptibility to invasion by non-natives.

To avoid a loss of forest functions, forest and conservation management needs to (i) evaluate the species’ potential to adapt to the changing environment, (ii) evaluate the artificial transfer of forest reproductive material from matching locations to speed up adaptation processes and reduce adaptation lags [101], (iii) rely on the potential of the other remaining native tree species to fill ecological gaps, or (iv) fill empty niches by means of assisted migration with new tree species which have desired characteristics. All passive measures raise the risk of uncontrolled penetration by tree species of an invasive character. Conservation management is faced with a decision between accepting such dynamics or counteracting them in advance or ex post facto. There have been positive and negative experiences with the active introduction of non-native tree species to riparian forests, such as Platanus orientalis, Juglans nigra, and Robinia pseudoacacia, as well as several non-native Populus species and their hybrids [102,103].

In contrast to our expectations, the deepening of rivers affected only two species considerably. We conclude that the models may be not sensitive enough to elevation modifications and highlight that soil parameters connected to water supply have not been modified. The periodic flooding favors riparian species adapted to such disturbances. However, it can be assumed that the models did not explicitly consider the suppression of riparian river dynamics, for instance by channel modifications or flood prevention dams, which can influence habitat suitability negatively.

4.2. Comparison of the Model Results and Predictor Importance with Other Studies

Our approach for the regional models targeted an optimal grain size as a trade-off between the necessity to describe riparian influences and the resolution of a forest stand level and model efficiency. The more realistic regional models generally explained less variation than the continent-wide models, in good agreement with previous studies [104,105].

The comparison of our future occurrence maps to existing SDM predictions was limited to the continent-wide model and a few species, as our target species are rarely investigated. Regarding Q. robur, very similar future distribution patterns were discovered by Schueler et al. [90] and Dyderski et al. [2]. Discrepancies were especially found in the southern European range, where our predictions indicate potential occurrences in the continental subtropical climates. Furthermore, a similar trend but with a more sensitive distribution pattern was predicted for the period 2061–2090 by Buras and Menzel [49], applying the same RCP scenarios. For F. excelsior, our prediction is comparable to Dyderski et al. [2] in RCP 8.5 across the full range, except in parts of the Scandinavian peninsula, where our occurrence probability prediction was significantly lower. Thurm et al. [51] modeled the future distribution of U. laevis using climate and soil data and found a significant increase in its distribution area to 227% across Europe up to 2070 by applying RCP 8.5. In contrast, our continent-wide model prediction of the same climate change scenario indicated a distribution change to the north-east, with major losses in Europe up to 2100. Encompassing further soil and river variables, our regional model also predicted much lower occurrence probabilities for the current and future climates. Reasons for these differences are related to a different variable set, differing origins of the input data, or the assumptions made in the modeling approach. More generally, such differences highlight the need for better datasets of species occurrences and environmental drivers for rare tree species with specific habitat requirements.

Regarding the importance of the model predictors in our regional models, soil variables were found to be the most important group, in agreement with other studies applying such environmental predictors [56,57,59]. Climate and river distance variables were found to have a lower impact. Velazco et al. [58] showed that on a worldwide scale the use of additional edaphic predictors increased SDM performance, and we concluded that soil variables could also increase the predictive ability of continent-wide models if aimed at this objective. Similarly to our regional model, a model for a mountainous landscape in Switzerland including climate and soil variables [59] found soil variables to be the most important for F. excelsior. By contrast, climatic variables for Q. robur were only slightly more important in our model, whereas they clearly dominated in the study by Walthert and Meier [59]. This is likely due to the different regional characteristics—Walthert and Meier [59] studied the distribution of Q. robur in Switzerland, including the cold-induced edge of its distribution, whereas our study analyzed oak populations close to the warm end of their current potential distribution. The variables of degree days were considered most important within the climatic group in the regional models for F. angustifolia, P. nigra, Q. robur, and U. laevis in contrast to the continent-wide model, where the maximum temperatures were the most important for six of the seven species. The importance of degree days was already highlighted by Dubuis et al. [106], who showed this predictor to be the most important in their models applying soil and climate data at a regional scale in Switzerland.

4.3. Data and Model Constraints

The distribution of tree species has been studied and modeled throughout the last decades [107,108,109], often focusing on zonal climax forest communities reflecting regional climate norms [59,109]. However, due to the limited availability of specific environmental predictors for their occurrence, riparian species have not been investigated to our knowledge.

Our models were based on two types of tree occurrence data sets and different environmental data, and they featured different grain sizes. The employed EU-Forest occurrence dataset includes the edges of continent-wide tree species distribution and is widely accepted for SDMs, but is limited to a resolution of 1 km. In contrast, collecting precise occurrence data as used in the regional model can be time- and cost-intensive, but strongly contributes to more realistic fine-scaled predictions. The applied bioclimatic data were similar for both developed model types. Superior to other available climate datasets, it provides 1 km spatial resolution (interpolated). The European Soil Database v2-Raster Library [79], with a spatial resolution of 1 km, matches the bioclimatic data, but cannot compete in accuracy with the other employed data and might limit the accuracy of our models. Forest communities can change within a few meters laterally or few decimeters vertically to the rivers, and with the regard to the precision of the occurrence, river, and elevation data, more accurate models could be built. Additionally, geographic proximity to a river increases the probability of flood events, with a limiting effect for many tree species. Consequently, our two river variables demonstrated an explanatory power nearly equivalent to bioclimatic variables over all investigated species. As the interaction between river level change and soil (water) parameters is complex, we did not modify the latter for our river level change scenario. In our models we could not differentiate between active and former floodplains where flooding was restricted by structural measures, as the applied systems for the delineation are different across the five countries and river systems. Alternatively to basic DEMs, flood area models could provide more reliable predictors for modeling riparian tree species distribution. In our case, such data were not available for the extent of our study area as it covered five countries and three river systems. Despite the influence of rivers on soil parameters and artificial limitation of flood areas, including river data enabled us to increase the resolution and accuracy.

For species distribution modeling, it is also important to be aware that the locations and species compositions of forests in Europe strongly reflect the results of anthropogenic influences [110,111,112]. Much of Europe’s forest area has undergone a significant transformation throughout history. In Hungary, for example, the age structure of forests reflects the growing periods of pines and poplars until the last “era of black locust”. Therefore, based on current distribution data, the natural distribution potential of tree species can only be partially reconstructed, especially in lowland areas. In SDMs, these limitations must be considered in order to interpret potential distributions. Furthermore, river levels have already changed significantly in the past and due to time lags in forest composition these changes might not be mirrored in the observed tree occurrences with which our model was trained. Therefore, a further change in river depth of 50 cm, as assumed in additional predictions of this study, might have greater consequences for habitat suitability than those expressed in our mapping.

We did not include any land cover information [2,60,61] besides the topographic river derivations, as we intended to predict occurrence probability independently from the current use of riparian zones to enable application for restoration projects. As an example of including land cover for areal habitat assessments, we hereby refer to Figure 6, where we cropped our results to forested areas. Although studies have shown the importance of differences in light regime between north- and south-facing slopes for species occurrence [62,108], we did not include this variable since our study area consisted of typical lowland river basins without considerable slopes.

In the field of ecology and the case of black-box methods used for predictions, a validation of the models not only statistically but also by existing ecological knowledge is important. Validating the predictions revealed limitations in the models’ behavior and indicated constraints based on the variable selection process. In the case of P. nigra, the model did not really find a variable that would negatively affect the distribution of the tree species as climate change progressed and the projections mainly indicate area gains. Similarly, we found constraints for the regional model of Q. robur, with future probability losses at the rivers and gains outside this area where distance as a predictor no longer had a significant effect, and the increase in degree-days and the increase in the annual heat-moisture index values became dominant and increased the probability. Thus, the predicted slight increase outside the vicinity of rivers should be treated with caution. In the Electronic Supplementary Material S2, we describe the functional relationship between the important predictors and predicted occurrence probabilities for both model types and all species, facilitating the correct interpretation of the predicted maps.

4.4. Practical Application of the Models

The continent-wide model easily assesses the climatic niche for the continent-wide population and can depict climate-change-driven modifications of occurrence probability on a continent-wide scale. With this model being limited to climate predictors, it is very sensitive to expected climate change. Thus, compared to the regional model, the predicted impact of climate change is higher and the resulting occurrence probability is in fact the bioclimatic suitability. Therefore, such continent-wide models might overestimate real occurrence changes but can be recommended for early warnings regarding climate change.

The precise occurrences applied in the regional model allowed us to model spatially complex patterns at a significantly smaller grain size. The proximity to the river could be efficiently used to understand the drivers for local variation in tree occurrence. Furthermore, the regional model focused on a subpopulation of the trees adapted to specific environmental conditions and does not capture the full populations’ environmental range. The future environment might be unsuitable for the local population but could potentially match the requirements of other (sub)populations. Therefore, predicted mismatches are not equal to a low occurrence probability for the species in general and do not necessarily exceed the population’s ecological amplitude. Derived ecological niches need to be seen in this context and are not spatially transferable [48,113]. In contrast, a general problem that increases with smaller geographic scale is that changing climate conditions are not encompassed by the training data. Nearly 60% of over 200 SDM studies reviewed by Porfirio [114] assessed climate change impacts at the regional scale and were potentially forced to extrapolate, whereas only 14% considered continent-wide scales. Our investigations highlight this problem, driven by serious climatic changes both in RCP 4.5 and RCP 8.5 scenarios, especially for the regional model. The predictions indicate major alterations in species distribution even without the extrapolation of trends in occurrence probability, and we assume that extrapolation would have further increased the magnitude of predicted changes.

Nevertheless, the complex spatial pattern that can be realized with regional models reaching down to stand level resolution, as well as the opportunity to precisely match the spatial extent to certain purposes or populations, predestines them to support decision-making by forest and conservation managers. Especially in terms of adapting to climate change, such realistic predictions are of increasing value. For the operative use, standard large-scale bioclimatic models alone are too coarse and frequently exclude relevant predictors. We recommend using multiple model predictions with different scales and various predictor groups [115] and to consider the uncertainties related to the underlying model assumptions [48] (pp. 126–127).

5. Conclusions

This study provides SDMs for seven major tree species in riparian forests to guide forest and conservation management in the five-country transboundary Biosphere Reserve Mura-Drava-Danube (TBR). Climate change is likely to force present foundation tree species out from their macro climatic niches. A. glutinosa, Q. robur, and U. laevis are expected to substantially suffer from changing climate conditions, whereas F. excelsior and U. minor could potentially benefit by the end of this century. The results of this study can support decision-making at a local level such as in tree species selection [46], prioritization of conservation and restoration activities [48,116,117], and raising awareness regarding the impacts of climate change [2,50,53].