Ecological Niche Modelling of Forest Tree Species in the Alpine Space: a Stacked SDM Approach at Regional Scale

UNIVERSITÀ DEGLI STUDI DI PADOVA Department of Land, Environment, Agriculture and Forestry Second Cycle Degree in Forestry and Environmental Science Ecological Niche Modelling of Forest Tree Species in the Alpine Space: a Stacked SDM Approach at Regional Scale Modellizzazione della nicchia ecologica delle specie arboree forestali in territorio alpino: un approccio stacked SDM a scala regionale Supervisor Prof. Francesco Pirotti Co-supervisors Dott. Michele Dalponte Dott. Davide Andreatta Dott. Harin Aiyanna Cheriyanda Raveendra Graduate candidate Damiano Oberosler Student ID 2090812 Academic year 2024 – 2025 2 3 TABLE OF CONTENTS Abstract ....................................................................................................................... 5 Riassunto ..................................................................................................................... 7 1 INTRODUCTION ............................................................................................. 9 1.1 What Are Ecological Niche Models? ...................................................................... 11 1.1.1 Ecological Niche Concept ............................................................................... 12 1.1.2 Ecological Niche Models vs. Species Distribution Models ............................ 14 1.2 A Brief History of Ecological Niche Modelling ...................................................... 15 1.3 Practical Applications of Ecological Niche Models ................................................ 18 1.4 Species Richness ...................................................................................................... 26 1.5 Knowledge Gap and Objectives ............................................................................... 27 2 MATERIALS AND METHODS ................................................................... 29 2.1 Study Area ............................................................................................................... 29 2.2 Methodological Framework: R Packages and Supporting Software ....................... 30 2.3 Species Occurrence and Environmental Predictor Datasets .................................... 32 2.3.1 Species Occurrence ......................................................................................... 32 2.3.2 Climatic Variables ........................................................................................... 35 2.3.3 Topographic Variables .................................................................................... 38 2.3.4 Soil Variables .................................................................................................. 38 2.3.5 Disturbance Variables ..................................................................................... 39 2.4 Model Framework .................................................................................................... 39 2.4.1 Ensemble Species Distribution Models ........................................................... 41 2.4.2 Stacked Species Distribution Models .............................................................. 45 2.4.3 Evaluation Metrics and Uncertainty Estimations ............................................ 46 2.4.4 Case studies ..................................................................................................... 47 3 RESULTS ......................................................................................................... 49 3.1 Projected Ecological Niches of 23 Forest Tree Species in Trentino ........................ 49 3.1.1 Evaluation metrics ........................................................................................... 61 3.2 Species Richness ...................................................................................................... 63 3.3 Case Studies ............................................................................................................. 69 3.3.1 Pinus cembra ................................................................................................... 70 3.3.2 Quercus pubescens .......................................................................................... 80 4 4 DISCUSSION .................................................................................................. 91 4.1 Interpretation of Distributional Changes .................................................................. 91 4.2 Species Richness Patterns and Community-Level Insights ...................................... 92 4.3 Environmental Variable Selection and Importance .................................................. 92 4.4 Performance and Uncertainty Assessment ............................................................... 94 4.5 Study Limitations and Future Developments ........................................................... 97 4.5.1 Latent Variables and Biotic Interactions .......................................................... 97 4.5.2 Occurrence Datasets and Human Influence ..................................................... 98 4.6 Management and Conservation Implications ......................................................... 100 5 CONCLUSIONS............................................................................................ 103 Supplementary Material ........................................................................................ 105 Bibliography ............................................................................................................ 105 Acknowledgments ................................................................................................... 116 5 ABSTRACT Climate change is deeply altering the structure and composition of forest ecosystems, particularly in mountain landscapes. Ecological Niche Modelling (ENM) provides a valuable tool to estimate the potential distribution of species under current and future climate scenarios. Although many studies have assessed tree-species distributions at continental or national extents with coarse spatial resolution, few have addressed regional, high-resolution modelling – despite its clear relevance for local forest planning and conservation. Key challenges for fine-scale modelling include limited field data availability and high computational and methodological complexity. This thesis implements Ecological Niche Modelling at 50-m resolution to estimate the potential distribution of 23 dominant forest tree species in the Autonomous Province of Trento (Italian Alps). The modelling was performed using the R package SSDM which combines an ensemble of statistical and machine-learning algorithms to produce individual species models and community-level outputs via Stacked Species Distribution Modelling (SSDM). Models were calibrated using downscaled, high-resolution climate projections (ECLIPS2.0), terrain metrics from a Digital Terrain Model (DTM), and detailed occurrence data from local forest inventories, remote sensing, and a citizen-science platform. Model performance was evaluated using metrics including Area Under the receiveroperating-characteristic Curve (AUC), Cohen’s Kappa, sensitivity, and specificity. Species richness maps were produced using the SSDM framework, providing a replicable protocol for modelling both individual species distributions and community composition. Additionally, future projections for the 23 species under two IPCC emission scenarios (RCP4.5 and RCP8.5) were generated to assess potential vegetation shifts across three time periods (2041–2060, 2061–2080, and 2081–2100). The analysis also included the assessment of feature importance and modelling uncertainty, which improves the interpretability and robustness of the projections. Species were classified according to projected changes in their distributional ranges. Taxa such as Fagus sylvatica L., Fraxinus ornus L. and Quercus ilex L. are predicted to exhibit range expansion (“gainers”), whereas Picea abies (L.) H.Karst., Pinus cembra L. and Fraxinus excelsior L. are anticipated to contract (“losers”). Other taxa, including Betula pendula Roth, Larix decidua Mill., and Populus tremula L., display non-linear temporal trend, with initial gains followed by subsequent losses of suitable areas. All the studied species will experience an upward elevation shift of their climate optimum. 6 Species richness is projected to remain largely stable under the RCP4.5 scenario, whereas a decline is expected under RCP8.5. Results can support adaptive forest planning and biodiversity conservation in the Alpine contexts. The approach adopted in this study contributes to bridging the gap between broad-scale modelling and local forest management needs. The work provides spatial insights into biodiversity patterns and could help to identify areas of high ecological significance or conservation concern. KEYWORDS: Niche Modelling; Stacked SDMs; Species richness; Occurrence mapping; Fine-scale modelling. 7 RIASSUNTO Il cambiamento climatico sta modificando sostanzialmente la struttura e la composizione degli ecosistemi forestali, con effetti particolarmente rilevanti nelle aree montane. La modellizzazione della nicchia ecologica (Ecological Niche Modelling, ENM) rappresenta uno strumento utile per stimare la distribuzione potenziale delle specie in scenari climatici presenti e futuri. Sebbene numerosi studi abbiano analizzato la distribuzione delle specie arboree su scale continentali o nazionali utilizzando dati a bassa risoluzione spaziale, poche ricerche hanno affrontato la modellizzazione regionale ad alta risoluzione, nonostante la sua evidente rilevanza per la pianificazione forestale locale e le politiche di conservazione. Le principali sfide per la modellizzazione ad alta risoluzione comprendono la limitata disponibilità di dati di campo e l’elevata complessità computazionale e metodologica. La presente tesi applica l’ENM a una risoluzione di 50 m per stimare la distribuzione potenziale di 23 specie arboree dominanti nella Provincia Autonoma di Trento (Alpi italiane). La modellizzazione è stata eseguita con il pacchetto R SSDM, che integra un insieme ensemble di algoritmi statistici e di machine learning per produrre modelli di specie individuali e output a livello di comunità mediante l’approccio di Stacked Species Distribution Modelling (SSDM). I modelli sono stati calibrati utilizzando proiezioni climatiche ottenute tramite downscaling ad alta risoluzione (ECLIPS2.0), metriche topografiche derivate dal Modello Digitale del Terreno (DTM) e dati di presenza dettagliati provenienti da inventari forestali locali, telerilevamento e una piattaforma di citizen science. La validazione dei modelli è stata effettuata mediante l’uso di diverse metriche di valutazione, tra cui l’area sotto la curva ROC (AUC), il coefficiente Kappa di Cohen, la sensibilità e la specificità. Sono state inoltre prodotte mappe di ricchezza specifica tramite l’approccio SSDM, fornendo un protocollo replicabile per la modellizzazione sia delle singole specie sia della composizione comunitaria. Le proiezioni future sono state generate per tre orizzonti temporali (2041–2060, 2061–2080, 2081–2100) e per due scenari di emissione dell’IPCC (RCP4.5 e RCP8.5). L’analisi include infine una valutazione dell’incertezza modellistica, volta a migliorare l’interpretabilità e la robustezza delle proiezioni. Le specie sono state classificate in base alle proiezioni di variazione dei rispettivi areali di distribuzione. Specie come Fagus sylvatica L., Fraxinus ornus L. e Quercus ilex L. sono 8 previste in espansione del loro areale («gainers»), mentre Picea abies (L.) H.Karst., Pinus cembra L. e Fraxinus excelsior L. sono attese in contrazione («losers»). Altre specie, tra cui Betula pendula Roth, Larix decidua Mill. e Populus tremula L., mostrano andamenti temporali non lineari, con guadagni iniziali seguiti da successive perdite dell’areale ottimale. La ricchezza specifica è prevista rimanere sostanzialmente stabile nello scenario RCP4.5, mentre ci si aspetta una diminuzione nello scenario RCP8.5. I risultati possono supportare la pianificazione forestale adattiva e la conservazione della biodiversità in contesti alpini. L’approccio adottato contribuisce a colmare il divario tra modellizzazioni su larga scala e le esigenze gestionali locali, offrendo informazioni spaziali dettagliate sui pattern di biodiversità e strumenti utili per l’individuazione di aree di elevato valore ecologico o di interesse conservazionistico. PAROLE CHIAVE: Modellizzazione della nicchia ecologica (ENM); Modelli aggregati di distribuzione delle specie (SSDM); ricchezza specifica; mappatura delle presenze; modellizzazione ad alta risoluzione. 9 1 INTRODUCTION Terrestrial ecosystems currently face an unprecedented challenge in both scale and pace: climate change (Parmesan et al., 2022). Anthropogenic greenhouse-gas emissions have elevated atmospheric CO₂ concentrations above 420 ppm – levels unmatched during the past 800,000 years (Lüthi et al., 2008) and by some estimates not seen for up to 14 million years (Hönisch et al., 2023). Global mean surface temperature has risen by approximately 0.9–1.2 °C relative to the 1850–1900 baseline (Eyring et al., 2021). Some regions experience larger changes; for instance, the Alpine region has warmed by about 1.8 °C since 1880 (Kotlarski et al., 2023). Anthropogenic forcing has also increased the magnitude, frequency, and spatiotemporal variability of extreme weather events (Seneviratne et al., 2021). Environmental changes have consistently shaped Earth's history, influencing the evolution and distribution of both plant and animal species (Huntley, 1991; Pitelka & Group, 1997). Over millions of years, habitats, including forest ecosystems, have maintained a dynamic equilibrium with environmental conditions, while species individually and gradually tracked their ecological niches in response to glacial and interglacial cycles (Davis & Shaw, 2001). Migration has thus emerged as a natural strategy for survival, alongside adaptation and genetic evolution. However, such processes typically occur over longer temporal scales and can entail significant costs and constraints, such as the loss of genetic variability or collateral effects on other adaptive traits (De Meester et al., 2018). Nevertheless, the current pace of climate transformation represents a historical anomaly. According to Hansen et al. (2023), the rate of global warming is exceeding 0.27° per decade, 10 times faster than an average postglacial era. This acceleration is pushing species beyond the limits of their spatial and temporal response capacities, contributing to the collapse of ecosystems that were historically considered stable (Canadell & Jackson, 2021). Forests, in particular, are proving to be highly vulnerable to these changes, primarily due to the longevity and relatively slow response times of tree species (Forzieri et al., 2021; Lindner et al., 2010; Vacek et al., 2023). Under stable conditions, forest ecosystems are known for their resilience and plasticity, enabling them to respond to environmental variations through local acclimation processes or by migrating toward more climatically favourable regions. However, the effectiveness of these strategies is increasingly compromised by the rapid advance of climate change, to the extent that, in the absence of a viable response, the risk of local extinction becomes tangible (Aitken et al., 2008; Loreto 16 in Beery et al., 2021) and Johnston (1924) studied the spread of invasive cactus species in Australia (cited in Guisan & Thuiller, 2005). By the 1950s, studies demonstrated that species respond individually to environmental changes, suggesting that modelling individual species is more appropriate than modelling entire communities (Elith & Leathwick, 2009). In 1963, Hintikka was among the first researchers to evaluate the distribution of various European species based on climatic variables (cited in Pearson & Dawson, 2003). Between the 1970s and 1980s, the first computer-based studies emerged (Guisan & Thuiller, 2005). Subsequently, more robust statistical models were developed, such as the Generalized Linear Models (GLMs), marking the beginning of regressionbased SDMs. These models introduced a new level of quantitative sophistication and realism, and are still widely used today (Elith & Leathwick, 2009). The advent and evolution of machine learning models has led to significant advancements in species distribution modelling. Machine learning models are adaptive algorithms capable of learning and self-modifying based on performance assessment criteria (Gobeyn et al., 2019; Figure 1.2), Figure 1.2. Chronological overview of the development of key statistical and machine learning methods used in Species Distribution Modelling. Source: Gobeyn et al., 2019 17 A significant advancement occurred in the late 1990s with the release of various highresolution global climate databases, including WorldClim – a global dataset of gridded environmental variables – which substantially enhanced the accuracy of species distribution predictions (Pecchi et al., 2019). Simultaneously, during the 2000s, the technology for collecting and digitizing geographic, topographic, and climatic data improved significantly. The development of Geographic Information Systems (GIS) enabled more effective handling of both presence data and environmental variables (Elith & Leathwick, 2009; Kearney & Porter, 2009). Driven by the exponential growth in computational power, current research is increasingly focused on the development of deep learning tools capable of providing more accurate information on species distributions through the use of machine learning algorithms (Bonannella et al., 2022). These algorithms can process large volumes of increasingly accessible field and environmental datasets and are better equipped to capture the interacting and non-linear relationships between species and their environments (Beery et al., 2021). Further advancements are also taking place in the field of meteorology, with the development of increasingly high-resolution maps for predicting climate evolution. Despite these advancements, the majority of studies have been conducted at broad spatial scales and low resolution (e.g. 1 to 50 km), which are not well-suited for practical forest management (Anselmetto et al., 2025). Maps with resolutions finer than 100 meters are more useful for policymakers, even at the expense of reduced prediction accuracy (Bonannella et al., 2022). Although data resolution has recently improved, the challenge remains to generate operational data at the local scale. This is largely due to the lack of high-resolution climate projections and the scarcity of accurate, well-georeferenced presence datasets, particularly in more remote areas. Limited progress has also been made at the regional scale regarding species richness. In managed forests, which constitute the majority of Europe’s forested areas, biological limitations such as competition and dispersal tend to be less significant, as these factors are largely governed by human intervention, whereby abiotic factors gain more importance in driving the species spatial distribution (Thurm et al., 2018). In this context, Ecological Niche Models, which are based on predicting the distribution solely on abiotic factors, can become particularly effective tools for silvicultural planning. 18 1.3 Practical Applications of Ecological Niche Models The forestry literature on Environmental Niche Modelling (or Species Distribution Modelling) is extensive. Numerous studies have been conducted on modelling the niches of species of major economic, ecological, social, and environmental importance across Europe. However, this non-exhaustive review of the available literature reveals a clear gap: a lack of local-scale studies that could provide concrete support for silvicultural planning. Studies conducted over large areas at continental or national scales are certainly valuable for understanding macro-trends in species distribution dynamics, but they are not suitable for forest management, which requires decision-making based on sub-compartment level information, typically at a spatial resolution of tens to hundreds of meters (Anselmetto et al., 2025). Furthermore, the use of local presence datasets, as opposed to more generalized or global ones, is considered to enhance the accuracy of predictions, as these datasets are less prone to collection biases and errors, but are subject to the niche truncation bias (Anselmetto et al., 2025). A brief review of the literature on species richness modelling is also provided. This synthesis highlights that few studies have applied Stacked species distribution models (SSDMs) to investigate these community-level parameters specifically within the forest sector. Starting with the single-species models, examples of continental-scale applications include Bonannella et al. (2022) for Europe, Henderson et al. (2014) for North America and Pérez Chaves et al. (2018) for South America. In these studies, the authors modelled potential current species distributions (at a high resolution of 30 meters) without projecting them in the future. The study by Bonannella et al. (2022), focusing on modelling both the Actual Natural Vegetation (ANV) and Potential Natural Vegetation (PNV) distribution of 16 of the most widespread forest species in Europe, currently represents one of the most comprehensive contributions in the field of SDM. The presence data used were sourced from heterogeneous databases including GBIF and national forest inventories. They cover the entire European continent and were subjected to spatial and qualitative filtering. Furthermore, the spatial resolution of the model is remarkably high (30 m), making this study particularly advanced from a methodological perspective. However, despite the overall high quality of the modelling framework, some limitations remain regarding its applicability at the local scale, especially in mountainous areas with complex morphology 19 such as the Autonomous Province of Trento (PAT). Specifically, for species with a limited or marginal distribution in the Alpine context, such as Quercus ilex L. (holm oak), the model demonstrates reduced predictive capacity. As shown in Figure 1.3, this species is not accurately modelled in areas where it is known to be stably present, such as the Val dei Laghi. This discrepancy is particularly evident in the case of potential distribution (Figure 1.3, A), suggesting that the model’s effectiveness may be lower in biogeographically peripheral areas relative to the species’ core range. Figure 1.4 illustrates that the distribution of holm oak presence points used to train the model is heavily concentrated in the Mediterranean area, with limited representation in Alpine regions such as Trentino. This imbalance may have constrained the model’s ability to accurately detect the species in marginal areas of its range. (A) 20 Figure 1.3. Maps of the actual distribution, Actual Natural Vegetation (ANV), (A) and potential distribution, Potential Natural Vegetation (PNV), (B) of Quercus ilex (holm oak), clipped for the Province of Trento. The spatial resolution is high, at 30 meters. It is evident that certain areas—such as Val dei Laghi (circled), where the species is endemic—were not accurately predicted by the study, particularly in the case of the actual vegetation map (A). Source: Bonannella et al. (2022, supplementary material), processed with QGIS. (B) 21 One of the most significant and recent contributions in Ecological Niche Modelling at the European level is that of Mauri et al. (2022), which examines the distribution of 67 major plant species, projecting their future distributions across three time horizons (2035, 2065, and 2095) under two emission scenarios from the IPCC AR5 (RCP 4.5 and 8.5; Figure 1.5). The spatial resolution used is 5 arc-minutes (approximately 9 km). The study employs six different algorithms implemented using the Biomod2 platform in R, and incorporates a natural dispersal model (MigClim). Figure 1.4. Distribution of presence points of holm oak used to train the models. The vast majority of the occurrence points are in the mediterranean area, thus the peripheral zones, as in the case of the Autonomous Province of Trento, are excluded (blue circle). This led to discrepancies in local-scale distribution studies. Source: Bonannella et al. (2022, supplementary material). 22 (A) (B) Figure 1.5. (A) Map of the potential distribution of silver fir (Abies alba Mill.) in Europe, under climate scenario RCP4.5 in 2035. (B) The same map clipped for the Autonomous Province of Trento. The resolution (5 arc-minutes) is not suitable for local-scale management. The full dataset, including raster files for all 67 species projected under the considered climate scenarios, is available in the study's repository. Source: Mauri et al. (2022, supplementary material), processed with QGIS. 23 Among the first European studies to model both the current and future distribution of up to 1,350 European plant species is that of Thuiller et al. (2005). The predictions in this study were generated using four algorithms and four climate scenarios (Third Assessment Report IPCC) for the period 2051–2080, with a spatial resolution of 50 km. Zimmermann et al. (2013) focused on more than 50 tree species in the Alps, projecting their distribution across four time periods (present, 1991–2020, 2021–2050, and 2051–2080). They employed six Regional Climate Models (RCMs) from the IPCC Fifth Assessment Report (AR5) and six of the most commonly used statistical models. Although the spatial resolution is not explicitly stated, the authors specified that a 100 m Digital Terrain Model (DTM) was used, which is relatively high compared to the previous studies. The models were combined using the ensemble SDM approach to produce habitat suitability maps, as illustrated in Figure 1.6. Figure 1.6: Habitat suitability maps for European beech (Fagus sylvatica L.) across four time periods: one current and three future projections. Red indicates a high level of confidence in habitat suitability, with more than 60% of models predicting species presence in those areas, whereas orange reflects uncertainty, with only 30–60% of models simulating species occurrence. The spatial resolution of 100 m may be particularly valuable for regional management purposes. Source: Zimmermann et al. (2013). 24 An important study at the European continental scale is that of Takolander et al. (2019), which models the current and future distribution of 14 of the most representative forest species. The study considers four emission scenarios (SRES) across two future time periods (2021–2050 and 2071–2100), using a spatial resolution of 10 arc-minutes. Eight algorithms are employed within the Biomod2 R platform, alongside a Dynamic Vegetation Model (LPJ-GUESS), which simulates tree population dynamics, competition, and biogeochemical cycles. Another European study, conducted by Buras & Menzel (2019), modelled the distribution of the 26 most abundant tree species, projecting them across three climate scenarios: a historical baseline (1961–1990) and two future projections, RCP4.5 and RCP8.5 (2061–2090; Figure 1.7). This study used a lower resolution (30 arc-minutes, approximately 55 km). Thurm et al. (2018) analysed the European distribution of 12 thermophilous tree species, along with Norway spruce and European beech as reference species, to assess whether these species could effectively replace those currently most widespread and utilized, but increasingly threatened by climate change (Figure 1.8). For their statistical analyses, the authors employed Generalized Additive Models (GAM) and Linear Regression Models (LM). The study considered two climate scenarios (RCP4.5 and RCP8.5) and two time periods (present and 2061–2080), at a spatial resolution of 1 km. The analysis was based on 7 climatic variables and 13 geological variables. Figure 1.7. Historical (1961–1990) and future (2061–2090) projections of the distribution of Scots pine (Pinus sylvestris L.) in Europe under two climate scenarios (RCP4.5 and RCP8.5). The spatial resolution is very low (approximately 55 km), allowing for analysis only at the continental scale. Source: Buras & Menzel (2019). 25 Other ENM studies notable in Europe are: Chakraborty, Móricz, et al. (2021), Dyderski et al. (2018, 2025), and Noce et al. (2017). At the Italian scale, the work of Noce et al. (2023) is particularly noteworthy. All these works modelled the future distribution of tree species using multiple algorithms at a spatial resolution greater than 1 km. For a thorough and up-to-date review of the literature on ENM applied to forest species, refer to Xu et al. (2025). In terms of estimating community parameters as species richness, explained in the following Section 1.4, fewer examples are found in the literature. Among a few others, Bedair et al. (2023) conducted a community-level study along the northwestern coastal strip of Egypt, an area known for high biodiversity and the presence of Mediterranean endemic species. The work employed the Stacked Species Distribution Modelling (SSDM) approach to model species richness using a spatial resolution of 30 arc-seconds (approximately 1 km). Predictions were generated using an ensemble modelling framework combining Random Forest and MaxEnt algorithms. The result is presented in Figure 1.9. Figure 1.8: (a) Current distribution of sweet chestnut (Castanea sativa Mill.) in Europe. (b) Chestnut distribution in 2070 under climate scenario RCP4.5. A general northward expansion of the species' distribution can be observed (probability > 0.5, colour scale from orange to red). However, the 1 km resolution does not allow for regional-scale analysis. Source: Thurm et al. (2018). 32 - ggspatial (Dunnington, 2023): Used to plot the habitat suitability maps. Throughout this work, QGIS (v. 3.32.3-Lima; QGIS.org, 2025) was employed to refine and display model outputs, while GitHub (2025) served as the platform for retrieving opensource code. Two distinct computational environments, provided by the CIRGEO department at the University of Padova, were used: a Windows computer equipped with 32 cores and 137 GB of RAM, and a High-Performance Computing (HPC) machine featuring 384 cores and 811 GB of RAM. Assistance with coding and technical writing was supported in part by OpenAI's ChatGPT version 3.5 and Google’s Gemini version 2.5 Pro. Additional materials included the results, the computational workflow and the codes are uploaded in Zenodo EU platform (link in Supplementary Material). 2.3 Species Occurrence and Environmental Predictor Datasets In this Ecological Niche Modelling study, local presence data for the 23 target species were collected, along with climatic data from the ECLIPS2.0 dataset (Chakraborty, Dobor, et al., 2021), topographic data derived from the Digital Terrain Model (DTM), soil-related features from the provincial Geological Map, and disturbance layers from the provincial Hazards Map. 2.3.1 Species Occurrence The main presence dataset used was derived from the Forest Unit (UFor) map, which represents the basic units for describing and interpreting the landscape and consists of homogeneous forest stand segments (PAT, 2016). The UFor polygons contain information on the percentage presence of the main forest species within each area. The complete file, provided by the Trentino Forest Service (Servizio Foreste della PAT) in shapefile format, includes all digitized polygons based on the local Forest Management Plans of nearly all managed forest properties in Trentino. For each species, the analysis selected UFor polygons in which the species accounted for at least 10% of the cover. Subsequently, a high-resolution (10 m) map of the current distribution of species and formations was created, using Sentinel-2 imagery, digital terrain models (DTMs), forest canopy height data obtained from LiDAR scans, and climatic information. The mapping included five types of vegetation. The UFor polygons were then overlaid onto the current vegetation map. This map was subdivided into 18 raster layers at a spatial resolution of 50 meters, each 33 containing presence and absence data for the respective species (absence data were not used in the study). If a polygon did not match the species indicated on the map, it was discarded; if a match was found, only the portion of the polygon corresponding to the species' presence was retained. Finally, a grid with 50 m × 50 m cells was generated, and grid nodes falling within the polygons of each species were extracted. The data were saved in CSV format, including the species name and X (easting) and Y (northing) coordinates in the EPSG:25832 reference system – ETRS89 / UTM zone 32N. For species with more than 6,000 points (n = 13), spatial thinning was applied using the thin_points function from the R package GeoThinneR. The process was carried out in three stages: • Preliminary reduction: points were reduced using the “brute force” method. • Block-wise thinning: the dataset was divided into chunks to optimize RAM usage and additional thinning was applied to each chunk. • Final aggregation: the thinned chunks were merged into a single spatial dataset. The parameters chunk_size and thin_dist were set manually, with various configurations tested to obtain a final number of points between 4,200 and 5,800 (Table 2.1) which were considered optimal to ensure accurate spatial predictions. The thinning step was essential for reducing the severe computational load in training the models, especially with over 6,000 occurrence points. Notably, even with access to a high-performance computing (HPC) infrastructure, the ensemble species distribution modelling phase (Section 2.4.1) remained highly time-consuming. This was largely due to the combination of fine spatial resolution, the relatively large extent of the study area, and the high number of target species, all of which imposed substantial computational demands. Another database used in this study was provided by the Edmund Mach Foundation (FEM), containing a set of highly accurate geolocated presence points (precision < 1 m) collected between 2007 and 2021. This dataset was used to supplement the Forest Unit data for seven minor species with fewer than 1000 points in the UFor (Betula pendula Roth, Prunus avium (L.) L., Acer pseudoplatanus L., Fraxinus excelsior L., Populus tremula L., Sorbus aucuparia L., Tilia cordata Mill.), which were selected, with several others, also due to their suitability for reforestation, as specified by the DGP No. 1343 of July 22, 2022 (Table 2.1). iNaturalist (iNat) database: iNaturalist is an online citizen science platform that collects observations of living species from around the world. Its data are widely used in 34 scientific research, including ENM studies (deboas, 2024), with more than 6,100 publications to date referencing this database (GBIF, 2025). The presence points for the seven minor species were downloaded using the R package rinat. The retrieved data, which included species names and observation coordinates, were clipped to the boundaries of the provincial territory, checked for taxonomic errors and duplicate records, and then saved in CSV format. These records were also used to supplement those derived from the Forest Unit dataset, which provides limited coverage for minor species. Table 2.1 Number of original presence points for each dataset (Edmund Mach Foundation, FEM; iNaturalist, iNat; and Forest Units, UFor), as well as the final number of points used after thinning for species with more than 6,000 points. The parameters Thin_distance and Chunk_size were adjusted during the thinning process using the R package GeoThinneR in order to obtain a thinned dataset ranging between 4,200 and 5,800 points. In bold the species authorized for the reforestation by the DGP No. 1343 of July 22, 2022. Species Presence dataset Number of points used (postthinning) Thin_distance (km) Chunk_size FEM iNat UFor Picea abies (L.) H.Karst. 225,830 5,832 15 20,000 Larix decidua Mill. 142,551 5,068 18 20,000 Fagus sylvatica L. 115,699 4,536 15 18,000 Ostrya carpinifolia Scop. 75,105 4,531 11 12,000 Fraxinus ornus L. 70,510 4,536 2 10,000 Abies alba Mill. 63,067 4,526 2 10,000 Pinus sylvestris L. 55,564 4,450 7 9,000 Quercus pubescens Willd. 44,610 4,523 7 6,000 Aria edulis (Willd.) M.Roem. 27,488 4,531 5 4,000 Pinus cembra L. 19,953 4,494 4 3,000 Pinus mugo Turra 16,745 4,466 4 2,500 Alnus alnobetula (Ehrh.) K.Koch 14,755 4,814 4 2,000 Pinus nigra J.F.Arnold 10,209 4,206 4 1,500 Castanea sativa Mill 4,095 4,095 Quercus petraea (Matt.) Liebl. 3,402 3,402 Quercus ilex L. 1,488 1,488 Betula pendula Roth 711 63 281 1,055 Prunus avium (L.) L. 311 42 275 628 Acer pseudoplatanus L. 454 131 206 791 Fraxinus excelsior L. 494 51 175 720 Populus tremula L. 686 45 158 889 Sorbus aucuparia L. 367 157 149 673 Tilia cordata Mill. 175 33 58 266 35 For each of the seven less-represented species, the three presence datasets were aggregated and initially organized by species. The CSV files were then standardized to retain only the columns containing the scientific species name and the X and Y coordinates of the presence points. 2.3.2 Climatic Variables This study utilised climatic data from the European CLimate Index ProjectionS (ECLIPS) dataset (Chakraborty, Dobor, et al., 2021), specifically from its high-resolution ECLIPS 2.0 version. This dataset includes spatially downscaled climate projections derived from the EURO-CORDEX ensemble, providing raster maps at a resolution of 30 arc-seconds (approximately 1 × 1 km), distributed in GeoTIFF format. From the 42 available climatic indices in ECLIPS 2.0, 11 variables were selected for analysis, based on their ecological relevance in modelling forest species distributions. The selected indices are detailed in Table 2.2, which includes their acronyms, units of measurement, and calculation methods. These variables capture key climatic dimensions such as temperature, precipitation, seasonality, and moisture availability. Table 2.2. Selected Climatic Indices for the Ecological niche modelling. The table shows the eleven climatic indices utilized as predictive variables for modelling both current and future forest species distributions. For each index, its unit of measurement and the methodology for its calculation are provided. Source: Chakraborty, Dobor, et al. (2021, supplementary material), modified. Acronym Variable name Unit Calculation MWMT Mean warmest month temperature °C MWMT = (Tmon.mean) MCMT Mean coldest month temperature °C MCMT = (Tmon.mean) TD Continentality °C MWMT - MCMT AHM Annual heat: moisture index °C/mm SHM Summer heat: moisture index °C/mm DDbelow0 Degree-days below 0°C °C DDabove18 Degree-days above 18°C °C NFFD Number of frost-free days - EMT Extreme minimum temperature °C MAP Annual total precipitation mm Summary of daily precipitation over a year Tmax_an Annual mean of maximum temperature °C Average of daily maximum temperature over a year 36 Climatic data were extracted for the baseline period 1991–2010 – referred throughout the work as present or current climate – as well as for three future time periods: 2041– 2060, 2061–2080, and 2081–2100. Future climate projections were obtained under two Representative Concentration Pathways (RCPs): RCP4.5 and RCP8.5, which represent alternative greenhouse gas concentration trajectories developed by the Intergovernmental Panel on Climate Change (IPCC) in its Fifth Assessment Report (AR5) (IPCC, 2014). The nomenclature of these pathways is directly linked to their projected radiative forcing, defined as the change in the net energy balance of the Earth-atmosphere system caused by an external factor, expressed in Watts per square meter (W/m²) (IPCC, 2014). Specifically, a positive radiative forcing indicates that more energy, primarily from solar radiation, is being absorbed by the Earth system than is being radiated back to space, leading to a warming effect (Gregor, 2025). Under this framework: • RCP4.5 assumes a moderate emissions and mitigations scenario, stabilising the radiative forcing at 4.5 W/m² by 2100; • RCP8.5 represents a high-emission scenario with no mitigation (business-asusual), resulting in a radiative forcing of 8.5 W/m² by the end of the century. Each scenario is associated with projected increases in global mean surface temperature, as summarised in Table 2.3 and illustrated in Figure 2.2. These projections provide a basis for evaluating potential impacts on forest ecosystems under different emission pathways. Table 2.3. Projected Global Mean Surface Temperature Increases Relative to 1986–2005 Baseline. The table illustrates the anticipated warming trends under various Representative Concentration Pathways (RCPs: RCP2.6, RCP4.5, RCP6.0, and RCP8.5) for two distinct future periods: 2046–2065 and 2081–2100. The 'Likely range' quantifies the uncertainty as the 5th to 95th percentile interval derived from climate model simulations. Source: IPCC (2014), modified. 2046–2065 2081–2100 Scenario Mean (°C) Likely range (°C) Mean (°C) Likely range (°C) RCP2.6 1.0 0.4 - 1.6 1.0 0.3 - 1.7 RCP4.5 1.4 0.9 - 2.0 1.8 1.1 - 2.6 RCP6.0 1.3 0.8 - 1.8 2.2 1.4 - 3.1 RCP8.5 2.0 1.4 - 2.6 3.7 2.6 - 4.8 37 To account for uncertainties in climate modelling, projections for each RCP and time period were derived from five Regional Climate Models (RCMs), whose specifications are listed in Table 2.4, along with their corresponding driving Global Climate Models (GCMs). For each combination of RCP and time period, the relevant raster maps from the five RCMs were averaged – using the raster package in R – to produce a consensus climate projection. Table 2.4. Regional Climate Models (RCMs) and their driving Global Climate Models (GCMs) employed in the study. These models contributed to the averaged climate rasters utilized for modelling current and future species distributions. Source: Chakraborty, Dobor, et al. (2021), modified. All raster layers were subsequently resampled to a spatial resolution of 50 m using the nearest neighbor algorithm via the R package terra. The coordinate reference system was also standardised by converting the spatial reference from EPSG:22232 to EPSG:25832, ensuring consistency with other geospatial datasets used in the study. Institute RCM Driving GCM Climate Limited-Area Modelling Community (CLM-Community) CLMcom-CLM4-8-17 CNRM-CERFACS-CNRM-CM5 Climate Limited-Area Modelling Community (CLM-Community) CLMcom-CLM4-8-17 MPI-M-MPI-ESM-LR Danish Meteorological Institute (DMI) DMI-HIRHAM5 ICHEC-EC-EARTH Royal Netherlands Meteorological Institute (KNMI) KNMI-RACMO22E MOHC-HadGEM2-ES Helmholtz-Zentrum Geesthacht, Climate Service Center, Max Planck Institute for Meteorology (MPI-CSC) MPI-CSC-REMO2009 MPI-M-MPI-ESM-LR Figure 2.2. Annual anthropogenic CO₂ emissions and global average surface temperature change under IPCC RCP scenarios, compared with historical data. Panel (a) shows historical and projected CO₂ emissions. Panel (b) shows projected temperature changes (2081-2100) relative to 1986-2005 for RCP scenarios. Shaded areas: likely range. Source: IPCC (2014). 38 2.3.3 Topographic Variables Topographic factors, alongside climatic variables, play a crucial role in determining plant‐species distributions. All indices were derived from a 10 m Digital Terrain Model (DTM) obtained via the WebGIS portal of the Autonomous Province of Trento (PAT, 2025b) and resampled to 50 m resolution using the nearest neighbor algorithm in QGIS. The selected variables are: - Slope: calculated in QGIS (degrees), it indicates terrain steepness and its influence on soil drainage and erosion. - Terrain Ruggedness Index (TRI): computed using QGIS, it quantifies the elevation difference between adjacent cells. - Heatload Index (HLI): obtained using slope, aspect and latitude, it quantifies the potential thermal energy received at a site in a year (Zelený & Lin, 2021). This index substitutes raw aspect, – which presents critical issues, for example, 1° and 360°, although numerically opposite, both correspond to north-facing – by providing a continuous, ecologically meaningful measure of thermal exposure. - Topographic Wetness Index (TWI): it estimates soil-moisture accumulation based on the topography. The last two variables were extracted using the R package spatialEco (Evans & Murphy, 2023). Elevation was excluded as it does not directly drive species' occurrences but is rather correlated with other variables, such as mean annual temperature. While future projection models might partially base habitat suitability on elevation, it is crucial to consider that past climatic conditions associated with a specific altitude may no longer persist, potentially introducing bias and undermining the model's predictive accuracy. 2.3.4 Soil Variables Soil-related features, important for understanding plant distribution, were derived from the provincial Geological Map (PAT, 2024). Specifically, soil pH was categorized into three binary raster layers. In these layers, a value of 1 signifies the presence of a specific soil property, while a value of 0 indicates its absence: - Acidic Soil: areas with predominant acidic soil conditions. - Basic Soil: areas with predominant basic soil conditions. 39 - Mixed Soil: areas where both acidic and basic soil types are present, or where conditions are neutral (pH ≈ 7). These raster layers were resampled to a 50 m resolution to align with other environmental predictors in the study. 2.3.5 Disturbance Variables Disturbance variables were incorporated to account for geomorphological processes that could influence the species distributions. Two specific disturbance types were included as raster layers, sourced from the WebGIS portal (PAT, 2025c): - Rockfall: A layer indicating areas prone to rockfall events; - Landslides: A layer indicating areas affected by landslide phenomena. These layers were initially rasterized and then binarized to represent hazard levels. A value of 0 was assigned to areas classified as "no hazards," "negligible," and "low" hazard, while a value of 1 corresponded to areas with "medium" and "high" hazard levels. 2.4 Model Framework The modelling framework adopted in this study relies on the integration of single-species distribution models (ESDMs) into community-level predictions using the R package SSDM. The process includes the preparation of environmental predictors and species occurrence data (Section 2.3), the construction of ensemble models for each species, and the aggregation of these into stacked models to generate spatial predictions of species richness. A schematic overview of the entire workflow is presented in Figure 2.3. The diagram outlines key phases including data preparation, modelling, spatial outputs, and the elaboration of case-study examples. 40 The initial phase of the workflow consisted of two main steps. First, all raster layers of environmental and climatic variables – named according to their acronyms (see Section 2.3.3 and Table 2.2) – were stacked using the load_var function from the R package SSDM. In ecological niche modelling studies, it is crucial that all raster files share the same coordinate reference system, spatial extent and resolution (Guisan & Thuiller, 2005; Guisan & Zimmermann, 2000; Hijmans & Elith, 2025). The SSDM package automatically manages discrepancies in spatial extent and resolution between input layers, ensuring compatibility across all variables (Schmitt et al., 2017). Figure 2.3. Workflow diagram of the ecological niche modelling process followed in the present work. The flowchart illustrates the main parts of the work’s process, starting from the occurrence and environmental data, to the final outputs. In italic are written the main functions from the R package SSDM that were used to build the models. The 23 species are illustrated in Table 2.1. The 2 future scenarios are: intermediate (RCP4.5) and business-as-usual (RCP8.5). The time periods are: 2041-2060, 2061-2080, and 2081-2100. 41 The second step of this phase involved a thinning process applied to species occurrence data, primarily to reduce computational effort. This step is described in detail in Section 2.3.1. 2.4.1 Ensemble Species Distribution Models Ensemble Species Distribution Modelling (ESDM) combines multiple statistical algorithms, both machine learning and regression-based, to generate robust predictions of species distributions. By integrating outputs from different modelling approaches, ESDMs mitigate algorithm-specific biases and enhance ecological reliability (Araújo & New, 2007; Benito et al., 2013). In this study current habitat suitability of the 23 target forest species was modelled using the ensemble_modelling function from the package SSDM. This function requires two primary inputs: (i) a raster stack of environmental predictors and (ii) a CSV file containing the X and Y coordinates of species presence points (Table 2.1). Absence data were not required, as the model function automatically generates pseudoabsences beyond the known presence locations (Schmitt et al., 2017). Eight algorithms were selected for ensemble modelling: • Random Forest (RF): A machine-learning method that constructs numerous decision trees and outputs the mode of predictions. It captures complex, nonlinear interactions, handles large datasets, and is robust to overfitting and noisy data (Friedman, 2001; Hastie et al., 2009); • Generalized Linear Model (GLM): Extends linear regression by allowing nonnormal error distributions. In ENM, it models species presence-absence relationships and provides interpretable coefficients for predictor influence (Elith & Leathwick, 2009; Friedman, 2001; Guisan & Zimmermann, 2000); • Generalized Additive Model (GAM): Enhances GLMs by incorporating smooth, non-parametric functions to capture nonlinear species-environment relationships, improving flexibility in ecological modelling (Hastie et al., 2009); • Classification Tree Analysis (CTA): A decision-tree algorithm that forecasts presence or absence using environmental thresholds, offering transparent rules for interpreting species responses (De’Ath & Fabricius, 2000; Hastie et al., 2009); 48 Quercus pubescens expected to show an expansion under moderate conditions (Section 3.3). This targeted approach provides a clear and concise example of the typical outputs generated by the SSDM package, including continuous habitat suitability maps, binary presence–absence projections, uncertainty estimates, variable importance plots, and algorithm performance metrics (e.g., AUC, Kappa). Their selection serves as a representative illustration of the broader ensemble modelling process. All the results, included the other species, are available in the Zenodo digital repository (link in Supplementary Material). 49 3 RESULTS 3.1 Projected Ecological Niches of 23 Forest Tree Species in Trentino This section presents the primary findings from the ecological niche modelling of 23 forest tree species within the province of Trento. The analysis is structured into three main sections. First, it examines the projected changes in suitable habitat area and altitudinal shifts for each species, comparing the current scenario with future projections under two Representative Concentration Pathways (RCP 4.5 and RCP 8.5) across three distinct time periods (2041–2060, 2061–2080, and 2081–2100). Second, the section analyses species richness at the community level as an aggregate indicator of forest biodiversity, evaluating its evolution in response to climate change. Finally, two detailed case studies are presented for downy oak (Quercus pubescens) and Swiss stone pine (Pinus cembra), selected for their contrasting projected responses to future climatic conditions. All the results are uploaded in Zenodo platform (link in Supplementary Material). The projected impacts of climate change on the suitable habitat area for the 23 studied forest tree species are presented in Figure 3.1. The responses vary significantly among species, allowing for their classification into three categories based on their projected trends: - "Gainers": This group includes predominantly thermophilous species that are expected to expand their suitable habitat. Species such as common whitebeam (Aria edulis), manna ash (Fraxinus ornus), European hop-hornbeam (Ostrya carpinifolia), black pine (Pinus nigra), holm oak (Quercus ilex), and downy oak (Quercus pubescens) show a consistent positive trend across future scenarios and time periods. Notably, the Mediterranean species holm oak is projected to increase in habitat suitability by an order of magnitude. The mountain pine (Pinus mugo) and the European beech (Fagus sylvatica) will also have a positive trend throughout the decades. The sweet chestnut (Castanea sativa), the wild cherry (Prunus avium), the sessile oak (Quercus petraea) and the small-leaved lime (Tilia cordata) can also fit into this category, despite their high suitability areas are expected to decline. 50 - "Losers": This category comprises species facing a stable reduction in suitable habitat area. The decline is especially dramatic for the two high-altitude, coldadapted species Norway spruce (Picea abies) and Swiss stone pine (Pinus cembra). For this last species in particular the projections indicate a near-total disappearance of its habitat under almost all future scenarios, highlighting its extreme vulnerability. A strong decline is expected also for the common ash (Fraxinus excelsior), while silver fir (Abies alba) and Scots pine (Pinus sylvestris) will face a prevalent reduction of their high suitable areas. - "Different behaviours": A third group of species, including sycamore maple (Acer pseudoplatanus), green alder (Alnus alnobetula), Silver birch (Betula pendula), European larch (Larix decidua), European aspen (Populus tremula), and rowan (Sorbus aucuparia) exhibits non-linear temporal trend, with initial gain offset by consequent losses, especially in the business-as-usual scenario RCP8.5. The gains in higher elevations are compensated by the losses at lower elevations. 51 52 Figure 3.1. Changes in suitable surface area (in km²) for the 23 studied species. across the two RCP climate scenarios (4.5 and 8.5) in the three future time periods (2041-2060, 2061-2080 and 2081-2100), compared with the present scenario. Responses to climate change are very diversified, with the extremes being the Pinus cembra that is projected to almost disappear in all the scenarios, and Quercus ilex that is projected to increase its potential habitat suitability by one order of magnitude. 53 The absolute differences in the suitable surface area of the future scenarios and time periods, compared with the present climate, are illustrated in Figure 3.2. The holm oak (Quercus ilex) is the species with the most important increase in potential suitable area, for all scenarios and timeframes, gaining, in the intermediate-emission scenario RCP4.5, about 2000 km2 by the end of the century. In absolute terms the major suitable area loss will be experimented by the Norway spruce (Picea abies) in all scenarios and time periods (-1500 km2 by the end of the century in the RCP4.5). Figure 3.2. Absolute difference in suitable surface area in km2 compared with the present climate. 54 In relative terms, for the holm oak the gain in relative potential suitable area is evident in all scenarios (Figure 3.3). By the end of the century for this species the potential suitable area, considering the medium-emission scenario, could increase by 20 times. The major loss will be experimented by the Swiss stone pine, that will lose 100% of its area in all the suitability classes by the end of the century (Section 3.3.1). Figure 3.3. Relative difference in suitable surface area compared with the present climate. 55 The analysis of altitudinal shifts – for the medium and high suitability areas – revealed a general upward migration trend for all the 23 species (Figure 3.4). For instance, the Norway spruce, currently abundant at 1200 - 1800 m a.s.l., is projected to shift its optimal suitability to 1600 - 2100 m a.s.l. under RCP4.5 by the end of the century while the holm oak, now with the optimum at 250 – 750 m a.s.l. will expand its suitable habitat beyond 1000 m a.s.l. 56 The average elevation shift for the medium and high suitable areas for all the 23 species is 83 m by the mid-century and 265 m by the end of the century, considering the RCP4.5 scenario (Table 3.1). High variance across the species is registered (the standard deviation Figure 3.4. Variation in the altitudinal range of 23 tree species across four temporal periods (Present, 20212040, 2041-2060, 2081-2100) and two emission scenarios (RCP4.5 and RCP8.5). Only medium and high suitability classes are included in the plot. Altitudinal distributions are presented using violin plots, which indicate the probability density for each species, scenario, and period. The mean is represented by the line, the median by the dot. 57 is 83 and 75 m respectively). Under the RCP8.5 scenario, altitude shifts are anticipated to double, comparing of those projected under the intermediate scenario. Table 3.1. Mean elevation gain (in m a.s.l.) for the two emission scenarios in the three time-periods referred to the present scenario. The overall mean indicate elevation shift across all the 23 species. RCP4.5 RCP8.5 Species Present 2041-2060 2061-2080 2081-2100 2041-2060 2061-2080 2081-2100 Abies alba 1247 18 180 284 189 456 809 Acer pseudoplatanus 1018 16 171 201 159 369 521 Alnus alnobetula 1898 56 154 232 162 367 667 Aria edulis 750 97 181 248 189 358 543 Betula pendula 1062 200 334 434 338 656 989 Castanea sativa 762 36 225 305 238 582 1011 Fagus sylvatica 1143 85 220 309 233 480 804 Fraxinus excelsior 917 14 173 237 184 369 651 Fraxinus ornus 725 62 140 195 142 289 467 Larix decidua 1708 29 82 151 112 263 556 Ostrya carpinifolia 733 66 145 208 154 290 473 Picea abies 1529 41 170 291 183 445 - Pinus cembra 1968 414 459 - - - - Pinus mugo 1489 58 160 233 155 342 509 Pinus nigra 592 50 168 210 155 338 520 Pinus sylvestris 857 9 177 237 168 348 625 Populus tremula 969 125 267 299 259 605 839 Prunus avium 681 -4 93 165 111 340 379 Quercus ilex 464 101 232 271 223 411 603 Quercus petraea 807 66 190 291 216 474 920 Quercus pubescens 689 89 185 226 185 350 520 Sorbus aucuparia 1272 150 278 387 306 587 880 Tilia cordata 720 126 281 418 296 516 Mean (m a.s.l.) 1043 83 203 265 198 420 664 Standard deviation 423 88 81 75 60 112 190 The relative importance of predictor variables differs notably among the analysed species (Figure 3.5). Climatic and topographic variables generally show the highest contributions. For some species, topographic variables were found to be the most important. The most prominent example is the Topographic Ruggedness Index (TRI), which is the most influential variable for six species (Aria edulis, Betula pendula, Castanea sativa, Fraxinus ornus, Ostrya carpinifolia, and Quercus petraea), with an importance constantly over 15%. The heatload index is the most important factor for Alnus alnobetula, 64 Valsugana Valley). Notably, the model predicted a strong reduction in species richness for area with predominant acid soil (e.g. Lagorai mountain chain area, east of the study area). (A) (B) 65 (D) (C) 66 (E) (F) 67 The projected impact of climate change on the spatial distribution of species richness varies significantly depending on the climate scenario considered (Figure 3.8). Under the intermediate scenario, RCP4.5, the distribution of species richness is projected to remain relatively stable throughout the 21st century. The overall shape of the curve, which peaks with a richness of two species occupying approximately 1,400 km², shows minimal deviation across the future time periods compared to the present-day baseline. This pattern suggests a certain resilience of the species assemblage to the climatic changes expected under this scenario. In stark contrast, the business-as-usual RCP8.5 scenario is projected to cause severe alterations to the species richness landscape. A clear trend of biodiversity decline is evident, characterized by two primary phenomena. Firstly, a progressive loss of areas with high species richness is projected; for instance, areas supporting six or more species are expected to shrink substantially and nearly vanish by the 2061-2080 period. Secondly, and most dramatically, a shift towards landscapes dominated by low species richness is predicted. This is highlighted by the sharp increase in the area suitable for only (G) Figure 3.7. Future species richness map under present (A) and two future climate scenarios in three time periods: RCP4.5 for 2041–2060 (B), 2061–2080 (C), and 2081–2100 (D); and RCP8.5 for 2041–2060 (E), 2061–2080 (F), and 2081–2100 (G). 68 a single species, which is projected to grow from approximately 900 km² in the present day to over 1,500 km² by the 2081-2100 period (Figure 3.9). Such a trend indicates a future of increasing biotic homogenization under a high-emissions trajectory. Figure 3.9. Distribution of species richness under two climate change scenarios. The plots show the total surface area (km²) occupied by a given number of species (species richness). The analysis is presented for the two RCP scenarios (4.5 and 8.5) across four time periods: Present, 2041-2060, 2061-2080, and 20812100. Figure 3.8. Changes in mean species richness in the entire study area: confront between the RCP4.5 and RCP8.5 scenarios. 69 The mean species richness error (mean absolute error per cell) is 3.14 species for the present (Std Dev = 2.19) and decreases under future scenarios – reaching 1.39 for RCP8.5 (2081–2100) – with intermediate values for the other time frames (Table 3.4). These values represent the average absolute difference between predicted and observed species richness per raster cell, computed from the retained occurrence data in the holdout cross-validation (75/25 split, 5 replicates); full details are reported in Section 2.4. Table 3.4. Mean species richness error of the stacked SDMs across two scenarios and four time periods. 3.3 Case Studies Among the 23 species modelled, Pinus cembra (Swiss stone pine) and Quercus pubescens (downy oak) were selected as representative examples to illustrate the outputs and interpretative potential of the ESDMs. These species were chosen due to their ecological relevance and contrasting projected trends under future climate conditions in the Province of Trento: The Swiss stone pine is expected to experience a severe reduction in suitable habitat, while the downy oak is projected to benefit from future climatic shifts. This section serves as a practical demonstration of how model outputs – developed in the Section 2.4.1 – can be analysed and interpreted for ecological and conservation purposes. Results from the ESDM comprise: habitat suitability maps, binary presence/absence maps, algorithm correlation heatmaps, and uncertainty estimations. The graphs showing the suitable area changes, the elevation shifts and the variable importance ranking were presented in Section 3.1. Mean species richness error Std Dev Present 3.14 2.19 RCP4.5 2041-2060 3.00 2.07 RCP4.5 2061-2080 2.55 1.93 RCP4.5 2081-2100 2.62 2.07 RCP8.5 2041-2060 2.84 2.04 RCP8.5 2061-2080 2.21 2.13 RCP8.5 2081-2100 1.39 2.04 70 3.3.1 Pinus cembra The habitat suitability maps for the Swiss stone pine project a severe and rapid spatial contraction of its climatic niche under all future scenarios (Figure 3.10, A – G). The current distribution shows extensive areas of optimal (dark red; 385 km2) and medium (orange/yellow; 221 km2) suitability at the highest elevations of the study area, as well as other 329 km2 of low-suitability areas (green). This trend of decline begins in the near-term and becomes progressively more pronounced over time. Even under the intermediate RCP4.5 scenario, optimal and medium-suitability habitats almost completely disappear by mid-century (-99.9%), replaced by fragmented patches of low suitability (reduced by 76.8% by the decades 2041 – 2060). The reduction is greatly accelerated and more profound under the business-as-usual scenario, leaving only small low-suitability areas (- 96% in 2041 – 2060). In these scenarios, the study area becomes entirely unsuitable for the species by the end of the century, indicating a high risk of regional extinction for Swiss stone pine due to climate change. (A) 71 (B) (C) 72 (D) (E) 73 (F) (G) Figure 3.10. Habitat suitability maps of Pinus cembra under present (A) and future (B–G) climatic scenarios. Suitability is represented as a continuous gradient from unsuitable (dark blue) to optimal (dark red). A CORINE land cover mask has been applied to exclude areas unlikely to be colonized. Future projections are shown for two emission scenarios across three time periods: RCP 4.5 for 2041–2060 (B), 2061–2080 (C), and 2081–2100 (D); and RCP 8.5 for 2041–2060 (E), 2061–2080 (F), and 2081–2100 (G). 80 3.3.2 Quercus pubescens The habitat suitability maps for the downy oak indicate that this thermophilous species is projected to be a significant beneficiary of future climate change in the region. Currently, its optimaland medium-suitability areas (514 and 379 km2 respectively) are largely confined to the low-elevation valley floors (Figure 3.14, A). Under all future scenarios, a notable expansion of suitable habitat is projected (Figure 3.14, B – G). This change is characterized by a significant upward altitudinal shift, with high-suitability zones expanding from the valley floors to occupy the mid-elevation slopes that are currently suboptimal or marginal. For RCP4.5 scenario in 2061 – 2080 is expected an increase of 57.8% for the high-suitability areas, +99.0% for medium-suitability areas and +48% for lowsuitability areas. This trend of habitat expansion becomes progressively more pronounced over time and is substantially more rapid and extensive under the high-emission RCP8.5 pathway. By the end of the century for both emission scenarios, optimal conditions for downy oak are projected to dominate much of the region's valleys, included the Inneralpine ones like Fiemme, Fassa and Sole valleys. (A) 81 (B) (C) 82 (E) (D) 83 (F) (G) Figure 3.14 Habitat suitability maps of downy oak under present (A) and future (B–G) climatic scenarios. Suitability is represented as a continuous gradient from unsuitable (Unsuitable, dark blue) to optimal (dark red). A CORINE land cover mask is applied in order to exclude unlikely colonizable areas (see Section 2.4.1 for the classes). Future projections are shown for two emission scenarios across three time periods: RCP4.5 for 2041–2060 (B), 2061–2080 (C), and 2081–2100 (D); and RCP8.5 for 2041–2060 (E), 2061– 2080 (F), and 2081–2100 (G). 84 Under the intermediate RCP4.5 scenario, the projected upward elevation shift for downy oak is +89 m by mid-century and +226 m by the end of the century (Table 3.1). These shifts are projected to be more than double under the business-as-usual RCP8.5 scenario, reaching +185 m and +520 m, respectively. To translate the continuous habitat suitability projections into binary presence/absence maps, a threshold of 0.70, calculated using the SES method, was applied. The resulting map for the current period shows the predicted presence of the downy oak is largely confined to the low-elevation valley floors and sun-exposed lower slopes (Figure 3.15). The projection for the 2061–2080 period under the RCP4.5 scenario reveals a significant expansion of the species' potential range. This expansion is characterized by a distinct upward altitudinal shift, with the predicted presence extending further up the valley slopes and into adjacent areas that are currently classified as unsuitable. This result provides a clear spatial representation of the gains in suitable area for Quercus pubescens, reinforcing the finding that it is a likely beneficiary of future climatic warming in the region. (A) 85 An evaluation of the individual algorithms reveals a high overall quality for the downy oak ensemble model (Table 3.6). The strength of the ensemble is supported by a mean AUC of 0.93 and a mean Cohen's Kappa of 0.71, indicating excellent discrimination and substantial agreement with observations, respectively (Table 3.6). The performance of the individual models varied. Random Forest (RF) stood out as the most consistently accurate algorithm for this species too, achieving the highest scores for both AUC (0.96) and Kappa (0.81), coupled with a strong calibration index of 0.90. In stark contrast, the Generalized Boosted Regression Model (GBM) presented a significant caveat. While its discriminative ability was high (AUC = 0.94), its severely negative calibration score (-0.29) renders its raw probability outputs untrustworthy. (B) Figure 3.15. Binary distribution maps of Quercus pubescens under present (A) and future (B–G) climatic scenarios. Red areas indicate predicted presence, while grey areas represent absence and simultaneously display the underlying Digital Terrain Model (DTM) to provide topographic context. Future projections refer to two emission scenarios across three time periods: RCP 4.5 for 2041–2060 (B), 2061–2080 (C), and 2081–2100 (D); and RCP 8.5 for 2041–2060 (E), 2061–2080 (F), and 2081–2100 (G). Even with a relatively high threshold (0.70), Quercus pubescens a is projected to increase their range by the year 2080. 86 Table 3.6. Performance evaluation of the eight individual algorithms comprising the ensemble model for Quercus pubescens. The table displays key metrics including the Area Under the Curve (AUC), omission rate, sensitivity, specificity, Cohen's Kappa, and a calibration index. All algorithms were successfully run and retained for the final ensemble (Kept model: 5/5 repetitions). For Quercus pubescens, the predictions from all eight algorithms show very high positive correlation under present-day conditions, with most pairwise correlations exceeding 0.85 (e.g., GAM and MARS at 0.95; Figure 3.16). In future projections, a decline in overall agreement is observed, again being most pronounced in the late-century RCP8.5 scenario. However, a key pattern for this species is that the decrease in concordance is disproportionately driven by the GAM. While most other algorithms remain highly correlated with each other, their agreement with GAM diminishes significantly. This is best illustrated in the RCP8.5 2081-2100 scenario: the correlation between many model pairs like CTA and GBM remains extremely high at 0.90. In stark contrast, the correlation between GAM and most other models collapses, becoming strongly negative in several cases (e.g., GAM vs. GLM at -0.34; GAM vs. CTA at -0.17). Algorithm evaluation Threshold AUC Omission rate Sensitivity Specificity Prop correct Cohen’s Kappa Calibr. Kept model ANN 0.60 0.86 0.16 0.92 0.75 0.84 0.67 0.93 5 CTA 0.64 0.91 0.14 0.86 0.86 0.86 0.72 0.95 5 GAM 0.81 0.95 0.11 0.89 0.89 0.89 0.69 0.78 5 GBM 0.52 0.93 0.14 0.86 0.86 0.86 0.71 -0.29 5 GLM 0.80 0.93 0.13 0.87 0.87 0.87 0.63 0.79 5 MARS 0.85 0.94 0.12 0.88 0.88 0.88 0.64 0.81 5 RF 0.65 0.96 0.09 0.91 0.91 0.91 0.81 0.90 5 SVM 0.77 0.95 0.11 0.89 0.89 0.89 0.77 0.82 5 ESDM Mean 0.70 0.93 0.13 0.88 0.86 0.87 0.71 0.71 5 Std Dev 0.12 0.03 0.02 0.02 0.05 0.02 0.06 0.41 0 87 (A) (B) (C) (D) (E) 88 The spatial distribution of model uncertainty for Quercus pubescens reveals a distinct altitudinal pattern (Figure 3.17, A – B). Under current conditions, areas of lowest uncertainty (i.e., highest model agreement) are found in the high-elevation alpine areas. Conversely, the greatest uncertainty is concentrated in the low-elevation floors of the main valleys and in clear bands at mid-elevations (Figure 3.17, A). These transitional zones, which represent the current lower and upper limits of the species' core habitat, show a maximum uncertainty value of 0.185. In the future projection for 2061–2080 (RCP4.5), this pattern shifts. The maximum uncertainty value across the landscape decreases slightly to 0.168 (Figure 3.17, B). Spatially, the high-uncertainty bands are no longer located at the same mid-elevations; instead, they have moved noticeably upwards in altitude. Figure 3.16. Pairwise correlation matrix of present predictions from the eight individual algorithms used in the ensemble model for Quercus pubescens. The values represent the Pearson correlation coefficient between the continuous habitat suitability predictions of each algorithm pair. The colour scale indicates the strength of the correlation, with dark blue representing high positive correlation (strong agreement) and red representing high negative correlation (strong disagreement). (F) (G) 89 (A) (B) Figure 3.17 Uncertainty of the Ensemble Species Distribution Model for Quercus pubescens. The maps illustrate the between-algorithm variance for the present (A) and a future scenario (B; RCP4.5, 2061–2080). Warmer colours indicate higher disagreement among the models. The maximum observed variance is 0.185 in the present map and decreases to 0.168 in the future projection, reflecting a growing consensus on the suitability of mid-elevation habitats as the climate warms. 96 skills – and thus fell outside the scope of this particular study – it represents a potential improvement on the workflow. A different consideration, however, applies to the availability of temporary computer memory (RAM). Particularly during the stacking of individual models and their subsequent projection into future scenarios, a minimum of 150 GB of RAM proved essential, a memory capacity not commonly available outside of dedicated HPC systems. Uncertainty maps from SSDM outputs highlight greater disagreement among the algorithms at the margins of the high and medium suitability areas, suggesting cautious interpretation in these areas. Future sensitivity analyses should test (a) multiple thresholding rules, (b) alternative stacking algorithms (e.g., probabilistic stacking vs. binary stacking), and (c) explicit dispersal scenarios (no dispersal vs full dispersal vs species-specific dispersal kernels) to quantify projection spread. It is essential to emphasize that outputs represent hypotheses of potential suitability rather than deterministic forecasts, requiring validation and adaptive updating. Key uncertainty sources included: 1) Data uncertainty from presence-only datasets. Forest Unit (Ufor) inventories may omit private lands. The Edmund Mach Foundation (FEM) and the iNaturalist data may be spatially clustered and biased toward accessible areas or subject to identification errors; 2) Model uncertainty from algorithmic variation and predictor collinearity decisions. While a collinearity test was considered, eventually a uniform predictor set across species was maintained for comparability; 3) Climate projection uncertainty – aggregating GCMs/RCMs into a mean scenario reduces variability but conceals range of plausible futures; downscaling errors also contribute. Comparing the uncertainty patterns between the cold-adapted specialist, Pinus cembra, and the thermophilic generalist, Quercus pubescens, offers a powerful insight into the challenges of forecasting species' futures. For the warmth-loving downy oak, there is a strong consensus on the direction of its response to climate change: an upward expansion in altitude. The uncertainty, which is concentrated at the leading edge of this expansion and largely driven by the outlier behaviour of the GAM algorithm, is therefore primarily about the rate and ultimate limit of this advance. In stark contrast, the future for the high-altitude Swiss stone pine is characterized by systemic and widespread uncertainty, but remains within moderate values (Figure 3.13). 97 As its climatic niche is projected to shrink and fragment, the models show broad disagreement, not just from a single outlier (Figure 3.12). This reflects a fundamental uncertainty about where, or even if, viable refugia will persist. 4.5 Study Limitations and Future Developments This section outlines the primary limitations of the study and the strategies employed to mitigate their potential impact on the findings. A comprehensive interpretation of this study's findings requires acknowledging its inherent limitations, which also pave the way for future research directions. These limitations span the ecological assumptions of the models, the nature of the input data, and the computational framework itself. 4.5.1 Latent Variables and Biotic Interactions A primary consideration is that the habitat suitability maps represent a species' potential climatic niche, not its realized niche or a guaranteed future distribution. The models are a necessary simplification of complex ecological realities and do not account for numerous external factors that can mediate species presence, such as human interventions (e.g., forest management), stochastic events (e.g., fires, storms), and the impact of pests and diseases. While this simplification is a limitation, it is also an advantage for interpretability, as it provides a clear baseline for understanding the powerful influence of climate and topography. It isolates these drivers from other confounding factors, offering a crucial "tile" in the larger mosaic of ecosystem dynamics. The models are based on climatic and topographic predictors but do not explicitly account for several latent variables that critically influence species distribution on the ground. These include: - Biotic Interactions: Factors like competition, facilitation, and predation are not modelled. For example, the projected expansion of Quercus pubescens may be limited in reality by competition from established species like Fagus sylvatica. - Soil Properties: While broad soil pH categories were included, more detailed factors like microbial communities, which enrich soil fertility, and specific soil characteristics were not considered. Regional-scale soil maps may lack the detail to capture the local variations that influence individual plant establishment and growth. Moreover, most soil maps are categorical in nature and exhibit numerous 98 soil type, rendering their incorporation into SDM frameworks challenging due to the requirement for binary encoding of categorical variables. However, alternative approaches using globally (Hengl et al., 2025) or regionally (Rota et al., 2024) available continuous datasets do exist. The projections should thus be interpreted with caution, as they may present a more favourable scenario than the actual ecological changes that will transpire (Dyderski et al., 2025), with potential suitability further restrained by biotic interaction and unsuitable soil characteristics. Integrating other ecologically and socio-economics significant variable such as the increasing demand for wood, the land cover changes and the demographic trends (Keenan, 2015) could enhance model robustness, always taking into account the risk of collinearity. 4.5.2 Occurrence Datasets and Human Influence The accuracy of any SDM is fundamentally dependent on the quality of the occurrence data. This study faces three main limitations in this regard: - Managed vs. Primary Forests: The occurrence data were collected almost exclusively from managed forests, as primary forests are not present in the province. The models, therefore, reflect species distributions as shaped by decades of human silvicultural practices, not just natural environmental filters. - Data Biases: The use of citizen-science data (iNaturalist), while valuable for supplementing records, carries potential biases related to spatial clustering in accessible areas and potential identification and georeferencing errors (Estien et al., 2024; López-Guillén et al., 2024). For the purposes of this study, however, these limitations were not considered particularly significant, as the proportion of potentially erroneous data points (due to species misidentification or inaccurate location) was deemed insufficient to affect the ESDM models. - Niche truncation: In this study, the use of local species data has produced a highly accurate model of current tree species distribution, which is invaluable for regional forest management. However, as demonstrated by Anselmetto et al. (2025), this fine-scale approach is susceptible to niche truncation when projecting into novel future climates. This limitation arises because the model is calibrated on only a partial segment of a species' full environmental tolerance, which can lead to an overestimation of the magnitude of future changes. Consequently, the model’s 99 predictions are a good reflection of the potential climatic niche but may not capture the full realized niche. Specifically, the most thermophilus species, including Q. ilex and Q. pubescens, are not currently at their thermal tolerance threshold in any area in the study region. Consequently, the model fails to adequately constrain their future suitability. This likely results in the overestimation of their suitability at low elevations under the extreme warming scenarios of RCP 8.5 during the 2081-2100 period. Furthermore, the current climatic conditions in the study area lack the continental warm valley environments that constitute typical habitat for Q. pubescens, as documented in the central and western Alps (de Pasquale & Spinelli, 2022). This absence prevents the model from accurately identifying continental valleys as potentially suitable habitat for the downy oak under projected warming scenarios. To create more robust and less biased input datasets, future studies could leverage more advanced data aggregation tools. For instance, the R package SPOCC (Owens et al., 2024) allows for the integration of multiple citizen science databases, which can help to crossvalidate records and reduce sampling bias. Future research should also aim to mitigate the niche truncation bias by adopting a data pooling methodology, integrating detailed local data with broader inventories to create more robust and reliable future projections. Many frameworks for “nesting” whole-range and sub-range distribution models have been proposed, and summarized in Guisan et al. (2025). Despite this constraint for future scenarios, the current model's fine-scale precision remains essential. It uniquely captures local anthropogenic signals and provides the detailed habitat information necessary for effective, on-the-ground conservation and management planning (Bonannella et al., 2022). Building on this work, future research should move toward more comprehensive ecosystem forecasts. The following steps are recommended: - Integrate Dynamic Processes: Future models should incorporate species-specific dispersal capacities and dynamic vegetation models (DVMs) to simulate transient responses to climate change, moving from potential distribution to projected actual distribution. - Incorporate Biotic Interactions: Explicitly modelling biotic factors, including competition among tree species and the influence of pests and pathogens (e.g., bark beetle dynamics for Picea abies), would significantly refine the accuracy of the projections. 100 - Address Niche Truncation: To mitigate the risk of niche truncation, a data-pooling methodology should be adopted, integrating the detailed local data from this study with broader, trans-regional inventories to capture a more complete representation of species' climatic tolerances. In addition, training models over a broader geographical extent – while retaining the same fine spatial resolution – and subsequently cropping them to the study area could help capture a wider portion of each species’ ecological niche, thus improving the reliability of the models. 4.6 Management and Conservation Implications There is rising interest and need in developing effective forest policies in all Europe (Nabuurs et al., 2019). Although continental and national policy frameworks are necessarily informed by coarse-resolution assessments, operational forest management is implemented at the local scale; therefore, high-resolution, management-relevant maps and models (here at 50 m) are needed (Anselmetto et al., 2025). This work can help the policy makers to identify the best suitable locations for each species, with the aim of effectively implement reforestation and assisted migration programs. For instance, Picea abies – the most important species in the province in terms of surface area and multifunctionality (PAT, 2025a), and one of the most threatened by the climate change (Dyderski et al., 2025; Section 3.1) – in most areas will need a substitute in order to guarantee a stable forest cover, that can be selected according to its future projections. Human intervention should be applied only where necessary, always prioritising natural local regeneration when feasible. Some species, selected and authorized for the reforestations (Section 2.3.1), confirmed their suitability for the future climate (e.g. Fagus sylvatica, Prunus avium, and Tilia cordata), while others must be taken into consideration carefully, since the locations where they are planted today might not be suitable in the future (e.g. Pinus sylvestris, Pinus cembra, and Fraxinus ornus). For the “Losers” and the “Different behaviours” species (Section 3.1) a thorough evaluation of the present and future suitability is advised, in order to find the best locations. The habitat suitability maps created in this study at 50 m resolution can help these reforestations policies. The implications of this research are threefold: - For forestry practice, the high-resolution suitability maps offer a valuable tool for long-term strategic planning, enabling forest managers to identify vulnerable species, 101 locate potential climatic refugia, and evaluate the suitability of tree species for future reforestation and assisted migration efforts. - For conservation policy, the findings provide a scientific basis for prioritizing conservation actions, particularly for "loser" species, and for designing resilient ecological networks. - For future research, this work establishes a methodological benchmark for regionalscale ENM and highlights the critical importance of integrating high-resolution local data. 102 103 5 CONCLUSIONS This thesis addressed a critical gap between broad-scale climate impact assessments and the operational needs of regional forest management in mountain ecosystems. The primary objective was to develop, apply, and validate a high-resolution (50 m) Ecological Niche Modelling framework to project the potential distribution of 23 dominant tree species in the Autonomous Province of Trento under future climate change scenarios. By employing a Stacked Species Distribution Model (SSDM) approach, the study successfully generated a replicable and computationally robust workflow, providing spatially explicit insights at a scale directly relevant to local conservation planning and silvicultural decisionmaking. The ecological findings reveal a significant, climate-driven restructuring of the region’s forest communities projected by the end of the 21st century. A clear divergence emerged between the future trajectories of thermophilous species, such as Quercus pubescens, which are anticipated to undergo substantial upward altitudinal expansion, and cold-adapted, high-elevation specialists, notably Pinus cembra, which face a severe contraction and high risk of local extinction. At the community level, this turnover is projected to drive an overall contained decline in mean tree species richness under the intermediate RCP4.5 scenario and a structural shift toward biotic homogenization and accelerated richness decline under the business-as-usual RCP8.5 scenario. Methodologically, the ensemble modelling approach proved to be a powerful tool for this regional application, with individual species models (ESDMs) achieving high predictive accuracy (mean AUC > 0.90), thereby validating the workflow's effectiveness for regionalscale ecological niche modelling. The stacking process (SSDM), on the other hand, requires careful thresholding, a finer tuning, and comprehensive uncertainty assessment, since the mean species richness errors were high for all the scenarios and time periods. The projections represent hypotheses of potential climatic suitability and do not incorporate dynamic ecological processes. Key factors such as species-specific dispersal rates, complex biotic interactions (e.g., competition, pest and disease outbreaks), and the potential for local genetic adaptation were not explicitly modelled. Furthermore, the reliance on presence-only data from managed forests may introduce a bias reflecting anthropogenic influences rather than the species' full fundamental niche, potentially overestimating the magnitude of future changes. 104 The ESDM/SSDM analyses produced 50 m habitat suitability and species richness maps that can directly inform forest policy, guide site-level reforestation and assisted migration, and prioritise targeted conservation areas. Future works should focus on linking these correlative outputs with dynamic or dispersal models, combining broader-scale occurrence data to reduce niche truncation, improving fine-scale predictor layers, and optimising the computational workflow to enhance robustness and operational utility. 105 Supplementary Material Reproducibility is a central requirement of modern science. All model outputs (maps and evaluation metrics), the R scripts used for data processing and modelling are publicly available in the EU research repository on Zenodo: https://doi.org/10.5281/zenodo.16781155. Bibliography Adler, C., P. Wester, I. Bhatt, C. Huggel, G.E. Insarov, M.D. Morecroft, V. Muccione, and A. Prakash (2022). Cross-Chapter Paper 5: Mountains. In: Climate Change 2022: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [H.-O. Pörtner, D.C. Roberts, M. Tignor, E.S. Poloczanska, K. Mintenbeck, A. Alegría, M. Craig, S. Langsdorf, S. Löschke, V. Möller, A. Okem, B. Rama (eds.)]. Cambridge University Press, Cambridge, UK and New York, NY, USA, pp. 2273–2318, doi:10.1017/9781009325844.022. Aitken, S. N., Yeaman, S., Holliday, J. A., Wang, T., & Curtis‐McLane, S. (2008). Adaptation, migration or extirpation: climate change outcomes for tree populations. Evolutionary Applications, 1(1), 95–111. https://doi.org/10.1111/j.1752-4571.2007.00013.x Anselmetto, N., Morresi, D., Barbarino, S., Loglisci, N., Betts, M. G., & Garbarino, M. (2025). Species distribution models built with local species data perform better for current time, but suffer from niche truncation. Agricultural and Forest Meteorology, 362. https://doi.org/10.1016/j.agrformet.2024.110361 Araújo, M. B., & New, M. (2007). Ensemble forecasting of species distributions. Trends in Ecology and Evolution, 22(1), 42–47. https://doi.org/10.1016/j.tree.2006.09.010 Barve, V., & Hart, E. (2022). rinat: Access «iNaturalist» Data Through APIs (R package version 0.1.9; Numero R package version 0.1.9). https://github.com/ropensci/rinat Bedair, H., Shaltout, K., & Halmy, M. W. A. (2023). Stacked machine learning models for predicting species richness and endemism for Mediterranean endemic plants in the Mareotis subsector in Egypt. Plant Ecology, 224(12), 1113–1126. https://doi.org/10.1007/s11258-02301366-6 Beery, S., Cole, E., Parker, J., Perona, P., & Winner, K. (2021). Species Distribution Modeling for Machine Learning Practitioners: A Review. Proceedings of 2021 4th ACM SIGCAS Conference on Computing and Sustainable Societies, COMPASS 2021, 329–348. https://doi.org/10.1145/3460112.3471966 Benito, B. M., Cayuela, L., & Albuquerque, F. S. (2013). The impact of modelling choices in the predictive performance of richness maps derived from species-distribution models: Guidelines 112 Noce, S., Collalti, A., & Santini, M. (2017). Likelihood of changes in forest species suitability, distribution, and diversity under future climate: The case of Southern Europe. Ecology and Evolution, 7(22), 9358–9375. https://doi.org/10.1002/ece3.3427 Odasso, M., Miori, M., & Gandolfo, C. (2018). I tipi forestali del Trentino: descrizione e aspetti dinamici (Servizio foreste e fauna & Provincia Autonoma di Trento, A c. Di). Owens, H., Barve, V., & Chamberlain, S. (2024). spocc: Interface to Species Occurrence Data Sources (R package version 1.2.3). https://CRAN.R-project.org/package=spocc Parmesan, C., M.D. Morecroft, Y. Trisurat, R. Adrian, G.Z. Anshari, A. Arneth, Q. Gao, P. Gonzalez, R. Harris, J. Price, N. Stevens, and G.H. Talukdarr (2022). Terrestrial and Freshwater Ecosystems and Their Services. In: Climate Change 2022: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [H.-O. Portner, D.C. Roberts, M. Tignor, E.S. Poloczanska, K. Mintenbeck, A. Alegria, M. Craig, S. Langsdorf, S. Loschke, V. Moller, A. Okem, B. Rama (eds.)]. Cambridge University Press, Cambridge, UK and New York, NY, USA, pp. 197–377, doi:10.1017/9781009325844.004. Parolo, G., & Rossi, G. (2008). Upward migration of vascular plants following a climate warming trend in the Alps. Basic and Applied Ecology, 9(2), 100–107. https://doi.org/10.1016/j.baae.2007.01.005 PAT, AA.VV. (2022). Stato di attuazione del Piano d’Azione per la gestione degli interventi di esbosco e ricostituzione dei boschi danneggiati dalla Tempesta Vaia, https://forestefauna.provincia.tn.it/content/download/16088/263412/file/Report_VAIA_Stato _finale_rid.pdf. (Last consultation 18/08/2025). PAT. (2016). Linee tecniche per la pianificazione forestale aziendale. Servizio Foreste e Fauna, https://forestefauna.provincia.tn.it/Documenti/pianificazione-forestale/Manuale-Lineetecniche-per-la-redazione-di-piani-di-gestione-forestale-aziendale-Versione-approvata-condeterminazione-del-Dirigente-del-Servizio-Foreste-e-fauna-n.-556-del-19-dicembre-2016 PAT (2024). Carta Geologica della Provincia Autonoma di Trento https://www.provincia.tn.it/News/Approfondimenti/Carta-Geologica-della-ProvinciaAutonoma-di-Trento. Last consultation 18/08/2025. PAT (2025a). Le foreste in Trentino. Servizio Foreste e Servizio Faunistico. https://forestefauna.provincia.tn.it/Foreste/Foreste-in-Trentino/Le-foreste-in-Trentino PAT (2025b). LiDAR - Soleggiamenti DTM - PAT 2006 / 2009. Geocatalogo. https://siat.provincia.tn.it/geonetwork/srv/ita/catalog.search#/metadata/p_TN:5c773de5e3a7-4024-ab4c-ea313d3ba901. Last consultation 18/08/2025. PAT (2025c). Pericolosità Frane: https://siat.provincia.tn.it/geonetwork/srv/ita/catalog.search#/metadata/p_TN:c8b19872aab0-40ce-aeac-6bd5b916ae83; Geocatalogo, Pericolosità Crolli: 113 https://siat.provincia.tn.it/geonetwork/srv/ita/catalog.search#/metadata/p_TN:f26dfcc1-ffc949ca-bc59-e6958e01cc4f. Geocatalogo. Last consultation 18/08/2025. Pearson, R. G., & Dawson, T. P. (2003). Predicting the impacts of climate change on the distribution of species: Are bioclimate envelope models useful? Global Ecology and Biogeography, 12(5), 361–371. https://doi.org/10.1046/J.1466-822X.2003.00042.X Pecchi, M., Marchi, M., Burton, V., Giannetti, F., Moriondo, M., Bernetti, I., Bindi, M., & Chirici, G. (2019). Species distribution modelling to support forest management. A literature review. Ecological Modelling, 411. https://doi.org/10.1016/j.ecolmodel.2019.108817 Pecchi, M., Marchi, M., Moriondo, M., Forzieri, G., Ammoniaci, M., Bernetti, I., Bindi, M., & Chirici, G. (2020). Potential impact of climate change on the forest coverage and the spatial distribution of 19 key forest tree species in italy under RCP4.5 IPCC trajectory for 2050s. Forests, 11(9). https://doi.org/10.3390/F11090934 Pérez Chaves, P., Ruokolainen, K., & Tuomisto, H. (2018). Using remote sensing to model tree species distribution in Peruvian lowland Amazonia. Biotropica, 50(5), 758–767. https://doi.org/10.1111/btp.12597 Peterson, A. T., & Soberón, J. (2012). Species distribution modeling and ecological niche modeling: Getting the Concepts Right. Natureza a Conservacao, 10(2), 102–107. https://doi.org/10.4322/natcon.2012.019 Pitelka, L. F., & Group, P. M. W. (1997). Plant Migration and Climate Change: A more realistic portrait of plant migration is essential to predicting biological responses to global warming in a world drastically altered by human activity. American Scientist, 85(5), 464–473. http://www.jstor.org/stable/27856854 Pottier, J., Dubuis, A., Pellissier, L., Maiorano, L., Rossier, L., Randin, C. F., Vittoz, P., & Guisan, A. (2013). The accuracy of plant assemblage prediction from species distribution models varies along environmental gradients. Global Ecology and Biogeography, 22(1), 52–63. https://doi.org/10.1111/j.1466-8238.2012.00790.x Prosser, F., Bertolli, A., Festi, F., & Perazza, G. (2019). Flora del Trentino (Edizioni Osiride, A c. Di). Fondazione Museo Civico di Rovereto. Puth, M. T., Neuhäuser, M., & Ruxton, G. D. (2015). Effective use of Spearman’s and Kendall’s correlation coefficients forassociation between two measured traits. In Animal Behaviour (Vol. 102, pagg. 77–84). Academic Press. https://doi.org/10.1016/j.anbehav.2015.01.010 QGIS.org (2025). QGIS Geographic Information System. QGIS Association. http://www.qgis.org Raveendra, H. A. C., Picco, L., & Pirotti, F. (2024). Using Earth observation to support the predictive accuracy of species distribution models in ecological restoration: a case study of Poland beavers (Castor fiber). International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences - ISPRS Archives, 48(3), 63–70. https://doi.org/10.5194/isprs-archives-XLVIII-3-2024-63-2024 114 R Core Team (2024). _R: A Language and Environment for Statistical Computing_. R Foundation for Statistical Computing, Vienna, Austria. <https://www.R-project.org/>. Rota, F., Scherrer, D., Bergamini, A., Price, B., Walthert, L., & Baltensweiler, A. (2024). Unravelling the impact of soil data quality on species distribution models of temperate forest woody plants. Science of the Total Environment, 944. https://doi.org/10.1016/j.scitotenv.2024.173719 Sanjit, L., & Bhatt, D. (2005). How relevant are the concepts of species diversity and species richness? Journal of Biosciences, 30, 557–560. Schmitt, S., Pouteau, R., Justeau, D., de Boissieu, F., & Birnbaum, P. (2017). ssdm: An r package to predict distribution of species richness and composition based on stacked species distribution models. Methods in Ecology and Evolution, 8(12), 1795–1803. https://doi.org/10.1111/2041-210X.12841 Seim, A., Marquer, L., Bisson, U., Hofmann, J., Herzig, F., Kontic, R., Lechterbeck, J., Muigg, B., Neyses-Eiden, M., Rzepecki, A., Rösch, M., Walder, F., Weidemüller, J., & Tegel, W. (2022). Historical Spruce Abundance in Central Europe: A Combined Dendrochronological and Palynological Approach. Frontiers in Ecology and Evolution, 10. https://doi.org/10.3389/fevo.2022.909453 Seneviratne, S.I., X. Zhang, M. Adnan, W. Badi, C. Dereczynski, A. Di Luca, S. Ghosh, I. Iskandar, J. Kossin, S. Lewis, F. Otto, I. Pinto, M. Satoh, S.M. Vicente-Serrano, M. Wehner, and B. Zhou (2021). Weather and Climate Extreme Events in a Changing Climate. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Masson-Delmotte, V., P. Zhai, A. Pirani, S.L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M.I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J.B.R. Matthews, T.K. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, pp. 1513–1766, doi:10.1017/9781009157896.013. Soberón, J., & Peterson, A. T. (2005). Interpretation of Models of Foundamental Ecological Niches and Species’ Distributional Areas. Biodiversity Informatics, 2, 1–10. https://doi.org/10.17161/bi.v2i0.4 Takolander, A., Hickler, T., Meller, L., & Cabeza, M. (2019). Comparing future shifts in tree species distributions across Europe projected by statistical and dynamic process-based models. Regional Environmental Change, 19(1), 251–266. https://doi.org/10.1007/s10113-018-1403x Thakur, S., Dhyani, R., Negi, V. S., & Bhatt, I. D. (2022). Water–energy, climate, and habitat heterogeneity mutually drives spatial pattern of tree species richness in the Indian Western Himalaya. Frontiers in Forests and Global Change, 5. https://doi.org/10.3389/ffgc.2022.1022082 115 Thuiller, W., Lavorel, S., Araújo, M. B., Sykes, M. T., & Prentice, I. C. (2005). Climate change threats to plant diversity in Europe. Proceedings of the National Academy of Sciences of the United States of America, 102(23), 8245–8250. https://doi.org/10.1073/pnas.0409902102 Thurm, E. A., Hernandez, L., Baltensweiler, A., Ayan, S., Rasztovits, E., Bielak, K., Zlatanov, T. M., Hladnik, D., Balic, B., Freudenschuss, A., Büchsenmeister, R., & Falk, W. (2018). Alternative tree species under climate warming in managed European forests. Forest Ecology and Management, 430, 485–497. https://doi.org/10.1016/j.foreco.2018.08.028 Vacek, Z., Vacek, S., & Cukor, J. (2023). European forests under global climate change: Review of tree growth processes, crises and management strategies. Journal of Environmental Management, 332. https://doi.org/10.1016/j.jenvman.2023.117353 Wickham, H. (2017). ggplot2: Elegant Graphics for Data Analysis. In Journal of Statistical Software (2; 2a ed., Vol. 77, Numero Book Review 2). Springer-Verlag. https://doi.org/10.18637/jss.v077.b02 Williams, M. I., & Dumroese, R. K. (2013). Preparing for climate change: Forestry and assisted migration. In Journal of Forestry (Vol. 111, Numero 4, pagg. 287–297). https://doi.org/10.5849/jof.13-016 Xu, W., Luo, D., Peterson, K., Zhao, Y., Yu, Y., Ye, Z., Sun, J., Yan, K., & Wang, T. (2025). Advancements in ecological niche models for forest adaptation to climate change: a comprehensive review. Biological Reviews. https://doi.org/10.1111/brv.70023 Zelený D, Lin T. (2021). heatload. Analysis of community ecology data in R. https://anadatr.davidzeleny.net/doku.php/en:customized_functions:heatload Zimmermann, N. E., Jandl, R., Hanewinkel, M., Kunstler, G., Klling, C., Gasparini, P., Breznikar, A., Meier, E. S., Normand, S., Ulmer, U., Gschwandtner, T., Veit, H., Naumann, M., Falk, W., Mellert, K., Rizzo, M., Skudnik, M., & Psomas, A. (2013). Potential Future Ranges of Tree Species in the Alps. In Management Strategies to Adapt Alpine Space Forests to Climate Change Risks (pagg. 37–48). InTech. https://doi.org/10.5772/56279 116 Acknowledgments Part of this work – the available HPC processing – was supported within the Agritech National Research Center and received funding from the European Union NextGenerationEU (PIANO NAZIONALE DI RIPRESA E RESILIENZA (PNRR)— MISSIONE 4 COMPONENTE 2, INVESTIMENTO 1.4—D.D. 1032 17/06/2022, CN00000022)—WP 4.3.