Likelihood of changes in forest species suitability, distribution, and diversity under future climate: The case of Southern Europe

Abstract Forest conservation strategies and plans can be unsuccessful if the new habitat conditions determined by climate change are not considered. Our work aims at investigating the likelihood of future suitability, distribution and diversity for some common European forest species under the projected changes in climate, focusing on Southern Europe. We combine an Ensemble Platform for Species Distribution Models (SDMs) to five Global Circulation Models (GCMs) driven by two Representative Concentration Pathways (RCPs), to produce maps of future climate‐driven habitat suitability for ten categories of forest species and two time horizons. For each forest category and time horizon, ten maps of future distribution (5 GCMs by 2 RCPs) are thus combined in a single suitability map supplied with information about the “likelihood” adopting the IPCC terminology based on consensus among projections. Then, the statistical significance of spatially aggregated changes in forest composition at local and regional level is analyzed. Finally, we discuss the importance, among SDMs, that environmental predictors seem to have in influencing forest distribution. Future impacts of climate change appear to be diversified across forest categories. A strong change in forest regional distribution and local diversity is projected to take place, as some forest categories will find more suitable conditions in previously unsuitable locations, while for other categories the same new conditions will become less suited. A decrease in species diversity is projected in most of the area, with Alpine region showing the potentiality to become a refuge for species migration.

preserving ecosystems against both natural hazards and human threats, in particular due to environmental resources' exploitation and land-use changes. Therefore, in delineating these strategies and plans, the concepts of biogeography and expected shifts in species habitats, as well as the modification in species regional distribution and local diversity under climate projections, become critical (Araújo, Alagador, Cabeza, Nogués-Bravo, & Thuiller, 2011;Araujo, Cabeza, Thuiller, Hannah, & Williams, 2004;Bakkenes, Alkemade, Ihle, Leemans, & Latour, 2002;Guisan et al., 2013;Kelly & Goulden, 2008;Serra-Diaz, Scheller, Syphard, & Franklin, 2015). For instance, while a recent study on forest presence hot spots demonstrated the local value of protected areas in southern Europe (Noce, Collalti, Valentini, & Santini, 2016), these areas do not guarantee the same benefits under future environmental scenarios .
It is therefore crucial to recognize the role of biogeography and of the factors affecting ecosystems' composition and biodiversity at scale of species ranges over time (Franklin, 2016).
Forest ecosystems are strategic for biodiversity conservation; on a centennial scale, they have evolved their resilience and adaptation capability to disturbances (e.g., droughts, fires, windstorms, pests, diseases, and invasive species), including migration as an option (Aitken, Yeaman, Holliday, Wang, & Curtis-McLane, 2008).
Under global change and, thus, with an altered intrinsic vulnerability due to modified average environmental conditions, forest communities have to face an additional challenge: coping with a quickly increasing variability of extreme events and disturbances (Seidl, Spies, Peterson, Stephens, & Hicke, 2016) as well as novel perturbations (e.g., new diseases; Pautasso, Schlegel, & Holdenrieder, 2015). Such complex transformations are occurring too fast to be accompanied by both evolutionary adaptation (Dale et al., 2001;Lindner et al., 2010;Trumbore, Brando, & Hartmann, 2015) and migration processes, although the latter seems the best option for species (Corlett & Westcott, 2013). It is widely demonstrated that the geographic distribution of forest species is strictly correlated with medium-to long-term climate conditions (Araújo & Pearson, 2005;Park Williams et al., 2012;Schimper & Fisher, 1903), topographic factors (Bellingham & Tanner, 2000) and, especially in the Mediterranean region, human activities (Barbero, Bonin, Loisel, & Quézel, 1990). These elements, when combined, realize the concept of niche (Guisan & Zimmermann, 2000).
Correlative Species Distribution Models (SDMs) that are also known as bioclimatic envelope models, correlative ecological niche models, or habitat suitability models, explore the relationships and the equilibrium between the geographical distribution of species and a set of environmental variables (Austin, 2002;Guisan & Zimmermann, 2000;Naimi & Araújo, 2016;Peterson et al., 2011).
SDMs, in conjunction with Geographic Information Systems (GIS) tools, are promising research tools to map and predict the potential spread of endemic or invasive species in the past and the future, respectively (Franklin, 2010). Thanks to these tools, scientists are in fact enabled to inform decision makers and stakeholders on how to formulate or prioritize biodiversity and biogeography conservation plans and to contribute to the maintenance of multiple services provided by forests, including the mitigation of climate change and its impacts (Franklin, 2013;Gama, Crespo, Dolbeth, & Anastácio, 2015).
Nevertheless, SDM approach has some limitations in particular related to unavoidable assumptions and uncertainties (Guisan & Thuiller, 2005;Watling et al., 2015). The most important element concerns the selection of environmental explanatory factors (or predictors) (Lexer & Hönninger, 2004). Often this choice is constrained by the availability of datasets (that also have their own uncertainty as stated in Bedia, Herrera, & Gutierrez, 2013), or by the suspicion of redundancy among predictors that seem correlated (collinear). Unfortunately, any mechanistic and ecological understanding of the "true" predictors still lacks of interpretability (Dormann et al., 2013). Even in case of dealing with collinearity, adopting one or another method presents different uncertainties (Dormann, Purschke, Márquez, Lautenbach, & Schröder, 2008). Moreover, the assumption that the relationships between predictors and species presence/absence assessed for the historical period will maintain in the future could be really uncertain (Gavin et al., 2014;Hijmans & Graham, 2006). Moreover, a modeled species is rarely observed in its full climate space, and usually model performances are tested just in a "restricted" climate space (Hannemann, Willis, & Macias-Fauria, 2016) and within boundaries smaller than physiological limits of species (Loehle & LeBlanc, 1996), especially in regional studies; this is also due to the limited availability of harmonized species distribution data. Another major limitation is that SDMs do not routinely consider relevant population and dispersal dynamics or intraspecific variation in climatic tolerances (Holt, 2009;Zurell et al., 2016). Summarizing, results are widely dependent on: (i) the reliability and accuracy of species occurrence data; (ii) the significance of the environmental variables selected; (iii) the quality of related data; and (iv) the parameterization or configuration of the applied models (Chakraborty et al., 2016;Nenzén & Araújo, 2011;Thuiller, 2003;Thuiller, Lafourcade, Engler, & Araújo, 2009). Given that all above elements cause a large variability in the predictions (Cheaib et al., 2012;Pearson et al., 2006;Thuiller et al., 2004), the Ensemble Forecasting approach has been developed and widely adopted (Araujo & New, 2007;Heikkinen et al., 2006;Komac, Esteban, Trapero, & Caritg, 2016;Marmion, Parviainen, Luoto, Heikkinen, & Thuiller, 2009). This approach combines individual SDM predictions to provide consensus predictions (Capinha & Anastácio, 2011), enabling more robust evaluations, that is, addressing the uncertainty related to SDMs.
After assessing the accuracy of the Ensemble Forecasting-SDM approach in reproducing the historical species distributions, results for future projections mainly depend on the updated values of the considered dynamic predictors (i.e., changing at the time scale of simulations). These are usually bioclimatic variables calculated from the outputs of climate model simulations that, in turn, are driven by scenarios on greenhouse gas (GHG) emissions or concentrations.
Considering a single projection is, however, not recommended, indeed, as highlighted by Lindner et al. (2014), the scientific community cannot still accurately forecast GHG emissions and the ways in which the climate will evolve.  (Friedlingstein et al., 2014;Taylor, Stouffer, & Meehl, 2012) due to the models' physics, initialization and/or configurations, increased by the consideration of four different representative concentration pathways (RCPs) (van Vuuren et al., 2011). Such variability in climate projections needs to be considered, quantified, and well managed, especially when used in concatenated impacts evaluations, so to explore a broad range of possible developments (Cheaib et al., 2012).
Further studies (Goberville, Beaugrand, Hautekèete, Piquot, & Luczak, 2015;Keenan, Maria Serra, Lloret, Ninyerola, & Sabaté, 2011;Lindner et al., 2008;Wang, Campbell, O'Neill, & Aitken, 2012) concerning climate change and impacts on European forests focused the attention on the large uncertainty and lack of knowledge related to future distribution scenarios. It is thus clear that climate change impact evaluations could be used for planning and decisions only after assessing the uncertainty or, to more directly inform stakeholders and policy makers, after treating the uncertainty by quantifying the likelihood of outcomes based on the agreement among multiple simulations.
In this work, we adopt the Ensemble Forecasting-SDM approach to predict the possible impacts of climate change in terms of geographic range shifts, over medium and long term, for ten forest categories (groups of species). We focus our attention on southern Europe, and the Mediterranean Basin in particular, as they are among the world's major areas for plant biodiversity and endemism (Medail & Quezel, 1997, 1999Myers, Mittermeier, Mittermeier, daFonseca, & Kent, 2000;Underwood, Viers, Klausmeyer, Cox, & Shaw, 2009). These regions are dominated by a Warm Temperate climate that ranges from dry and warm to hot summers (Alessandri et al., 2014;Kottek, Grieser, Beck, Rudolf, & Rubel, 2006), and they are projected, with a high degree of consistency among different projections, as hot spots of climate change under IPCC-AR4 scenarios (Giorgi, 2006). In these areas, raising temperature and decreasing summer precipitations will lead to an increase in summer droughts (Giorgi & Lionello, 2008;Mariotti et al., 2008). In terms of average climate, the robust assessment (i.e., derived by probabilistic ensemble evaluation) in Alessandri et al. (2014) suggests that this climate is expected to experience a northward shift under IPCC-AR5 intermediate (RCP4.5) emission scenario, leaving space for more arid conditions as also evinced in Santini and di Paola (2015). Previous impact studies like those in Santini, Collalti, and Valentini (2014) also show the consensus in predicting an increase of water scarcity and fire disturbances in this region, with negative consequences for forest ecosystems.
Other studies (Duputié, Zimmermann, & Chuine, 2014;Zimmermann, Jandl et al., 2013), and the JRC "Tree Species Maps-Species Habitat Suitability" dataset (http://forest.jrc.ec.europa.eu/ download/data/species-distribution/) in particular, give useful information about the expected forest suitability in the European area. Our work aims to complement these efforts by merging/harmonizing two datasets of forest species presence (Noce et al., 2016), and exploiting the latest projections under CMIP5, guaranteeing robustness in terms of model consensus in predicting future distribution of each forest category. We use the IPCC terminology on likelihood to treat the uncertainty on future outlooks (Mastrandrea et al., 2011). The approach we adopt is what we call a cascade ensemble system, which concatenates the ensemble forecasting approach of SDMs to a sub-ensemble of CMIP5 climate projections.
Further, we investigate, for the whole study area and for its main bio-geographical subregions, the projected future evolution in terms of forest spatial arrangement (distribution) and local species diversity, which is part of the wider biodiversity concept (Boudouresque, 2011).
We then compare and discuss the overall importance, averaged across model, of bioclimatic predictors at species level, corroborated with some ecological and mechanistic considerations based on knowledge of historical habitat distribution. The results presented here can be helpful for further ecological and biodiversity conservation evaluations as well as support of medium-to long-term protection or restoration strategies that account for climate change and its uncertainty range.

| MATERIALS AND METHODS
In the following paragraphs data and methods are described, while

| Species distribution modeling
The potential distribution of the considered tree forest categories under climate change was projected through the "BIOMOD2" package v3.3-7 (https://cran.r-project.org/web/packages/biomod2/index. html; Thuiller, Georges, Engler, & Breiner, 2016) implemented in R (R Core Team, 2015). The BIOMOD2 platform uses individual species representation and a community-based approach, enabling ensembles of models to be evaluated (calibration phase) with the historical conditions (represented by environmental predictors) and then re-applied to project the potential spatial distributions (habitat suitability driven by climate) of species along future time horizons (Naimi & Araújo, 2016).
Ten models are embedded in BIOMOD2 and include the following: two regression methods, which build linear or nonlinear relationships between species occurrence and their environmen- BIOMOD2 also offers a set of rules and metrics for SDM evaluations, based on the confusion matrix and its elaborations (Miller, 2010) into three metrics: Area Under the Curve of the Receiver Operating Characteristic plot (ROC), Cohen's K and True Skills Statistics (TSS) (Duque-Lazo et al., 2016). ROC is a threshold-independent model evaluation indicator (Franklin, 2010), and it is also independent of prevalence (i.e., the frequency of occurrence) of target species. It plots the commission error (1-specificity; false positives) against omission error (sensitivity; true positives). ROC ranges between 0.5-1, where 1 represents a perfect discrimination between presence and absence, and 0.5 represents a random fit. Both Cohen's K and TSS are threshold-dependent measures of model accuracy. They both range from −1 to +1, with +1 indicating perfect agreement between predictions and observations, and 0 or less indicating an agreement no better than a random classification (Zhang et al., 2015). Cohen's K has been criticized as being strongly influenced by the species prevalence in the data, and the TSS has been introduced to solve this problem (Allouche, Tsoar, & Kadmon, 2006).
Regardless of the evaluation metrics, model should always be evaluated by means of independent data or adopting data splitting procedures. In BIOMOD2, a user-defined proportion of the original data can be used for training the models, while the remaining proportion is used for model evaluation. The package also allows conducting n data splitting runs, providing an n-fold cross-validation procedure. Only runs that meet criteria in terms of evaluation metrics are included in the ensemble, and weighted according to their performances before to generate the final binary (presence/absence) map.
In addition to model evaluation, the importance of selected and tested predictors can be analyzed. Once the models are trained (i.e., calibrated), a standard prediction is made. Afterwards, one of the variables is randomized, and a new prediction is made. The correlation score between the new prediction and the standard prediction is considered to give an estimation of the variable importance in the model. A good correlation score between the two predictions suggests that the randomized variable has little importance. On the contrary, a low correlation means a significant difference in the predictions, making that F I G U R E 1 Graphical overview of methods variable important for the model. Importance is expressed as 1 minus correlation and converted to percentage for easier interpretation.
SDMs predict species occurrences in space and time based on simple to complex relationships between presence-only or presence-absence point records of species and environmental variables.
Aiming to explore the impacts of climate change on the geographic distribution of forest species in southern Europe, we investigated a large set of available spatial datasets on forest distribution (Trombik & Hlásny, 2013) and produced a database that expressed the fraction of presence for ten forest categories (i.e., grouping several species) as described in Noce et al. (2016). This database was created through forest category harmonization and spatial overlay procedures by merging two existing Pan-European datasets: (1) the European Forest Institute (EFI) "Tree species maps for Europe" (Brus et al., 2012) and (2) the Joint Research Centre (JRC) "Novel Maps for Forest Tree Species in Europe" (Köble & Seufert, 2001). The forest category's distribution layers within this Source Database (hereafter SDb), originally at 30 arc-sec resolution, have been upscaled to 0.25°, and the fractional presence value has been averaged within the new coarser cell by means of ESRI ArcGIS 10.1.
The 0.25° resolution was selected after several considerations.
The most important one, as stated by Kriticos and Leriche (2010), is that maintaining a fine resolution without compatible resolution of climatic data leads to inaccurate projections and redundancy in data sampling. It should be taken into account that the original resolution of climate projections data by General Circulation Models (GCMs) from CMIP5 used here (see below) ranges from around 1° to 2.5° (see, e.g., Scoccimarro, Gualdi, Bellucci, Zampieri, & Navarra, 2013).
Furthermore, we conducted a regional-scale study and, as demonstrated by Guisan et al. (2007), very high resolution would not improve the accuracy of the SDMs. Conversely, too coarse resolution (approximately 0.5° or more) can identify inappropriate regions (Seo, Thorne, Hannah, & Thuiller, 2009). The adopted resolution was a compromise among the resolutions of the input datasets.
We subsequently created the Species Occurrence Database (SODb) by transforming continuous (percentage of presence) data in SDb into the presence/absence data as needed by SDMs. We assumed the values above 2% could be defined as "presence," while lower values were considered as "absence." This threshold allows excluding simulations of small and isolated populations that may represent relicts of past climate change (see Petit, Hampe, & Cheddadi, 2005) or small plantations. The ten forest categories considered in our work, with the associated (group of) species, are described in Table 1. Thus, SODb contained ten layers, each comprising 4,352 presence/absence locations (map units or pixels). Hereafter, the forest category name will be used instead of genus names.
F I G U R E 2 Envelope of study area (black boundaries) with forestlands and shrublands according to FAO Global Land Cover SHARE

| Environmental predictors and SDM calibration
Two categories of environmental predictors (EPs) have been considered for the calibration phase of SDMs. Nine (9) topographic predictors (TPs) and nineteen (19) 1960-1990Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), resampled up to 10 arc-min (10′, approx. 0.16°) resolution. They were further upscaled in ESRI ArcGIS to 0.25° resolution with the help of bilinear method to fit the resolution of forest category presence data. Similarly, maps of TPs were produced starting from the Digital Elevation Model, as processed from the Shuttle Radar Topographic Mission data (SRTM; Farr et al., 2007) to create other topographic derivatives (slope, aspect) and descriptive statistics (mean, minimum, and maximum values). The 90-m resolution of the original dataset was resampled to 0.25° resolution. Statistics for elevation and slope were calculated for the 90-m resolution data within the new 0.25° pixel.
Several BPs can be highly correlated (or collinear), as they come from the same variables (temperature or precipitation) just operated at different scales (year, season, month), or even combined one another; however, the issue of whether or how treating collinearity in SDMs is still highly debated and unsolved (Dormann et al., 2013). First, collinearity is surely most important for linear models (De Veaux & Ungar, 1994) that are a very limited part of BIOMOD2 package mostly made of nonlinear and/or machine-learning algorithms. Furthermore, collinearity is more noteworthy when the main goal is inferring the influence of predictors over occurrence, that is, for process understanding and/or interpretation purposes, while it seems not impacting predictions (De Veaux & Ungar, 1994).
However, also in case of predictions, collinearity can be ignored only if the SDMs are applied over the same space and time, that is, if the collinearity structure remains constant (Dormann et al., 2013). While in our study the spatial domain is the same between calibration and projections, we also verified, by quantifying the collinearity among variables through the most common metrics (pairwise correlation coefficient, r, and the variance inflation factor, VIF; Naimi & Araújo, 2016), that the collinearity structure is not statistically different, at T A B L E 1 Name of forest categories and related species Although mechanistic ecological knowledge for variable selection as well as more objective data reduction methods, such as principal component and factor analysis, can be used to reduce multicollinearity, this is often at the expense of interpretability (Miller, 2010).
Given all that above, we agree with Mateo, Felicísimo, Pottier, Guisan, and Muñoz (2012) that with no a priori reason for removing some variables, we can keep all of them for the analyses. The consequent risk of overfitting, likely for species with low number of occurrences, does not apply in our case, reducing the possibility of predicting more around the know presences (Mateo et al., 2012).  Zhang et al. (2015).
SDM ensemble projections in terms of probability of presence for each forest category were weightily averaged according to their performances and converted into a binary (presence/absence) prediction map using a threshold that maximized the evaluation statistics considered for this purpose, that is, ROC.
For these final maps, the analysis of performances was conducted via some metrics based on the confusion matrix.

| Future projections and likelihood quantification
To project future forest categories' distribution, the topographic predictors were considered unchanging, while for bioclimatic predictors we used those extracted by five bias-corrected CMIP5 GCMs: GDFL-CM3, HadGEM2-CC, MIROC5, INM-CM4.0, and CSSM4 (for acronyms and description see Scoccimarro et al., 2013) driven by RCP 4.5 and RCP 8. After producing maps of potential spatial distribution (namely climate-driven habitat suitability, hereafter also simply "suitability") of each forest category under each CEM, the spread in results due to GCMs and RCP scenarios was addressed by adapting the approach and terminology to treat the uncertainty and communicate the "likelihood" of outcomes, as proposed by IPCC-AR5 (Mastrandrea et al., 2011). If the suitability is predicted only by 1 CEM, the outcome is considered as extremely unlikely; if by 2 or 3 CEMs, as unlikely; if by 4 to 6 CEMs, as about as likely as not; if by 7 or 8 CEMs, as likely; and if by 9 or 10 CEMs, as extremely likely (Table 3).

| Changes in distribution and diversity
Two approaches were adopted to evaluate aggregated changes in terms of regional distribution and local diversity of forest categories, by averaging presence/absence results of the five GCMs for each RCPs and time horizons, called hereafter sub-CEMs.
First, the re-arrangement of forest categories' frequency distribution was assumed a proxy of modifications in large-scale forest composition (regional distribution) across the future new suitable areas for forests.
Second, the number of coexisting forest categories at spatial unit (pixel) level was considered representative of small-scale forest composition (local diversity), and the updated frequency distribution of pixels with different presence (number) of forest categories was analyzed.
In both cases, the statistical significance of the differences in frequency distribution across RCPs and time horizons was assessed through eu/data-and-maps/data/biogeographical-regions-europe-3).
By combing the two approaches, it is possible to evaluate changes in both large-scale distribution and local diversity seeing, for example, if the general forested areas increase/decrease but the large-scale composition remains unaltered, or if the diversity of the spatial units remains "quantitatively" similar under future outlooks but it is due to a mosaic of different forest categories.

| Ensemble SDM performance and predictors' importance
As described in Miller (2010), there are several measures to test the predictive performance of SDMs.

| Future projections
For each forest category, 10 maps (5 GCM by 2 RPCs) were produced and then aggregated into a final map (one for each time horizon) that Concerning the changes in large-scale distribution and diversity, Alps and the Carpathians) these against a spread of QuercusSP in the northern side of study area (Fig. 4)  the currently most frequent mix of three forest categories is substituted by the highest frequency of map units made of two categories (medium term or RCP4.5) up to one category alone in the long-term RCP8.5 (Fig. 6).
Looking at biogeographical region level, some peculiarities emerge, with the Atlantic (Fig. S21), Continental (Fig. S22), and Mediterranean in the long-term RCP8.5 scenario, the areas with two and three categories become very frequent.
The Mediterranean region presents the most affected trends in either long term or RCP8.5 sub-CEMs, significantly different from those under medium term or RCP4.5 sub-CEMs, and especially from historical conditions. This fact suggests accelerated dynamics of suitability loss for the considered forest categories and a re-arrangement in distribution. The region is projected at risk of progressive disappearing of Betula, Larix and Picea; a reduction of more than a half of Fagus, PinusSylv and QuercusRP; a weaker reduction in Abies, Castanea and PinusPin; and a confirmed dominance of QuercusSP.
The smoothly diversified landscape, which appears largely represented by map units with two forest categories, loses diversity, and the zones with a single forest class seem dominating under all future projections.
A diverging trend is shown by the Alpine region (Fig. S24). After a slight overall increase of forest presence under either RCP4.5 or medium-term sub-CEMs, with significant re-arrangements of landscape forest composition, an even slighter reduction is projected toward the re-establishment of the "historical" presence (but with a noticeable reduction in Betula and Picea) and a spread of the two Quercus categories, although the dominance of Fagus seems to be maintained. In terms of local diversity, no significant difference exists across sub-CEMs, suggesting that, under this preserved "quantitative" diversity, new species seem finding suitable conditions in the Alpine areas.

| DISCUSSION
The findings of this work should be regarded as "habitat suitability" or simply "suitability," referring to favorable habitats for the analyzed forest categories rather than to their expected future spatial distribution under the impacts of climate change, because of three main points.
First of all, the investigated forest categories represent a signifi- Furthermore, to interpret the suitability directly as an indicator of forest presence/absence, species migration should be considered free from ties, and constraining peculiarities of the study area should not be taken into account. Such peculiarities comprise the high degree of fragmentation of habitats owing to intensive anthropic land use, including the dense network of infrastructure barriers, and the topographic limits related to high topographic heterogeneity (Saltré, Duputié, Gaucherel, & Chuine, 2015). In this regard, the capacity of individual species to colonize new favorable areas or to adapt in existing ones depends partially on their competitive and dispersal capabilities (Boulangeat, Gravel, & Thuiller, 2012;Corlett & Westcott, 2013;Schiffers, Bourne, Lavergne, Thuiller, & Travis, 2013;Zhu, Woodall, & Clark, 2012).
Moreover, the identified likely transitions from a forest category to another are "potential" and we can imagine they cannot take place by natural succession in a short timeframe like the next 50-60 years, as expansion rates of species are much slower than the expected rate of warming: to bridge 100 km, an acorn-bearing species like oaks or beech may need 500 years or even more (Saltré et al., 2013), while wind-distributed tree species, like birch or some conifers, may reach much faster expansion rates. For this reason, the transition to a new equilibrium has to be supposed including ecological imbalances, like the appearance of unstable communities dominated by windpropagated, pioneer tree species or even intermittent deforestation of wide landscapes (Bussotti, Pollastrini, Holland, & Brüggemann, 2015).
Although climatic changes are still too rapid to allow a natural rebalancing of forest ecosystems (Corlett & Westcott, 2013), given the time horizons of this work, the species migration dynamics are considered more facilitated to adapt to new climatic conditions than evolutionary ones (Wiens, 2011).
However, the produced maps of future suitability allow summarizing three kinds of impacts on the potential re-distribution of forest categories under alternative climate outlooks: • Minor or no impacts on the historical distribution while new areas become suitable. PinusPin and QuercusSP belong to this group, as they seem to be positively affected by changing conditions, and the expected increasing temperature in particular (Giorgi & Lionello, 2008 (Falk & Hempelmann, 2013;Hickler et al., 2012;Kramer et al., 2010;Saltré et al., 2015).
• Negative impacts on the most part of historical distribution. According to our results, extremely likely and likely suitability for Betula will be guaranteed only in few areas (the Western, Central, and Dinaric Alps) with a strong reduction in the remaining range (particularly in France). Likewise, especially long-term impacts on Castanea are significant and largely widespread, matching with the results of Goberville et al. (2015). Particularly, suitability sensibly decreases in areas (Southern Italy), which are already affected by serious diseases (Vettraino et al., 2009 Meanwhile, the expected spreading of QuercuSP can be explained by the variety and the heterogeneity of this categories, but special attention must be paid to the auto-ecological traits of the species within; in particular the high resilience of these Mediterranean species to drier and warmer conditions (e.g., Quercus ilex L., Quercus pubescens Wild, and Quercus cerris L.) as described in Baldocchi et al. (2010), Morán-López, Poyatos, Llorens, and Sabate (2014), and Vaz et al. (2010).
Recent studies regarding a mountain area in Slovakia (Hlásny et al., 2017) reported the crucial importance of annual temperature and of temperature in the warmest and coldest months for productivity of some Abies species as reported also by Lebourgeois (2007)  Local scale diversity, analyzed in terms of potential presence in the same map unit, shows an appreciable reduction, especially on the long term, throughout the study domain. This alarming trend affects almost all bio-geographical subregions as described in other studies (Sala et al., 2000;Thuiller, Lavorel, Araujo, Sykes, & Prentice, 2005). The Alpine subregion, nevertheless, seems to be less affected, perhaps due to the combined effect of reduction of mountainous species (e.g., Picea abies) and new colonization of mesophilic species (Quercus robur and Q. petraea) that might find more suitable climatic conditions. Similar effects on alpine vascular vegetation have also been described by Pauli, Gottfried, Reiter, Klettner, and Grabherr (2007). Although observed altitudinal shift in the Swiss Alps have been associated more to land abandonment rather than to climate change (Gehrig-Fasel, , our study considers potential suitability, so the space left or not by other land uses does not influence the outcomes. The topic of diversity, even limited to the (group of) species level, is of strategic importance to increase ecosystem productivity because resources are better shared among neighboring species and are thus potentially more available (Loreau et al., 2001). Moreover, Jactel and Brockerhoff (2007), Lebourgeois, Gomez, Pinto, and Mérian (2013), Pretzsch, Schütze, and Uhl (2013) and Zhu et al. (2000) demonstrated that diversity increase resilience to both biotic stresses such as insect pests or diseases and abiotic stresses as droughts, especially for drought-prone areas as southern Europe. This point is crucial given the expected increase of drought frequency for the Mediterranean area (Vicente-Serrano et al., 2014).

| CONCLUSION
In this work, we coupled a set of SDMs with alternative climate projections under different global circulation models and emission scenarios.
We aimed first at investigating the likelihood of future changes in the climate-driven habitat suitability, or simply "suitability," including the modification in (regional) distribution and (local) diversity, for some forest categories in Southern Europe. Then, the cascade ensemble system, allowed to treat together both SDMs and climatic projections' uncertainty, providing probabilistic and thus more robust information about potential future changes in forest habitats under changed environmental conditions.
Results in terms of forest shift, regional composition, and local diversity confirmed how climate change is likely to led to significant modifications in the future, affecting forests with different degrees of magnitude across species and with various levels of uncertainty within species also due to the spatial heterogeneity, the auto-ecological traits, and the adaptation strategies.
While some forest categories will find more suitable conditions in previously unsuitable locations, for other categories, the same new conditions will become less suited. A decrease in local species diversity is projected in most of the area, with Alpine region showing the potentiality to become a refuge for species migration.
Results, flagged with the related likelihood, represent useful and immediate information to be communicated to stakeholders and policy makers for supporting their design of future forest conservation as well as protection strategies and plans.