Species distribution modelling
For the modelling of the species distributions, we performed all models
in R version 4.0.3 (R Core Team 2020) with MaxEnt (version 3.4.1)
(Phillips et al. 2006) and dismo (version 1.3-3) (Hijmans et al. 2020)
and each model was replicated four times, using the cross-block
validation method in the ENMeval package version 2.0.3 (Muscarella et
al. 2014). Model evaluation measures were averaged across these 4
models. In total, we developed 68,288 models (see Appendix S1 in
supporting information) and an overview can be found in figure 1 and
table S5 following ODMap (Zurell et al. 2020). We ran all models in
parallel using the package snow version 0.4-3 (Tierney et al. 2018) and
the package parallel in base R version 4.0.3 (R Core Team 2020). The
performance was compared with the Area Under the Curve (AUC) of the
Receiver Operating characteristics curve, which is a threshold
independent evaluation method of the performance of the calibration and
evaluation dataset (Phillips et al. 2006; Elith et al. 2006). We also
looked at the percentage contribution of the variable to the model gain,
which is a common measure of variable importance (Bradie & Leung 2017).
The degree of overfitting was used to select the appropriate
regularization multiplier value of 5 (Appendix S2). The evaluation AUC
of the different resolutions were ranked per species to reduce
variability between species and find the optimal model settings for the
different groups of bees. We made prediction maps of the models with the
biotic variable at 1 km for the calculation of the corrected Akaike
information criterion (AICc; Burnham & Anderson 2002) and Continuous
Boyce Index (CBI; Hirzel et al. 2006). Presence and absence maps were
made based on the maximum training sensitivity and specificity threshold
(De Barros et al. 2012). The AICc and Continuous Boyce Index are
evaluation measures that include a penalty based on the number of
variables that is used in the models, unlike the AUC. The CBI, AICc and
evaluation AUC were compared between models including biotic factors and
excluding biotic factors with a one-sample Wilcoxon signed rank
(Wilcoxon 1945). We also compared the percentage contribution of the
different variable classes (land use, climate, soil, and biotic
variables) and individual variables, using a Kruskal-Wallis H test
(Kruskal & Wallis 1952) with a post-hoc Nemenyi test (Sachs
1997).