Altitudinal shifting of major forest tree species in Italian mountains under climate change

Noce, Sergio; Cipriano, Cristina; Santini, Monia

doi:10.3389/ffgc.2023.1250651

ORIGINAL RESEARCH article

Front. For. Glob. Change, 08 September 2023

Sec. Forest Management

Volume 6 - 2023 | https://doi.org/10.3389/ffgc.2023.1250651

Altitudinal shifting of major forest tree species in Italian mountains under climate change

$\r\nSergio Noce,$ Sergio Noce^1,2^*

Cristina Cipriano^1,2

Monia Santini¹

¹Impacts on Agriculture, Forests and Ecosystem Services (IAFES) Division, Fondazione Centro Euro-Mediterraneo sui Cambiamenti Climatici, Viterbo, Italy
²National Biodiversity Future Center (NBFC), Palermo, Italy

Climate change has profound implications for global ecosystems, particularly in mountainous regions where species distribution and composition are highly sensitive to changing environmental conditions. Understanding the potential impacts of climate change on native forest species is crucial for effective conservation and management strategies. Despite numerous studies on climate change impacts, there remains a need to investigate the future dynamics of climate suitability for key native forest species, especially in specific mountainous sections. This study aims to address this knowledge gap by examining the potential shifts in altitudinal range and suitability for forest species in Italy's mountainous regions. By using species distribution models, through MaxEnt we show the divergent impacts among species and scenarios, with most species experiencing a contraction in their altitudinal range of suitability whereas others show the potential to extend beyond the current tree line. The Northern and North-Eastern Apennines exhibit the greatest and most widespread impacts on all species, emphasizing their vulnerability. Our findings highlight the complex and dynamic nature of climate change impacts on forest species in Italy. While most species are projected to experience a contraction in their altitudinal range, the European larch in the Alpine region and the Turkey oak in the Apennines show potential gains and could play significant roles in maintaining wooded populations. The tree line is generally expected to shift upward, impacting the European beech—a keystone species in the Italian mountain environment—negatively in the Alpine arc and Northern Apennines, while showing good future suitability above 1,500 meters in the Central and Southern Apennines. Instead, the Maritime pine emerges as a promising candidate for the future of the Southern Apennines. The projected impacts on mountain biodiversity, particularly in terms of forest population composition, suggest the need for comprehensive conservation and management strategies. The study emphasizes the importance of using high-resolution climate data and considering multiple factors and scenarios when assessing species vulnerability. The findings have implications at the local, regional, and national levels, emphasizing the need for continued efforts in producing reliable datasets and forecasts to inform targeted conservation efforts and adaptive management strategies in the face of climate change.

1. Introduction

The findings of the most recent Italian National Inventory of Forests and forest Carbon Pools (Inventario Nazionale delle Foreste e dei Serbatoi forestali di Carbonio INFC-2015, Gasparini et al., 2022), published in September 2021, highlight a consistent increase in the forested area across Italy. Italian forests now cover 15 million hectares, which accounts for about one-third of the national territory. This data emphasizes the crucial importance of these areas from a land management perspective. Following estimates by the Italian National Institute of Statistics (ISTAT), the forested area was 5.7 million hectares in 1954, shortly after the Second World War, when significant deforestation occurred due to energy and wartime demands. However, the forested area gradually expanded to nearly 9 million hectares by 1985, as documented by the first National Inventory of Forests and forest Carbon Pools, and has continued to grow since then. Italian forests are steadily increasing in surface, recolonizing areas that were previously abandoned due to human activities. Recent data from the Food and Agriculture Organization (FAO) in their Forest Resources Assessment (FAO, 2020) positions Italy among the top ten countries globally in terms of the rate of expansion of forested areas. Such a consistent increase is likely influenced by various factors, including changes in land use and management practices, but is in particular dominated by natural processes (Agnoletti et al., 2022). This trend can have positive implications for biodiversity providing new habitats or corridors for wildlife. However, mountain forests in particular face heightened vulnerability to the impacts of climate change, primarily as a result of temperature limitations and their increased susceptibility to warming (Albrich et al., 2020). Italian wooded areas are of paramount importance, contributing significantly to both economic and non-economic sectors and providing multipurpose services in terms of production (timber, firewood, and their derivatives), protection and recreation, among others. Indeed, the concept of Forest Ecosystem Services (FES), intended as the direct and indirect contributions to human wellbeing by forest ecosystems, was developed and well-described by the Millennium Ecosystem Assessment in 2010 (Millenium Ecosystem Assessment, 2010). The multifaceted contributions of Italian wooded areas underscore their significance as vital assets, crucial for sustaining both the economy and the overall wellbeing of society, providing industries and individuals with raw materials and sustainable energy sources—if managed according to modern sustainable forestry criteria (Buonincontri et al., 2023; Testolin et al., 2023). Similarly, non-wood forest products (i.e., mushrooms, chestnuts, truffles, seeds, etc.) represent additional sources of value as they supply human nutrition, renewable materials, cultural and experiential services, creating job and income opportunities in rural areas (Weiss et al., 2020). Furthermore, forests provide invaluable contributions to soil creation and preservation, serving as essential regulators of the hydro-geological cycle, as well as influencing water availability and quality. Forests also support and enhance biodiversity, fostering the existence of a wide range of species and contributing to the overall richness and ecological balance of ecosystems. Lastly, the recreational benefits offered by forests have gained recognition as essential for human wellbeing (Anderson et al., 2023). Hence, Italian forests offer valuable aid to local economies, even in cases where wood may not be fully utilized or economically optimized as a resource. Finally, forests are the most efficient and cheapest means of carbon dioxide removal from the atmosphere.

The rapid changes in climate, as highlighted by the Intergovernmental Panel on Climate Change (IPCC, 2022), have raised significant concerns regarding the health and functionality of forest ecosystems. Over the past few decades, there has been a notable intensification of disturbance regimes, posing challenges to the provision of various FES, particularly in terms of biodiversity protection. Extreme events such as droughts and storms are becoming more frequent, prolonged, and intense, significantly impacting the resilience of Italian forests (i.e., the “Vaia” storm in the North-Eastern Alps, in November 2018). Over several decades, widespread reports of forest mortality and decay have emerged in both the Apennines and plain environments [e.g., leading to the deterioration of oak species (Conte et al., 2019)] as well as in Alpine regions [notably with Scots pine forests in the North-West (Vacchiano et al., 2012)]. Additionally, the climate crisis has undeniably contributed to an increase in the number, intensity, and relative risk of forest fires (Bacciu et al., 2012). These mounting pressures are driving substantial changes in the species composition of historic Italian forest stands (Di Pasquale et al., 2020; Pecchi et al., 2020; Sferlazza et al., 2023). The rapid pace of these complex transformations surpasses the potential for evolutionary adaptation (Lindner et al., 2010; Trumbore et al., 2015) and migration processes. Fully understanding these dynamics is crucial to try defining hypotheses about future scenarios.

Species Distribution Models (SDMs)—also referred to as Correlative Species Distribution Models, bioclimatic envelope models, correlative ecological niche models, or habitat suitability models—are computational models used to examine the relationships and equilibrium between the geographic distribution of species or species groups and a set of environmental variables (Guisan and Thuiller, 2005; Austin and Niel, 2011; Noce et al., 2017, 2019). SDMs, in conjunction with Geographic Information Systems (GIS) tools, offer promising approaches for mapping and predicting the potential range expansion of endemic or invasive species in both historical and future contexts, respectively (Franklin, 2010; Ali et al., 2021; Sofi et al., 2022). By integrating species presence-location data with geospatial information, SDMs rely on various algorithmic approaches, including machine learning and regression-based methods, among others, to develop non-linear and discontinuous relationships between species and their environmental conditions (Heumann et al., 2013). These models provide valuable insights into the potential distribution ranges of species under different environmental scenarios, aiding in the assessment of species' responses to changing environmental conditions. Through the utilization of SDMs, a better understanding of how species distributions may shift in response to factors such as climate change, land use change, or other environmental drivers can be obtained (Salinas-Ramos et al., 2021; Jamwal et al., 2022). Such information is crucial for effective conservation and management strategies, as it allows for proactive measures to be implemented to mitigate potential negative impacts or to promote the preservation of endangered species.

This study aims to project the potential suitability of selected target forest species in Italy into the medium-term future, considering the foreseen environmental changes associated with the ongoing climate crisis, especially in mountain areas. We seek to assess the potential shifts in the distribution areas of these species under two distinct socio-economic forcing scenarios: a moderate scenario (RCP 4.5, Thomson et al., 2011) and a high-emission scenario (RCP 8.5 Riahi et al., 2011). By examining the variations in potential distribution areas, both in terms of areas lost and gained, we aim to provide insights into the potential impacts of different future trajectories on these forest species. This analysis will contribute to a better understanding of the potential consequences of the climate crisis on Italian mountain forest ecosystems and will assist in developing informed strategies for their conservation and management.

2. Materials and methods

For comprehensive methodological reporting, we adhered to the ODMAP (Overview, Data, Model, Assessment, Prediction) protocol v1.0 (for detail see Supplementary Table 1), as proposed by Feng et al. (2019) and Zurell et al. (2020). Further details can be found in the Supplementary material. All analyses were conducted in ESRI ArcGIS Pro 3.1.1, ESRI ArcMap 10.8.2 and SAGA-GIS 7.8.2. NetCDf files were processed with the Climate Data Operators 2.0.6 collection.

2.1. Study area and species occurrence data

To encompass the entire Italian territory, the study area was delimited by national administrative boundaries. Presence data for forest species were obtained from the second National Inventory of Forests and forest Carbon Pools - INFC 2005 (https://www.inventarioforestale.org/it/). This inventory, freely available for research purposes, provides presence data (Figure 1A) for both native and allochthonous forest species across Italy. The data was collected using a random sampling method based on a 1 km x 1 km grid system aligned with meridians and parallels (Gasparini and Tabacchi, 2011). From this dataset, we derived a subset consisting of the 20 most representative species in terms of coverage across the entire territory, including 13 broad-leaved and seven needle-leaved species (Table 1). The characteristics of the survey that generated the INFC 2005 inventory led us to classify the occurrence dataset used for the SDM approach as a presence-only dataset. At the time of our analysis, the most recent data from the third National Inventory of Forests and forest Carbon Pools—INFC 2015 (published in a geospatial compatible format in October 2022) was not yet available, thus necessitating the use of the INFC 2005 dataset. The final rarefied occurrence dataset was obtained removing spatially autocorrelated occurrence points with a 20 km distance.

FIGURE 1

Figure 1. (A) Locations of INFC2005 samples. (B) Sections considered from the “Italian Ecoregion Map” dataset (Blasi et al., 2014, 2018).

TABLE 1

Table 1. Species selected for the analyses, with their Latin and Common names.

2.2. Environmental predictors

To ensure optimal modeling performance, we incorporated Very High Resolution (VHR) climate data into our analyses. The VHR-REA_IT dataset (Raffa et al., 2021), with a resolution of 2.2 km, covers the entirety of the Italian territory. This dataset was obtained by downscaling the ERA5 reanalysis dataset, which has a native resolution of ‘31 km (Hersbach et al., 2018), to a resolution of around 2.2 km for the reference period 1981–2020, and using the Regional Climate Model (RCM) COSMO-CLM (Rockel et al., 2008).

From the VHR-REA_IT dataset, we selected four variables (maximum, minimum, mean temperature, and precipitation) at a native temporal resolution of 1 hour. These variables were then converted into monthly mean values for the period 1991-2020, referred to as the “historical” period. Subsequently, the Climate Tools Library in SAGA-GIS 7.8.2 (https://saga-gis.sourceforge.io/saga_tool_doc/7.7.0/climate_tools.html) was employed to process these variables and derive a series of 19 bioclimatic indicators following the definitions provided by Worldclim (https://www.worldclim.org/data/bioclim.html).

For the analysis of altitude and slope, the Digital Elevation Model (DEM) over Europe (EU-DEM v1.1) was used, as available from the Copernicus Land Monitoring Service (https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1). This DEM combines data from SRTM and ASTER GDEM sources, is at 25 m resolution and covers the EEA39 countries and provides information on altitude within the study area, allowing to derive other topographic derivatives through GIS tools. In this case, the slope was calculated. The environmental predictors considered for the SDM analyses comprise the 19 bioclimatic indicators and the 2 topographic variables.

2.3. Model fitting and tuning

To tune the modeling process, we used the SDMtoolbox package v2.5 (Brown et al., 2017) based on MaxEnt v3.4.3 (Merow et al., 2013; Phillips et al., 2017) (http://biodiversityinformatics.amnh.org/open_source/maxent) to develop spatial models with a logistic output of historical suitability. The presence data for the 20 target species, obtained from INFC2005, were used for model training. The MaxEnt machine learning algorithm was selected due to its several advantages over other algorithms, particularly its requirement of presence-only data (Chiang and Valdez, 2019).

During the tuning phase, various settings were explored to train the MaxEnt models, including the number of predictors, background data selection, model complexity, and threshold selection. To address multicollinearity among the predictors, a preliminary analysis was conducted as indicated by Dormann et al. (2013). Additionally, to prevent overfitting (Elith et al., 2010), we performed correlation analyses with pairwise Pearson excluding predictors higher than defined thresholds (0.7, 0.8, 0.9).

For the background data, different selection types (Minimum Convex Polygon, Buffer Distance from Observation Points) and selection distances ranging from 20 to 500 km were considered. As it was not possible to calibrate the model on independent data as suggested by Araujo et al. (2005), in each run, the presence data were divided into three groups to train both spatially segregated and non-spatially segregated models.

In the MaxEnt settings, “logistic” was set as the output format, the replicated run type was selected as “crossvalidate”, and the random test percentage was set to 20. The Replicates number was set to 5. Five feature classes were included: linear, quadratic, product, hinge, and threshold. Moreover, a combination of regularization multipliers (0.2, 0.5, 1, 1.5, 2, 5, 10) was employed to fine-tune the models. Response curves were generated to analyze the relationships between predictor variables and habitat suitability for the target species. Additionally, we assessed the importance of the predictor variables through jackknifing (Baldwin, 2009) with minimum occurrence points set to 15. The best model was selected considering AUC values and then the Omission Error Rate (OER). Additionally, the True Skill Statistic (TSS) (Allouche et al., 2006) was assessed for the best models. The minimum number of occurrence points to model distribution of species was set to 5.

We generated both continuous and binary outputs to assess the habitat suitability, or probability of occurrence, for each species. However, this study will focus solely on the discussion of the continuous outputs, whereas the binary outputs (presence-no presence) will not be addressed. To convert the continuous data into binary format, we employed two threshold methods: the 10th percentile training presence (PTP) and the maximum test sensitivity and specificity logistic (MTSS).

Subsequently, we employed a model selection approach (Zurell et al., 2020) to identify the best model. This approach involves comparing different model structures and settings to choose a single optimal model or a set of best-performing models. The selection of the best model is based on the need to enhance prediction accuracy by reducing the variance of predicted values or to facilitate interpretation (Hastie et al., 2009). In our study, the best model from MaxEnt was chosen based on a combination of evaluation measures, including omission error, area under the receiver operating characteristic curve (AUC), and prediction rate.

2.4. Future projections

Following the tuning phase and selection of the best model, we proceeded to use it to generate maps depicting the future land suitability for the target species for the future time period of 2021–2050. To obtain future projections, we recalculated the 19 bioclimatic indicators using the VHR-PRO_IT dataset (Raffa et al., 2023). This dataset is a downscaled version of the COSMO-CLM simulation over Italy, previously produced at 8 km (Bucchignani et al., 2016; Zollo et al., 2016) resolution and driven by the CMCC-CM General Circulation Model (Scoccimarro et al., 2011).

The VHR-PRO_IT dataset covers the time window of 1981-2070, for 1981-2005 under the historical greenhouse gas forcing, and for 2006-2070 under the Representative Concentration Pathways (RCPs) 4.5 and 8.5 (also “scenarios” hereafter). RCP 4.5 is a stabilization scenario where total radiative forcing is stabilized, shortly after 2100, to 4.5 Wm^-² (650 ppm CO₂-equivalent) by employing technologies and strategies to reduce GHG emissions, whereas RCP 8.5 is a business as usual scenario and it is characterized by increasing GHG emissions and high GHG concentration levels, leading to 8.5 Wm^-² in 2100 (1,370 ppmv CO₂-equivalent). In our calculations, we specifically considered the time span of 2021–2050, aligning with the guidelines provided by the Intergovernmental Panel on Climate Change (IPCC) in their 6th Assessment Report (IPCC, 2022) for evaluating climate change. Additionally, we adhered to the recommendations of the World Meteorological Organization (https://public.wmo.int/en/media/news/updated-30-year-reference-period-reflects-changing-climate) and the National Oceanic and Atmospheric Administration (NOAA, Bates et al., 2016) for statistical analyses of climate data, which recommend using a 30-year period as a representative measure of climate conditions in a given area.

The VHR-PRO_IT dataset is a scenario-driven simulation, therefore a correction of model bias was required to ensure comparability with the historical dataset VHR-REA_IT. We applied a constant anomaly correction based on the difference (for temperatures) or ratio (for precipitation) between VHR-REA_IT and VHR-PRO_IT during the overlapping period of 1991–2020, following the approach outlined in Maraun and Widmann (2018). This corrected version of the VHR-PRO_IT dataset was denoted as VHR-PRO_IT-ac and was calculated according to the following equations:

\begin{array}{l} T x^{VHR-PRO_IT-ac} (f u t) = T x^{VHR-PRO_IT} (f u t) \\ + (T x^{VHR-REA_IT} (o v e r) - T x^{VHR-PRO_IT} (o v e r)) & (1) \end{array}

\begin{array}{l} \begin{array}{l} P x^{VHR-PRO_IT-ac} (f u t) = P x^{VHR-PRO_IT} (f u t) * \frac{P x^{VHR-REA_IT} (o v e r)}{P x^{VHR-PRO_IT} (o v e r)} \\ where \frac{P x^{VHR-REA_IT} (o v e r)}{P x^{VHR-PRO_IT} (o v e r)} < 4 \end{array} & (2) \end{array}

and

\begin{array}{l} \begin{array}{l} P x^{VHR-PRO_IT-ac} (f u t) = P x^{VHR-PRO_IT} (f u t) * 4 \\ where \frac{P x^{VHR-REA_IT} (o v e r)}{P x^{VHR-PRO_IT} (o v e r)} \geq 4 \end{array} & (3) \end{array}

Where fut is the 2021–2050 time period; over is the 1991–2020 overlapping period; T x are temperatures (mean, max, min); P is precipitation. The threshold for the correction factor of precipitation was set to 4, simplifying the approach in Sperna Weiland et al. (2010) and considering analyses previously done for the domain under study.

To address the potential challenges arising from differences in environmental conditions between the historical and future periods, we employed the clamping option in MaxEnt during the future projection phase (Phillips et al., 2006; Radosavljevic and Anderson, 2014). This approach helps mitigate the effects of environmental discrepancies by constraining the model's response to values within the range observed during the calibration phase. By setting the clamping option, we aimed to enhance the reliability and accuracy of future projections, enabling better comparisons and interpretation of the results.

2.5. Elevation analyses

At first, the EU-DEM dataset was classified into 150-meter bands. Afterward, zonal statistics were conducted on the three suitability datasets within five distinct and homogeneous mountainous biogeographical regions defined by the “Italian Ecoregion Map” (Blasi et al., 2014, 2018). Specifically, we focused on five sections representing two regions in the Alps (Western and Central-Eastern) and three regions in the Apennines (Northern-Northwestern, Central, and Southern) (refer to Figure 1B for a visual representation of the selected sections). This division allowed for a more focused analysis of the suitability data within these specific mountainous regions, taking into account their unique characteristics and ecological dynamics.

3. Results

3.1. Model performance

The tuning phase, aimed at optimizing the model performance based on AUC values, resulted in the selection of specific MaxEnt settings. These settings include the exclusion of highly correlated predictors with a threshold of 0.9, background selection using the Minimum Convex Polygon method with a distance of 500 km, 5 replicates, and the consideration of Linear, Quadratic, Hinge, and Threshold feature classes with a regularization multiplier of 0.5. These settings were carefully chosen to ensure that a consistent and effective set of parameters maximized the performance of all 20 species simultaneously. The AUC results for these settings can be found in Table 1, providing a comprehensive overview of the model's discriminatory power for each species. We found higher AUC values for Swiss stone pine (0.963), Silver fir (0.928), and Maritime pine (0.922), lower values for Downy oak (0.771), Manna ash (0.804) and Aleppo pine (0.810). The TSS results can be found in Supplementary Table 2.

3.2. Elevation analyses

The generated suitability maps for the historical (4.5 and 8.5 scenarios) and future periods have a spatial resolution of approximately 2.2 km. Figure 2 displays as an example the predicted suitability specifically for the European beech under RCP 8.5 scenario.

FIGURE 2

Figure 2. Example of raw data, representing the output of our modeling procedure. Map of expected suitability for European beech 2021–2050, under the RCP 8.5 scenario. Visualization threshold 0.2.

3.2.1. Western Alps

This section encompasses a total area of 1,794,031 hectares and is bordered by the Maritime Alps to the South−West and Lake Maggiore to the Eastern boundary. The altitude within this section ranges from 26 to 4,790 meters. Notable changes in suitability are observed for various species within this area as shown in Table 2. The Silver fir is projected to experience a significant loss in suitability, ranging from −44% under the RCP 4.5 scenario to −20% under the RCP 8.5 scenario. Similarly, the Common hazel, European beech, and Sessile oak are also expected to undergo substantial suitability losses (−27% to −42% for Common hazel, −27% to −9% for European beech, and −12% to −17% for Sessile oak). A moderate loss in suitability is anticipated for the Chestnut (−8% to −7%). On the other hand, the Field maple, European larch, and Turkey oak are expected to experience significant gains in suitability (+26% to +32% for Field maple, +33% to +37% for European larch, and +15% to +38% for Turkey oak). A moderate gain in suitability is projected for the Pedunculate oak (+6% to +10%), Norway spruce (+18% to +7%), and European hornbeam (+3% to +17%). Divergent projections are observed between the two scenarios for the Manna ash (−6% to +29%), Hop hornbeam (−4% to +53%), Swiss stone pine (+1% to −18%), Maritime pine (−29% to +5%), and Downy oak (−24% to +19%).

TABLE 2

Table 2. Results for Section 1—Western Alps.

Regarding altitudinal shifts, noteworthy observations include the upward shifting of the maximum suitability range for the European hornbeam and Turkey oak, with an increase from 2 bands (300 m) under the RCP 4.5 scenario to 3 bands (450 m) under the RCP 8.5 scenario. The Chestnut (Figure 3A) is expected to experience a reduction in suitability at lower altitudes in both scenarios, with an upward shift of 2 bands (300 m) under RCP 8.5. Similar projections apply to the European beech, but only under the RCP 8.5 scenario, where a gain in suitability at higher altitudes with an upward shift of 3 bands (450 m) is expected. The Manna ash exhibits no changes in suitability under the RCP 4.5 scenario, but a substantial gain is projected above 600 meters above sea level (a.s.l.) under the RCP 8.5 scenario. The European larch (Figure 3B) and Hop hornbeam display similar patterns with a strong gain in suitability across all higher altitude bands. For the Norway spruce (Figure 3C), a gain in suitability is expected at altitudes above 1200 m, whereas no changes are anticipated below this threshold. The optimal range for the Swiss stone pine appears to shift upwards by approximately 450 m (3 bands). Unclear or divergent signals are observed for the other species. Supplementary Figure S1 contains all graphs not included in the main text.

FIGURE 3

Figure 3. Suitability in altitudinal bands for the Western Alps section. In gray, historical (1981–2020); light blue, future RCP 4.5 (2021–2050); orange, future RCP 8.5 (2021-2050). (A) Chestnut (B) European larch (C) Norway spruce.

3.2.2. Central and Eastern Alps

This section covers a total area of 3,656,143 hectares and extends from the Eastern shore of Lake Maggiore to the Julian Alps. The altitude within this section ranges from 25 to 3,950 meters above sea level. Noteworthy findings (Table 3) include a significant reduction in suitability for the Silver fir (-43% to -30%), Common hazel (-46% to -57%), and Sessile oak (-23% to -31%). A moderate reduction in suitability is projected for the Chestnut (-9% to -11%), European beech (-24% to -5%), and Swiss stone pine (0 to -13%), whereas gains are expected for the European larch (+39% to +40%) and Norway spruce (+21% to +7%). Divergent projections are observed between the two scenarios for the Field maple (-20% to +11%), European hornbeam (-9% to +23%), Scots pine (+6% to -8%), Manna ash (-12% to +42%), and Pedunculate oak (-34% to +5%).

TABLE 3

Table 3. Results for Section 2—Central and Eastern Alps.

Regarding altitudinal shifts, significant reductions in suitability are predicted across the entire altitudinal range for the Silver fir (Figure 4A), Common hazel, and Sessile oak. Minor reductions in suitability are expected for the Chestnut, particularly at altitudes below 1,000 meters. For the Field maple, a shift of 2 bands (300 meters) upwards is anticipated, particularly under the RCP 8.5 scenario. Similar shifts, but under both scenarios, are projected for the two Hornbeams and the Manna ash. The European beech (Figure 4B) shows an upward shift of 1 band (150 meters). Strong gains in suitability, even above the current tree line, are expected for the European larch and Norway spruce, particularly above the 1,350–1,500 meter band. Projections for the Stone and Scots pines (Figure 4C) indicate no major changes, whereas the Downy oak shows divergent results. Supplementary Figure S2 contains all graphs not included in the main text.

FIGURE 4

Figure 4. Suitability in altitudinal bands for the Central and Eastern Alps section. In gray, historical (1981-2020); light blue, future RCP 4.5 (2021–2050); orange, future RCP 8.5 (2021–2050). (A) Silver fir (B) European beech (C) Scots pine.

3.2.3. Northern and Northwestern Apennines

This section covers an area of 3,880,014 hectares and extends from the Ligurian Apennines in the North to the Tuscan−Romagna Apennines in the South, including the hills known as the “Colline Metallifere”. The altitude within this section ranges from 10 to 2,142 meters above sea level. Our projections (Table 4) indicate a reduction in suitability for all species, except for the Pedunculate and Sessile oaks, under the RCP 4.5 scenario. Specifically, there is a strong reduction in suitability for the Silver fir (−43% to −33%), Common hazel (−5% to −42%), European hornbeam (−26% to −56%), Hop hornbeam (−38% to −14%), Chestnut (−23% to −26%), European beech (−32% to −22%), Manna ash (−21% to −28%), Norway spruce (−21% to −38%), Aleppo pine (−59% to −13%), Maritime pine (−48% to −31%), and Downy oak (−26% to −22%), as well as a moderate reduction for the Field maple (−18% to −15%) and Turkey oak (−11% to −12%). Divergent results are observed for the Sessile oak (+5% to −13%) and Pedunculate oak (+5% to −13%).

TABLE 4

Table 4. Results for Section 3—Northern and Northwestern Apennines.

Regarding altitudinal shifts, a reduction in suitability for the Silver fir is expected at lower altitudes (below 1,500 m a.s.l.). The Field maple, European beech (Figure 5A), and Maritime pine (Figure 5B) are projected to shift upward by one band (150 meters), whereas the European hornbeam, Scots pine, and Sessile oak are anticipated to shift up by two bands (300 meters). Similar projections are expected for the Chestnut and Turkey oak (Figure 5C), with a reduction in suitability up to 800 m a.s.l. and a gain in suitability above. There is a significant reduction in suitability for the Holm oak (across the entire altitudinal range) and Downy oak (up to 600 m a.s.l.), whereas divergent results are observed for the Hop hornbeam. Supplementary Figure S3 contains all graphs not included in the main text.

FIGURE 5

Figure 5. Suitability in altitudinal bands for the Northern and Northwestern Apennines section. In gray, historical (1981-2020); light blue, future RCP 4.5 (2021–2050); orange, future RCP 8.5 (2021–2050). (A) European beech (B) Maritime pine (C) Turkey oak.

3.2.4. Central Apennines

This particular section covers an area of 2,639,776 hectares and extends from the Umbria-Marche Apennines in the North to the Mainarde and Maiella mountains in the South. The altitude within this section ranges from 0 to 2,850 meters above sea level. As reported in Table 5, we expect a strong reduction in suitability for the Silver fir (−30% to −24%), Aleppo Pine (−60% to −18%), and the Holm (−49% to −25%) and Downy (−40% to −16%) oaks. There is a moderate reduction for the Chestnut (−6% to −4%), European beech (−14% to −8%), Manna ash (−14% to −16%), European (−15%) and Hop (−20% to 0%) hornbeam, Field maple (−6% to −12%), and Maritime pine (−13% to 0%), whereas a gain in suitability is expected for the Turkey oak (+8% to +1%) and Pedunculate oak (+26% to +17%). Divergent projections are observed between the two scenarios for the Common hazel (+16% to −19%) and Sessile oak (+15% to −13%).

TABLE 5

Table 5. Results for Section 4—Central Apennines.

Regarding altitudinal shifts, a reduction in suitability is projected across the entire altitudinal range for the Silver fir (particularly under the RCP 4.5 scenario), Aleppo pine, and Holm and Downy oaks (Figure 6A). Our projections indicate an upward shift (two bands or 300 meters) for the Maritime pine. Similar patterns are observed for several species, with a reduction in suitability up to specific altitudes and a gain in suitability above. For example, the Field maple shows a reduction up to 600 m a.s.l. and a gain above, the Chestnut (Figure 6B) shows a reduction up to 900 m a.s.l. and a gain above the tree line, the European beech shows a reduction up to 1,500 m a.s.l. and a gain up to the tree line, and the Turkey oak (Figure 6C) shows a reduction up to 600 m a.s.l. and a gain up to the tree line. We also anticipate an increase in suitability for the Pedunculate oak in the 500–1,000 m a.s.l. range. In conclusion, projections for the Manna ash, Hop hornbeam, Common hazel, and Sessile oak show divergent results. Supplementary Figure S4 contains all graphs not included in the main text.

FIGURE 6

Figure 6. Suitability in altitudinal bands for the Central Apennines section. In gray, historical (1981-2020); light blue, future RCP 4.5 (2021-2050); orange, future RCP 8.5 (2021-2050). (A) Downy oak (B) Chestnut (C) Turkey oak.

3.2.5. Southern Apennines

Encompassing a total area of 1,943,464 hectares, the Southern Apennines section stretches from the Matese massif in the North to Pollino in the South. Altitude ranges from 32 to 2,250 meters above sea level. Projections (Table 6) indicate a significant decrease in suitability for the Silver fir (−39% to −39%), Aleppo pine (−56% to −26%), and Holm (−11% to −31%) and Downy (−26% to −20%) oaks. The Chestnut (−3% to −12%), European beech (−3% to −12%), Manna ash (−12% to −8%), and European (−13% to −20%) and Hop (−12% to −8%) hornbeams are expected to experience a moderate reduction in suitability. On the other hand, the Turkey (+13% to +5%) and Pedunculate (+35% to +36%) oaks, as well as the Maritime pine (+46% to +15%), are projected to gain suitability. Divergent outcomes are observed between the two scenarios for the Common hazel (0% to −32%), Field maple (+5% to −1%), and Sessile oak (+34% to −23%). Regarding altitudinal shifts, reductions in suitability are expected across the entire altitudinal range for the Silver fir (Figure 7A), European hornbeam, Common hazel (especially under the RCP 8.5 scenario), Manna ash, Aleppo pine (under RCP 4.5), and Downy oak. The results indicate minor reductions in suitability for the European beech up to an altitude of 1,400 meters, followed by an increase above this altitude, primarily around the current tree line. Strong gains in suitability, accompanied by an upward shift, are expected for the Maritime pine (especially from 750 meters to the current tree line—Figure 7B) and Turkey oak (from 600 meters above sea level to the current tree line—Figure 7C). Additionally, an increase in suitability for the Pedunculate oak is projected within the 450–1,050-meter bands, similar to the previous section. Divergent results were obtained for the Hop hornbeam. Supplementary Figure S5 contains all graphs not included in the main text.

TABLE 6

Table 6. Results for Section 5—Southern Apennines.

FIGURE 7

Figure 7. Suitability in altitudinal bands for the Southern Apennines section. In gray, historical (1981–2020); light blue, future RCP 4.5 (2021–2050); orange, future RCP 8.5 (2021–2050). (A) Silver fir (B) Maritime pine (C) Turkey oak.

4. Discussion

The findings of this study provide insights into the potential future dynamics of climate suitability for key native forest species in Italy.

Despite exhibiting strong performance values during the training phase, the modeling results revealed variations in performance among species. These discrepancies can be attributed to the ecological characteristics of the species. The performance values align well with the findings of Tsoar et al. (2007) and more notably, with McPherson and Jetz (2007), who observed that predictions tend to be more accurate for species with smaller range sizes and higher habitat specificity (as exemplified by the Swiss stone pine in our study) compared to more generalist species with broader ranges (such as the Downy oak or Manna ash).

Our focus was on five mountainous sections, including two in the Alps and three in the Apennines. Although with many overlapping points, future trajectories reveal diversified impacts among species and scenarios, with the RCP 4.5 scenario showing slightly worse overall outcomes. Most species are expected to experience a contraction in their altitudinal range of suitability, but some show a propensity to extend beyond the current tree line, as observed in previous studies (Cudlin et al., 2017; Beniston et al., 2018). Among the mountain sections considered, the Northern and North-Eastern Apennines exhibit the greatest and most widespread impacts on all species. Overall, it is difficult to unambiguously define successful or unsuccessful species, except for the Silver fir, which is projected to be highly vulnerable across all altitudinal bands and mountainous sections. Two species worth noting are the larch in the Alpine region and the Turkey oak in the Apennines, as they show potential gains and could play significant roles in maintaining wooded populations. However, our results regarding the European larch differ from previous studies, which found negative variations (Dyderski et al., 2018; Mamet et al., 2019; Pecchi et al., 2020), whereas our findings suggest the opposite. If confirmed, this species could become crucial for future planning activities, particularly at higher altitudes (1,300–1,500 m). Our results for the Turkey oak contradict some studies (Vitale et al., 2012; Pecchi et al., 2020) but align with others (Noce et al., 2017). The expected increase in suitability at very high altitudes for some major species, both in the Alps and, particularly, the Southern Apennines, implies a potential upward shift of the tree line, consistent with previous research (Harsch et al., 2009; Greenwood and Jump, 2014).

This may result in various consequences. Firstly, it could lead to a loss of diversity as specialized species with limited niche tolerance might disappear due to competition from more widely distributed species (Jump et al., 2012). Additionally, it may lead to a reduction in the number and available area of high mountain ecosystems, such as nival vegetation or alpine grasslands (Moiseev and Shiyatov, 2003; Garćıa-Romero et al., 2010), thereby altering the balance between various ecosystem services provided (Peng et al., 2009). For a comprehensive examination of the consequences of treeline shifts, please refer to the valuable work by Greenwood and Jump (2014).

The European beech, a keystone species in the Italian mountain environment, shows evident impacts as reported in Buonincontri et al. (2023). We observe an upward shift in its distribution within the Alpine arc and Northern Apennines, whereas good future suitability is expected at higher altitudes (above 1500 meters) in the Central and Southern Apennines. The Maritime pine emerges as a promising candidate for the future of the Southern Apennines, potentially expanding its presence and altitudinal distribution in association with other broad-leaved species. These findings align with a study in Spain (Barrio-Anta et al., 2020). However, the increased flammability of the Maritime pine (Núñez-Regueira et al., 1996) raises concerns, particularly considering the expected rise in fire risk in the region (Spano et al., 2020).

Regarding the predicted tree line upward shifting, it is plausible to attribute this phenomenon to the expected temperature increase in mountainous areas, which is known to impact the thermal limitations on the altitudinal distribution of species, including freezing tolerance and growth requirements (Körner, 2021). However, it is crucial to emphasize that our study focused solely on modeling climatic suitability, without considering other relevant factors that regulate the establishment of stable forest populations, such as soil availability and the potential intensification of high-altitude winds, including extreme wind events, which are expected to increase in frequency for Southern Europe (Outten and Sobolowski, 2021).

However, it is important to acknowledge certain limitations associated with the approach we have adopted and therefore the results obtained. The SDM approach is known to be subject to various assumptions and uncertainties (Guisan and Thuiller, 2005; Watling et al., 2014; Santini et al., 2021). Moreover, the assumption that the relationships between environmental variables and species presence observed in the historical period will remain consistent in the future introduces a considerable degree of uncertainty (Gavin et al., 2014). Furthermore, the accuracy and quality of the occurrence dataset represent a crucial factor that can influence the reliability of the results (Bloom et al., 2018). Although the INFC dataset is robust and based on a systematic sampling scheme across Italy, the privacy restrictions associated with providing only the coordinates of the southwestern corner of the 1 km grid introduce a certain level of uncertainty. However, given the coarser resolution of the climate data (2.2 km), this uncertainty is considered negligible, as demonstrated by Marchi and Ducci (2018). Another source of uncertainty stems from the temporal mismatch between the INFC 2005 survey, which took place over a couple of years starting in 2003, and the historical climatic data used to calibrate the model, which refers to the period 1991–2020. Despite this temporal discrepancy, it is important to consider that changes in forest species composition occur over longer timescales. The reliability of the results is also linked to the quality of the environmental variables used . In this regard, the VHR datasets employed in our study allow us to provide useful information at different scales, from local to regional and national. This highlights the importance of considering different emission scenarios to capture the range of possible outcomes and effectively plan for future conservation and management strategies. Lastly, it should be noted that the results obtained are highly dependent on the parameterization or configuration of the applied models (Hallgren et al., 2019). Adhering to the best practices defined by the ODMAP scheme has allowed us to address this aspect with a logical, transparent, and reproducible methodology (Fitzpatrick et al., 2021). The aforementioned aspects, among others (Jarnevich et al., 2015), highlight the significance of uncertainty in these types of studies. Consequently, all our results are presented as anomalies between the different simulations (historical and future).

The divergent projections observed between scenarios suggest varying impacts of climate change on suitability for the species under consideration, which can be seen as both a limitation and a strength of our study. These differences can be interpreted as the upper and lower bounds of the projected outcomes. Our results highlight the complex and dynamic nature of possible climate change impacts, emphasizing the need to consider multiple factors and scenarios when assessing species vulnerability and planning conservation actions. In conclusion, we expect significant and far-reaching impacts on mountain biodiversity, particularly in terms of forest population composition. The rapid pace of climate change in mountainous regions appears incompatible with the adaptive capacities and dynamics of arboreal plant organisms. This work also emphasizes the importance of using very high-resolution climate data, which is essential for formulating hypotheses about future forest dynamics and providing valuable information across different scales. Our findings have implications at the local, regional, and national levels and provide information that can improve future woodland management strategies. In further detail, this study can be useful in identifying priority locations for conservation, offering valuable guidance for multiple aspects of forest management and restoration. Firstly, providing insights that can assist in the selection of suitable species for future reforestation policies, considering their potential success in specific areas. Moreover, our results can aid in promoting certain species over others in silvicultural choices, which is crucial for optimizing ecosystem benefits and promoting biodiversity and resilience in managed forest stands (Testolin et al., 2023). Accelerating the onset of new species compositions can be particularly important in the context of ongoing environmental changes and the need to enhance forest ecosystem adaptability. Furthermore, our study has practical applications in guiding silvicultural decisions, including the number, frequency, and intensity of treatments. This aspect is essential for sustainable forest management, as it ensures proper stand development and growth dynamics. Our research also sheds light on areas where certain species are currently uncommon, offering options for potential successful species introductions and diversification. Conversely, we identify species that are already highly present in certain areas, suggesting the need to exclude them to avoid ecological imbalances and support biodiversity conservation efforts. Moreover, our findings suggest that processes aimed at making conditions more favorable for upward migration and tree line elevation on mountain ranges should be undertaken accompanied by soil protection actions and recurring stand health surveys and monitoring as described in Tomback et al. (2022). Simultaneously, it is essential to minimize anthropic disturbances, including the expansion of ski facilities or the construction of new slopes, particularly in proximity to the current upper limit of forest stands across all the mountain sections under consideration (Maliniemi and Virtanen, 2021). As demonstrated in previous studies (Jactel et al., 2017; González de Andrés, 2019), our findings support the transition toward silvicultural practices that favor mixed stands. This approach has been shown to promote ecosystem stability, increase resistance to disturbances, and improve overall forest health. In summary, our research has wide-ranging implications for conservation and forest management strategies. It provides valuable information to guide decisions related to species selection, silvicultural practices, and compositional arrangements in forest stands. By prioritizing these considerations, we can foster more resilient and diverse forests, enhancing their ecological value and ensuring their ability to adapt to changing environmental conditions.

Furthermore, they highlight once again the importance of considering different emission scenarios to encompass the full range of potential outcomes and effectively plan for future conservation and management strategies. It is important to acknowledge that this type of study is characterized by considerable uncertainty, and continued efforts are required to produce increasingly reliable datasets and forecasts. Understanding the climatic vulnerability of different species can assist in prioritizing conservation efforts and implementing targeted management strategies.

In the near future, our goal is to update our analyses using data from the latest national forest inventory, enlarging the set of species considered. Additionally, the inclusion of human-introduced, allochthonous, and invasive species could be a further step in enhancing our understanding of future forest dynamics.

Data availability statement

The datasets generated for this study are available in GIS format, primarily in raster form, from the corresponding author upon reasonable request. All datasets related to the RCP 8.5 scenario, including bioclimatic indicators and species suitability, are openly accessible in NetCDF format at https://dds.highlanderproject.eu/app/datasets/land-suitability-for-forests.

Author contributions

SN: conceptualization, formal analysis, and visualization. SN, CC, and MS: methodology and writing—review and editing. SN and CC: data curation and writing—original draft preparation. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially funded by the HIGHLANDER (HIGH performance computing to support smart LAND sERvices) project (Connecting European Facility—CEF—Telecommunications sector under agreement number INEA/CEF/ICT/A2018/1815462) and by the European Union-NextGenerationEU in scope of the National Biodiversity Future Center through Italian Ministry of Education, Universities and Research MIUR PNRR Mission 4.

Acknowledgments

The authors wish to thank Giulia Capotorti from Sapienza University (Rome, Italy) for providing the last version of the Italian Ecoregion Map dataset.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2023.1250651/full#supplementary-material

References

Agnoletti, M., Piras, F., Venturi, M., and Santoro, A. (2022). Cultural values and forest dynamics: The italian forests in the last 150 years. For. Ecol. Manage. 503, 119655. doi: 10.1016/j.foreco.2021.119655