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