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