Figure captions:
Figure 1: Schematic overview of the modelling workflow. The elements are represented in boxes and consist of models, information from databases or variables (brown), variables, indices, evaluation measures or experimental conditions (e.g. no data is used or the resolution of the used variables). The arrows indicate the information flow. The evaluation measures are: (i) the corrected Akaike information criteria (AICc; Burnham & Anderson 2002), (ii) the continuous Boyce index (Boyce index; Hirzel et al. 2006) and (iii) area under the curve of the receiver operating characteristics curve (AUC) applied to species distribution models (Elith et al. 2006).
Figure 2: The effect of the resolution and taxonomic level on the contribution of the biotic variable to the model, expressed as the ranking of the biotic variable contribution per species (from high to low: 1-5; 2A) and the difference in variable contribution between the species that the modelled species interacts with (interacting species) added at species and genus taxonomic level per species (2B). The arrows indicate the direction, where the variable contribution is the highest for the respective taxonomic level. The resolution is the scale in longitudinal and latitudinal direction at which the interacting species is observed. The gray area is the standard deviation.
Figure 3: The differences between models including host plant or parasitic host interactions and models with only land use, climate and soil variables. Evaluation measurements include Area Under the receiver operating characteristic Curve (AUC) value of the evaluation dataset (figure 3A), Continuous Boyce Index (CBI) of the evaluation dataset (figure 3B), Aikake Information Criteria for small sample sizes for both the evaluation and calibration data (AICc; Burnham & Anderson 2002; figure 3C) and AICc of the evaluation data only (figure 3D). Host plants and hosts of parasites were either included at the species or genus level. The difference in evaluation metrics for models with and without biotic factors, or difference from zero, is tested for significance with a One-Sample Wilcoxon Signed Rank Test (p<0.05*; p<0.01**, p<0.001***). For the AICc both the calibration and evaluation dataset were included, because 66 modelled bee species did not have enough evaluation datapoints to calculate the AICc.
Figure 4: The different boxplots represent the summed contribution of the five climate variables, the sixteen land use variables, the eight soil variables and the single biotic variable, averaged over the modelled species in the functional groups. The biotic variable is averaged over the species and genus taxonomic level of the visited plant or host bee. The different letters above the boxplots indicate significant differences between variable groups within the functional trait groups (p<0.05).
Figure 5: The comparison of the biotic interaction models to models with random interactions, described as null models, with plants (for the oligolectic and polylectic bees) or bees (for the cleptoparasitic bees). Figure 4A shows the distribution of the performance of the biotic interaction models, expressed as the rank of the evaluation AUC among all interaction models divided by the total number of models. The y-axis represents the total number of modelled species that fall within the performance threshold on the x-axis. For example, the performance in evaluation AUC of the known interaction was compared to the other 306 plant species and ranked based on the position. If the known interaction was the third best performing model, the modelled species would have the value of 0.98% (the percentage rank would be 3/307 * 100 = 0.98%) and fall within 0-2.5% best performing models. The two lines indicate the threshold of 5% and 25% best performing models. Figure 4B summarizes the results, comparing the percentage of modelled species that fall within the 5% best performing ranks, indicating a significant difference from the null models with p < 0.05 (5% best performing models), and 25% best performing ranks. Although the percentage of models that fall within the 5% best performing models is higher for the oligolectic bees and cleptoparasitic bees, the polylectic bees show a high percentage of performance within the 25% best performing models, showing a more general preference of biotic interactions. The number of random interactions for every null model are 306 interactions for the flower visiting bees with the interacting species at species level, 160 interactions for the flower visiting bees with the interacting species at genus level, 99 interactions for the cleptoparasitic bees with the interacting species at species level and 100 or 15 interactions for the cleptoparasitic bees with the interacting species at genus level (see Appendix S3 in supporting information).
Figure 6: The results of the Generalized Linear Models (GLMs) show the relation between flower specialization (Shannon-Wiener index of number of plants genera interacted with) and the contribution of the biotic variable to the models of the oligolectic and polylectic bees (Shannon 1948; figure 6A). Figure 6B shows effect of distribution of the most visited genus on the contribution of the biotic variable to the model. Figure 6C shows the relation between the distribution of the host bees and the contribution of the biotic variable to the models of the cleptoparasitic bees.