Species distribution models
To work out species distribution models (SDMs), we relied on the method proposed by Brambilla et al. (2022). Such an approach involves building maximum entropy models using MaxEnt (Phillips et al. 2006) in the software R (R Core Team 2020) by combining different packages. We used only MaxEnt because of the many advantages it offers advantages over alternative methods: it is the commonest algorithm for SDMs, it limits the potential undesired effects of false absences (Jiménez-Valverde et al. 2008, Elith et al. 2011), often performs better than other methods or ensemble modelling (Kaky et al. 2020), and has been already used to model distribution and potential changes for other species on similar same mountain ranges (Brambilla et al. 2022). We scattered 79393 background points (the highest possible number of independent locations) within a 10-km buffer drawn around all bumblebee records, to ensure that background points are placed in areas actually sampled or close to sampled ones, to adequately represent sampled environmental conditions (Brambilla et al. 2020).
By means of the ‘checkerboard 2’ method of the ENMeval package (Muscarella et al. 2014), occurrence data of each species were partitioned into two spatially independent datasets. In case of records of the same species overlapping within the same grid cell, they were considered as duplicates and only one was retained for the analyses. Training datasets included occurrences from three partitions, and testing datasets those from the remaining fourth one. We only fitted linear and quadratic relationships to reduce possible overfitting. The regularization multiplier was first selected by testing 0.5-increase values between 0.5 and 5, and that leading to the lowest AICc was chosen to build a base model.
Then, all the variables showing lambda equal to 0 (i.e., no tangible effect on species distribution) were discarded. A variable selection procedure was then carried out by leaving out one variable at a time according to increasing value of permutation importance, until the model’s AICc increased. We thus identified a most supported model, which was then subject to further tuning. Linear and quadratic features and the value of the regularization multiplier were checked again, if needed changed (always according to AICc), and a final model was thus produced and used for model evaluation and distribution prediction. The Area Under the Receiver Operating Characteristic (ROC) Curve (AUC), as well as the True Skills Statistics (TSS), were computed over training and testing data sets and used for model evaluation, together with the computation of the omission rates over the test dataset, both at the 10th percentile and at the minimum training presence, both computed on the training presence dataset. The two omission rates should be close to 0.1 and to 0, respectively, whereas AUC and TSS should be similar over the training and testing datasets as their absolute value is poorly informative (Lobo et al. 2008). The final model was used to predict a species’ environmental suitability according to its cloglog outputs. For further details, see Brambilla et al. (2022).
The distribution models obtained according to the above procedure were then used to predict environmental suitability over current and future conditions. From each map of predicted suitability, we derived a potential range by considering as suitable for a species all cells with an environmental suitability value higher than the tenth percentile threshold, considering the cloglog-transformed output. We selected such a threshold over the possible other ones as its use led to the results most consistent with the known actual distribution of the target species (Brambilla et al., 2022). To estimate the extent of the relative range under current and future conditions, the potential distribution so obtained had been overlapped with the extent of the mountain ranges respectively occupied by each studied species. Since mountain ranges were not available in public repositories as shapefiles, each mountain range was identified by selecting areas above 300 m a.s.l. (elevation threshold was taken from Kapos et al., 2000) within the commonly recognized geographic boundaries of Alps, Apennines and Pyrenees.
We calculated in-situ and ex-situ refugia from the distribution inferred from the above models and projections. We defined “refugia” the areas that probably will be suitable “in-situ” or “ex-situ” for each target species and region, following Brambilla et al. (2017): (1) “in-situ” areas suitable under current and all future conditions, (2) “ex-situ” potential areas, currently unsuitable for a species, but suitable under all future models. While “in-situ” refugia are fundamental for population resistance, “ex-situ” ones are key sites for future redistribution and resilience. Moreover, the refugia were overlapped with the current Protected Area network, obtained by merging Natura 2000 sites with the European inventory of nationally designated PAs (Nationally designated areas; CDDA), updated in 2020 (https://www.eea.europa.eu/data-and-maps/data/nationally-designated-areas-national-cdda-15; accessed 2 Feb 2021).