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

$$LS = intercept + \beta\_1 \cdot CM + \beta\_2 \cdot TH + \varepsilon \tag{1}$$

where *CM* represents the different climate simulation model we used, β<sup>1</sup> and β<sup>2</sup> 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.
