**3. 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).

**Figure 1.** True Skill Statistic values (TSS) obtained during the cross-validation in the Species Distribution Modelling procedure for each involved algorithm (**left**) and for each species (**right**).

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.

**Figure 2.** (**a**) Spatial pattern of the standard deviation of the LS difference between future (2050s) and current (1981–2010) values for each species and using all the future climate realizations (133 layers) and (**b**) PCA analysis run on the same data (i.e., LS variation at pixels level).

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.

**Figure 3.** Density distribution (histogram) of each "maximum GCMs and RCM" obtained in this study when using the maximum land suitability value for each pixel within the 19 analyzed species.

**Table 1.** Results of the linear model to determine the most important climate scenario between those used. The statistical significance is expressed as follow: 0 ≤ \*\*\* < 0.001.


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.

**Figure 4.** Maximum suitability values grouped by altitudinal envelopes (100 m) across the whole country (**left**) and on the 7272 INFC2005 inventory plots only (**right**). In the last two pictures on the bottom, boxplots were colored according to the average value if below (red) or above (green) zero, expressing on average a decrease or increase of LS values, respectively, for the total forested area in the studied country.


**Table 2.** Maximum number of pixels that exceed of different threshold level and values of mean, standard deviation, minimum, maximum of altitude.


**Table 2.** *Cont.*

**Table 3.** Results of linear model function on the INFC 2005 domain. In this table, the variable fortype indicates the forest type category (i.e., beech forest, silver fir forest, etc.), the Gtot variable indicates the total basal area in m2 and finally, the variable TypeFor considered the typology of forest (if pure or mixed, this characteristic is established on the base of basal area of different species) and the forest tree species (if tree was coniferous or broadleaves). The statistical significance is expressed as follow: 0 ≤ \*\*\* < 0.001 ≤ \*\* < 0.01.


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.

**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.
