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 provide a range of ecosystem services essential for human wellbeing. In a changing climate, forest management

According to the provided evidence, many uncertainties are still masked under the predictions generally provided by researchers in their studies, with climate as one of the main drivers. At the current time the SDM technique has been successfully used in Italy for some occurring tree species using National Forest Inventory data [27] or wide-range projections and broad spatial distribution data and forest categories [53], or spatial analysis on species richness [54]. However an extensive study on the whole country for a wide range of tree species and evaluating the modelling uncertainties is still missing The aim of this paper is to evaluate the uncertainties behind an SDM procedure in the Mediterranean environment (Italy) to support future SFM strategies. In this work, several future scenarios for 19 species, among the main forest tree species in Italy, were realized using six GCMs and one RCM, quantifying the discrepancies between them and within species when different climatic data are used. Suitability maps were obtained for Italy to provide indications to forest planners regarding the possible consequence and impact of climate change in Italian forest systems. Then adaptive forest management strategies were proposed dealing with potential impacts of climate change and uncertainties detected behind the modelling efforts.

Spatial Data and Climatic Scenarios
Forest inventory plots represent one of the main input data for SDM procedures, given their ability to provide tree-level information which allow a refinement of modelling steps. Among the 263 tree species detected in the framework of the last available national forest inventory (INFC 2005) 19 forest tree species were considered in this study and selected as the most interesting and relevant for Italy under economic, ecological, and aesthetic aspects. Their ecological requirements were previously studied by Pecchi et al. [25]. INFC 2005 was based on a three-phase sampling procedure resulting in a total of 7272 sampling plots, spatially distributed according to a probabilistic sampling scheme [55] and with associated data for 230,874 trees measured in the field [56]. In this framework, statistical inferences on the realized ecological niche of the 19 considered tree species was possible due to the probabilistic sampling scheme.
In order to derive the climatic niche of target species and to project its spatial distribution into the future conditions, current climate data (1981-2010 normal period) were firstly retrieved from the downscaled E-OBS climatological maps. This dataset is available for the whole Italy at 1 km of spatial resolution as a result of a downscaling procedure [46,57]. Such data were then used to generate the set of 19 Worldclim's bioclimatic variables to be used as predictors in SDM. This set is format by a series of biological important variables that better describe the annual and seasonality trends and the extreme and limiting factors [58]. These variables are generated using dismo, a package available for R statistical language [59] using the bioclim function. This step was done to compare the current climate condition with six GCMs we downloaded from the WorldClim website with 30 arc-sec of spatial resolution. The selected GCMs are those elaborated by the fourth version of Community Climate System (CCSM) here and for the following models CC, the Hadley Centre Global Environment Model version 2 family (HADGEM2 2-AO, 2-CC, 2-ES), respectively, HD, HE, and HG, the Max Planck Institute for Meteorology Earth System Model (MPI-ESM-LR) hereafter MP and the Meteorological Research Institute climate model (MRI-CGCM3) MG. To avoid potential biases that originated from different climate data sources (i.e., WorldClim portal and E-OBS data), the WorldClim future projections were recalculated as anomalies from the 1961-1990 climatic normal period, currently distributed as WorldClim version 1.4 [60,61]. Once anomalies were calculated, these were added to the same climatic normal period we obtained from E-OBS for Italy , using spatial reprojection to realign the two grids. An additional climate dataset was then added to this study and provided by the Institute of Bio-Economy (IBE) of Italian National Research Council (CNR), representing the RCM we used in this study. The RCM model is here represented by the output of COSMO-CLM climate model hereafter, COSMO, the climate version of operational weather forecast model COSMO-LM, developed by the German weather service [62]. This RCM was selected for its acknowledged ability to characterize the Italian climate conditions [63]. All climatic scenarios were referred to RCP 4.5 of AR5 for 2050s.

Species Distribution Ensemble Modelling
According to the existing literature, the ensemble forecasting model from different SDM techniques is recognized as the most powerful, stable, and well-referenced method to analyze the potential impact of climate change on tree species [31,64]. An ensemble (or sometimes consensus) modelling is based on the idea that each different modelling output represents a possible state of the real distribution. With this technique, single-model projections are combined into a final surface where the predictions are averaged. In this paper, the ensemble technique was used as predictive method for each of the 19 forest tree species to estimate their potential land suitability under current (i.e., 1981-2010) and future climate conditions (i.e., 2050s, RCP 4.5). The averaging technique was represented by the weighted mean of single model projections using the True Skill Statistic (TSS) indicator [65] calculated with a cross-validation procedure using 75% and 25% for training and testing as weight [42,66]. Furthermore, in order to account for the potential uncertainty that originated from different SDMs, nine algorithms were used for modelling tree species distributions. Fifty replications were performed for each algorithm for a total of 450 single-model projections for each investigated species. The algorithms implemented here include general linear model (GLM), generalized additive model (GAM), classification tree analysis (CTA), artificial neural network (ANN), flexible discriminant analysis (FDA), multivariate adaptive spline (MARS), random forest (RF), and maximum entropy (MAXENT). Codes are available in the biomod2 package [67] in the R statistical language [68].
To avoid collinearity problems amongst the predictors, a principal component analysis (PCA) was performed on the complete set of climatic variables [69]. PCA transforms the original predictors in uncorrelated (i.e., orthogonal) features by preserving the whole variability of the analyzed ecological system (i.e., the ecological variability of the Italian environment). The PCA-derived features were then used as input for the SDMs. Among all the NFI points a threshold of 15% for basal area share was used to filter NFI plots to generate presences (i.e., all the plots where the target species was representing more than 14.99% of total basal area) according to a previous investigation [25]. Afterwards, 10 different pseudo absences datasets (PA) with an equal number each of pseudo-absences than presences were generated with the Surface Range Envelope method [70]. Indeed, even if potentially available from the NFI dataset and detectable from tree-level information, the use of all the plots where the species has not been detected as absences can drive the models to biased predictions, even if setting prevalence to 0.5 [27]. The main reason behind this issue is that, in a managed environment, while the presence is objectively defined, the absence can be due to both inhospitable environment or forest management decision (selective logging, forest management, etc.) and no information is available to confirm any of the above-mentioned possibilities in the NFI data. This generated the final dataset composed by 4500 different single-algorithm single-PA prediction for the consensus model calculation. No soil information was added in the model as it was considered almost stable in the considered time period.

Suitability Maps Analysis and Uncertainties Quantification
From each single modelling cycle with species and climate scenario as cyclers, an ensemble map of land suitability was generated reporting the probability of occurrence of a given tree species in each pixel. A total of 133 future Land Suitability maps (LS) were obtained in addition to 19 current distribution LS maps. A difference in suitability values between future and current distribution maps was calculated for each species and used as input data for a further analysis where the connection between combined use of species and GCM/RCM was evaluated. The variability within GCM/RCM was then studied, with the aim of quantifying the climatic uncertainties in our study as well as the most likely effect of climate change in the Italian environment. To achieve this the 133 LS maps were grouped according to the used climatic simulation and, for each group, the maximum LS value for each pixel was calculated. A single map for each climatic scenario was then obtained representing the probability of a specific location (pixel) to be populated in the future (2050s) by at least one of the 19 considered species. These maps were processed using several LS thresholds, ranging between 51% and 90%, used to transform continuous values in binary predictions (1 or 0). Information on changes in the suitable envelope (i.e., all pixels equal to or higher than the threshold) were derived and especially concerning the total number of pixels (i.e., total forested area in the future) and altitudinal/latitudinal shift (i.e., extension/reduction/movement of the suitable envelope) to determine whether a spatial movement of the suitable envelope could be recognized. A simple linear model was then fitted to examine the influence of different thresholds and climate projections: where CM represents the different climate simulation model we used, β 1 and β 2 were the model coefficients of the fitted model, and TH is the threshold (between 51% and 90%) with ε as error term. Finally, after uncertainty assessment, the most influencing climate change scenario (i.e., the projection calculating the higher differences when referred to current situation) was used to study the most potentially dangerous impacts of climate change on the currently forested areas in Italy. Firstly, the raster of the "maximum pixels value" (i.e., the maximum LS value among all the tested tree species for each pixel) was calculated for both current and most variable future scenario. Then all the INFC 2005 inventory plots were superimposed on the raster and the plot-level LS variation extracted and modelled as a function of plot's attributes. Among these the spatial coordinates (latitude, longitude) the altitude, the forest type (i.e., beech forests, silver fir forests), the admixture level (i.e., pure, mixed), the admixture type (i.e., conifer and broadleaves or the opposite), the main species, and the other components of the forest stand obtained from the INFC2005 dataset were used as predictors in a model. Finally, a Tukey test was used to rank the LS change for each species in order to detect those whose climate change might be more dangerous in the framework of the Italian forest system.

Results
The spatial prediction for the 19 investigated forest tree species showed a wide variability between both algorithms and species. Concerning models, the best results were obtained with RF (average value of TSS 0.844 ±0.092) while the worst performances were observed for MAXENT (average value 0.752 ±0.121). TSS values were more variable amongst species ranging between an average value of 0.647 (±0.113) for Pinus pinea and 0.922 (±0.087) for Pinus cembra ( Figure 1).
Forests 2020, 11, x FOR PEER REVIEW 5 of 20 each pixel was calculated. A single map for each climatic scenario was then obtained representing the probability of a specific location (pixel) to be populated in the future (2050s) by at least one of the 19 considered species. These maps were processed using several LS thresholds, ranging between 51% and 90%, used to transform continuous values in binary predictions (1 or 0). Information on changes in the suitable envelope (i.e., all pixels equal to or higher than the threshold) were derived and especially concerning the total number of pixels (i.e., total forested area in the future) and altitudinal/latitudinal shift (i.e., extension/reduction/movement of the suitable envelope) to determine whether a spatial movement of the suitable envelope could be recognized. A simple linear model was then fitted to examine the influence of different thresholds and climate projections: where CM represents the different climate simulation model we used, β1 and β2 were the model coefficients of the fitted model, and TH is the threshold (between 51% and 90%) with ε as error term. Finally, after uncertainty assessment, the most influencing climate change scenario (i.e., the projection calculating the higher differences when referred to current situation) was used to study the most potentially dangerous impacts of climate change on the currently forested areas in Italy. Firstly, the raster of the "maximum pixels value" (i.e., the maximum LS value among all the tested tree species for each pixel) was calculated for both current and most variable future scenario. Then all the INFC 2005 inventory plots were superimposed on the raster and the plot-level LS variation extracted and modelled as a function of plot's attributes. Among these the spatial coordinates (latitude, longitude) the altitude, the forest type (i.e., beech forests, silver fir forests), the admixture level (i.e., pure, mixed), the admixture type (i.e., conifer and broadleaves or the opposite), the main species, and the other components of the forest stand obtained from the INFC2005 dataset were used as predictors in a model. Finally, a Tukey test was used to rank the LS change for each species in order to detect those whose climate change might be more dangerous in the framework of the Italian forest system.

Results
The spatial prediction for the 19 investigated forest tree species showed a wide variability between both algorithms and species. Concerning models, the best results were obtained with RF (average value of TSS 0.844 ±0.092) while the worst performances were observed for MAXENT (average value 0.752 ±0.121). TSS values were more variable amongst species ranging between an average value of 0.647 (±0.113) for Pinus pinea and 0.922 (±0.087) for Pinus cembra ( Figure 1).  When the standard deviation between projection maps was calculated (Figure 2, left) the central part of Italy was acknowledged as the most uncertain, with spatial projections poorly in agreement. The observed geographical pattern was also partially connected to the spatial shape of the Apennines chain between Latium, Tuscany, and Emilia-Romagna regions. Conversely, a general agreement was observed in flat areas such as the Po valley, spatially next to the central Apennines chain and currently characterized by farms, artificial Populus spp. plantations, agroforestry systems, and agricultural lands. According to the PCA analysis the within-species variability was more influential than the within-scenarios variability. Higher eigenvalues were obtained for factors expressing the between-species variability (e.g., COSMO, CC, HE, HD labels in Figure 2) than those obtained between scenarios which stressed the importance of a species-specific SDM approach. Among the climatic scenarios, the COSMO RCM was the most independent with all the GCMs (i.e., CC, HE, HD etc. labels in Figure 2) partially overlapping with some species and sharing the proportion of explained variability. When the standard deviation between projection maps was calculated (Figure 2, left) the central part of Italy was acknowledged as the most uncertain, with spatial projections poorly in agreement. The observed geographical pattern was also partially connected to the spatial shape of the Apennines chain between Latium, Tuscany, and Emilia-Romagna regions. Conversely, a general agreement was observed in flat areas such as the Po valley, spatially next to the central Apennines chain and currently characterized by farms, artificial Populus spp. plantations, agroforestry systems, and agricultural lands. According to the PCA analysis the within-species variability was more influential than the within-scenarios variability. Higher eigenvalues were obtained for factors expressing the between-species variability (e.g., COSMO, CC, HE, HD labels in Figure 2) than those obtained between scenarios which stressed the importance of a species-specific SDM approach. Among the climatic scenarios, the COSMO RCM was the most independent with all the GCMs (i.e., CC, HE, HD etc. labels in Figure 2) partially overlapping with some species and sharing the proportion of explained variability. In agreement with the PCA results, the histogram analyses of "maximum suitability rasters" reflected the COSMO climate scenario as the most divergent from the others and from current climatic conditions (Figure 3). In agreement with the PCA results, the histogram analyses of "maximum suitability rasters" reflected the COSMO climate scenario as the most divergent from the others and from current climatic conditions ( Figure 3).
While all the other GCMs used in this study showed a density plot mainly cumulated on the right side of the image with values of pixels comprised between 900 and 1000, two distinct peaks were found for COSMO, with the most important between values of pixels between 400 and 600, much lower than those observed for the other GCMs as well as the current scenario too.
The results of statistical model we ran on SDM prediction are reported in Table 1. According to this table, the number of pixels for a specific threshold was substantially similar between GCMs and generally higher than the COSMO model. Then the COSMO was also the most important predictor in the model, i.e., the prediction explaining most of the variability of the system. While all the other GCMs used in this study showed a density plot mainly cumulated on the right side of the image with values of pixels comprised between 900 and 1000, two distinct peaks were found for COSMO, with the most important between values of pixels between 400 and 600, much lower than those observed for the other GCMs as well as the current scenario too.
The results of statistical model we ran on SDM prediction are reported in Table 1. According to this table, the number of pixels for a specific threshold was substantially similar between GCMs and generally higher than the COSMO model. Then the COSMO was also the most important predictor in the model, i.e., the prediction explaining most of the variability of the system. Once the COSMO scenario was acknowledged as the most variable and different among the different climate projections, an assessment of LS change along an altitudinal gradient was calculated over the entire country (Figure 4, left) and only forested areas (i.e., the INFC2005 inventory plot, Figure 4 right side). A potential gain in terms of LS was predicted by the ensemble SDM especially at high altitude but only in the case of the whole Italian country. However, this gain was not able to compensate the global loss of LS, variable according to the threshold we used for binary transformation of the maps but comprised between +4% with 500 as threshold and −81% with 900 and both referred to COSMO modelling. Conversely, only a decrease in LS was found on the INFC2005 domain (i.e., forested areas).  Once the COSMO scenario was acknowledged as the most variable and different among the different climate projections, an assessment of LS change along an altitudinal gradient was calculated over the entire country (Figure 4, left) and only forested areas (i.e., the INFC2005 inventory plot, Figure 4 right side). A potential gain in terms of LS was predicted by the ensemble SDM especially at high altitude but only in the case of the whole Italian country. However, this gain was not able to compensate the global loss of LS, variable according to the threshold we used for binary transformation of the maps but comprised between +4% with 500 as threshold and −81% with 900 and both referred to COSMO modelling. Conversely, only a decrease in LS was found on the INFC2005 domain (i.e., forested areas).
In combination with the histogram analysis and the models described above, the use of a threshold for evaluating the total suitable forested area in Italy stressed the elevation as an important driver ( Table 2).
When such changes were modelled as a function of forest stand characteristics, the altitude variable intercepted the higher proportion of explained variance, close to the 45%. Latitude was highly relevant too, with about 35% of the total variance ( Table 3). The forest category was the last relevant predictor (11%) while the total basal area of the stand and admixture type were much less important than the other variables with values of explained variance of 0.3% and 0.4%, respectively. In combination with the histogram analysis and the models described above, the use of a threshold for evaluating the total suitable forested area in Italy stressed the elevation as an important driver ( Table 2).   A large variability between tree species was detected by the multiple comparison test we ran (Tukey HSD test) where the single-species predictions were analyzed in terms of LS change.
According to our models, the laricio pine (Pinus nigra subsp. laricio), Douglas fir (Pseudotsuga menziesii), and arolla pine (Pinus cembra) were characterized by a positive mean value of LS change, indicating a sort of possible expansion for the three species. Such values were +275.62 for laricio pine, +330.99 for Douglas fir, and finally +460. 16 for Arolla pine that represented the highest value among analyzed species. Such values were statistically significant too (alpha < 0.05) and classified as three different groups ("c", "b", and "a" letters) of statistical similarity where the group are represented by the different letters on the right end of the figure (Figure 5). All the remaining species were characterized by negative average values expressing a decrease of LS, sometimes included in a unique group such as cork oak (Quercus suber) and silver fir (Abies alba) whose means were −63.59 and −63.89, respectively (letters "g"). The second group is composed downy oak (Quercus pubescens) with a predicted decrease of −157.11 and holm oak (Quercus ilex) with −158.17 (letters "k"). Norway spruce was the most stable species with a mean value of −44.47, while the worst projection was calculated for the Mediterranean cypress (Cupressus sempervirens) and stone pine (Pinus pinea) with values of loss of −380.76 and −457.08, respectively. All other species were intermediate and comprised between −100 and −150. However, and despite average values which were just indicative, a wide range of uncertainty was clearly detectable and expressed by the wide range of variability. None of the 19 studied species were placed totally below or above the zero indicating that increasing and decreasing LS values were detected for all the species across the whole study area.

Species-Specific Requirements against a Changing Climate
The species-specific ecological requirements of forest tree species are one of the main drivers for ecological modelling. While similar output can be obtained with species sharing the same climatic Figure 5. Potential absolute variation in suitability values of considered forest tree species. The only spatial distribution map realized with COSMO climate scenario was used for this analysis. The spatial variation was calculated using Tukey test. On the x-axis the potential variation in habitat suitability values is reported.

Species-Specific Requirements against a Changing Climate
The species-specific ecological requirements of forest tree species are one of the main drivers for ecological modelling. While similar output can be obtained with species sharing the same climatic envelope (i.e., silver fir and European beech), different projections are instead calculated for species that are highly differentiated (e.g., European beech and holm oak). Even if just one RCP scenario was used in this study, large differences were found between RCMs and GCMs. Our results underly how the uncertainty on climate change projections have a great impact on spatial model simulation. The use of different types of climatic data (GCMs or RCMs) can lead to very different SDM projections and with potential impacts on SFM decisions [36,71,72]. The use of RCMs with respect to GCMs generally leads to better final climate projection and also to a systematic reduction of bias [73]. This aspect can represent a fair improvement especially for mountainous areas where the use of coarse data can only partially capture the effect of orography [74]. The results we obtained also highlighted the difference in the use of GCMs versus RCMs which are probably optimized scenarios for local areas but very complex and whose calculation is time consuming [61,75,76]. Unfortunately, the use of local data is not still very common and ensemble models are lacking in literature [77]. While several GCMs are sometimes used and then averaged, the use of a single average layer causes the loss of variability with no information on the range of all the potential predictions made by the same SDM procedure. For this reason, an uncertainty assessment should be always mandatory when forecasting climate change impacts. Some papers have also introduced the consensus method to assess the uncertainty in different climate scenarios [78,79], but the use of more GCMs, RCMs, and RCP projections seems to be necessary.
Concerning the mathematical structure of SDM, the importance of the quality of data sources is confirmed as well as its relationship with the uncertainty in species occurrence data and the different statistical technique used to predict the species distribution [45]. Uncertainty in species occurrence data can have a negative effect on the accuracy of a model and any possible correction might bring a potential reduction of the total number of records, removing the uncertain or filtering possible outliers [27]. However, this effect can have different impacts on the SDM according to the modelling technique. Even if MAXENT is the most used in scientific literature and acknowledged as able to provide high accuracy despite the use of occurrence data [31,80], this algorithm was the worst in this study. The reasons might be found in the low number of absences we used (i.e., the background points for MAXENT), probably too few to allow the model to work properly [70]. As a consequence a real and powerful SDM should be based on high-quality data, representative of the phenomena and without any prejudice on the modelling algorithm to be used, with the unbiased comparisons as the unique technique to assess their predictive power [42].
The above-mentioned differences between algorithms, climate projections, target species, LS thresholds to be used for binary transformation etc. (i.e., modelling uncertainties), might heavily impact the operational use of SDM as a decision support system to support strategic sustainable forest management. One of the main uses of SDM is the possibility to identify candidate tree species (genotypes) and provenance types (genotyping) which may be more adapted to future climate conditions in a specific area [9,81]. Provenance selection has the potential to support Assisted Migration strategies (AM) and in-situ or ex-situ conservation efforts to improve the resilience of forest systems [3,82]. While AM represents a possible action for a quick response to climate change threats, this should be realized carefully [10,83]. Such action is probably the most expensive, extreme, and potentially dangerous for ecosystems in case of biased SDM. In fact, despite the advantages attributable to this operation, linked to the avoiding of extinction of species and to supporting economic activity such as timber production, there are many potential disadvantages connected with AM operations that are related to a series of biological risks (the maladaptation or the introduction of invasive species or pests and disease) as well as ethical problems that are connected to the different points of view with respect to the relationship between nature and human and, therefore, the conflict among anthropocentric and eco-centric positions [84,85]. Consequently, AM must be driven by reliable models, averaging different models and GCMs outputs in a framework of statistical probability. The higher the uncertainty in the modelling steps is, the more dangerous and biased the efforts could be, with the probability of failure which is proportional to the magnitude of disconnection between what is projected and what is likely to occur.

SDM as a Tool for Forest Management Options in the Italian Framework
According to the provided evidence, the altitudinal gradient will play a very important role in Italy determining different patterns of species distributions in future climate conditions. This parameter already influences the shape, structure, and specific composition of forests worldwide with a direct effect on a series of important processes, such as water availability, temperature, and soil properties [51,86,87]. The tendency in altitudinal shift of different organisms, both animal and plant, is often confirmed by many research papers [8,12,88,89] with the altitudinal shift generally occurring at much lower speed than latitudinal [83]. If the velocity of colonization of new areas is too low when compared to expected climate change scenarios, then AM might be planned. In this case most of the studies are focused on the upper elevational limit, sometimes also called the leading edge, while the lower elevational limit or rear edges is less investigated even if it is fundamental to plan adequate conservation scenarios for threatened species [28,89,90]. According to Lenoir et al. [8] an average trend shift of 29 m in upward sense for a decade seems to be a reliable value for forest tree species in southern France considering the variation in optimum climate of species in two different periods, that is 1905-1985 and 1986-2005. A confirmation of this process regarding Italian mountains can be found in Rogora et al. [91], where a progressive thermophilization process of climate and a progressive natural introduction of typical species of lower altitudinal strip both for Alps and Apennine has been detected. According to our results, the altitudinal movement of the forested areas with the worst scenario (COSMO) seemed to be lower and around 18 m per decade, demonstrating a possibility of Italian forest tree species to colonize new lands. In this sense, the higher sensitivity to climate change of pure broadleaf stands is one of the main results of our modelling efforts. This result confirms the recent literature where a general contraction of broadleaf species, especially those species that are adapted to cold and wet conditions, was studied [92,93].
According to the provided results, forest management will play a fundamental role in a changing climate. Silvicultural practices in Italy should be aimed at increasing the species richness and favoring hardwoods currently growing as dominated species under conifer canopy, stimulating the natural regeneration, gene flow, and supporting (spatial) migration processes. The spatial variation we found in our models confirms the results of previous studies that establish for Mediterranean areas a general tendency to a loss in habitat suitability as consequence of decreasing precipitation amount and increase in temperature and in frequency and severity of drought period [53,74,94,95]. The possible consequence of climate change may also be accompanied by an increased wildfire and safety risk. This issue has also been acknowledged in many other research studies and mainly connected to extreme climatic events [96,97]. However, uncertainty assessment has also been detected as fundamental in this case too when predictive models are generated [98].
Our results highlight that only three species seem to be favored by climate change phenomenon: Arolla pine, Douglas fir, and laricio pine. The scientific literature confirms these results. As example, Casalegno et al. [99] indicated an increment of spatial distribution of Arolla pine as consequence of the progressive abandonment of pastures in the Alps. Instead, Douglas fir is indicated as a tolerant species versus drought events [100,101] and its habitat suitability is indicated in increment in Europe in future time periods by Dyderski et al. [102]. In the end, laricio pine is indicated as tolerant to heat and drought events [103]. Considering broadleaf species and the oaks group, our outcome is partially in agreement with the existing knowledge. A negative variation in habitat suitability was also predicted in Perkins et al. [104] while a negative decrease in habitat suitability can be read in Kim et al. [105] for cork oak in the Mediterranean area. Finally, a decrease in habitat suitability during the future period was predicted for holm and turkey oak by Vitale et al. [106] in the same environment we studied. Conversely a contrasting result with existing literature was found for downy oak. In this sense, our outcome highlights a possible decrease in habitat suitability while a potential increment was calculated by Vacchiano and Motta [88]. However, the difference might be attributed to the spatial extent they studied, a small region of Northern Italy where our model predicted an increase too. With attention to other broadleaf species, a possible negative variation in suitability values is expected for European beech that confirms the hypothesis by Noce et al. [53] and especially in the center and south of Apennine. With attention to conifer species such as European larch and Norway spruce, a reduction of habitat suitability is a possible event. A negative variation in habitat suitability of Norway spruce represents a focal point for the high economic value of timber and it reported by different previous works [92,107]. A negative variation in habitat suitability for European larch is also confirmed by Dyderski et al. [102] and Mamet et al. [108]. Additionally, silver fir loss in habitat suitability values is in agreement with Vitasse et al. [109]. Given the economic relevance of this group of species, SFM in Italy should take particular care in their management and supporting local enterprises and avoiding species substitution, maybe using different provenances and genotypes [110,111]. Finally, despite being considered as typical Mediterranean species, a possible decrease in suitability was predicted for species such as Italian cypress and stone pine. This possibility was confirmed in Klein et al. [112] with attention to Italian cypress and in Freire et al. [113] if stone pine is considered. A decrease in suitability values is also a possibility for Aleppo and Maritime pine and finally for black pine as previously reported by Silvério et al. [114] and by Buras and Menzel for black pine [74]. All the cited literature reports as main causes of the decrease in habitat suitability a high sensitivity towards drought events, an increase of wildfire events, and in the end an increase of pests and pathogens. Even if probably in agreement with literature, the low occurrence across INFC2005 might be the main shortcoming of our model for these tree species, owing to a possible underestimation of their potential ecological niche which will be the real niche responding to climate. In this framework, only monitoring efforts and provenance trials will support the solution of the issue in the next decades and allowing models to consider phenotypic plasticity [9].

Conclusions
Climate change will probably affect the spatial distribution of forest tree species worldwide and many research groups are currently working to adapt GCMs to local contexts. Anyway, the uncertainty is still wide. Many factors are involved with physical and anthropogenic processes on one hand and all the possible adaptive processes of forest systems to deal with climate change scenarios on the other, which are only partially known in a long-term period. With this study, an initial framework of the possible consequences of climate change phenomenon in Italian forest was proposed under the Fifth Assessment Report projections, trying to understand the different dynamics between different variables and not merely describing the potential expected species geographical shift. While any model can be built with any data coming from different sources, a real uncertainty assessment is fundamental to support useful and effective SFM strategies. Dealing with uncertainties and working with self-updating procedures seems to be the main path to address climate change effects properly, mitigating the negative effects and maintaining the delivery of ecosystems services from forests. Anyway, only monitoring networks and species-specific analysis will be able to certify or confute this tendency. Such new data will be fundamental to test current SDM and adjust projections properly. Additional results may be then provided using the new climate change pathways provided by the new IPCC projection in the Sixth Assessment Report. Funding: This study was partially supported by the PhD grant provided by the University of Florence to Matteo Pecchi. The MDPI Article Processing Charge fee was fully covered by Maurizio Marchi using the voucher he obtained from forests journal as one of the four best reviewers selected by the journal with the "2019 Outstanding Reviewer Award".