Introduction:
Over the last thirty years, there has been a massive increase in the number of published studies using species distribution models (SDMs) (reviewed in Lobo et al. 2010 and Melo-Merino et al. 2020). Species distribution models are used to identify areas of potentially suitable habitat by linking species occurrences to environmental variables (Loyola et al., 2012; Silva et al., 2014). These predictions of suitable habitats have many applications (Elith & Leathwick 2009), including: the estimation of potential distribution under different climate change scenarios (Marshall at al. 2018; Lima et al. 2020), the estimation of suitable areas for a species (Suzuki-Ohno et al. 2017) and assessing the potential invasiveness of an exotic species (Srivastava et al. 2019).
A key to understanding distribution patterns in insects lies in understanding their relationships with other organisms, e.g. pollinators and their floral resources (Willmer 2011). This biotic information is rarely included directly in distribution models, despite the fact that biotic interactions can theoretically improve our understanding and predictions of the distribution of species through different mechanisms (Wisz et al. 2013). Previous studies showed an improvement in the statistical performance of spatial models when including parasitic (Mathieu‐Bégné et al. 2021), facilitative (Heikkinen et al 2007), resource-consumer (Kissling et al. 2007; Bateman et al. 2012; Roslin et al. 2017; Atauchi et al. 2018; Herrera et al. 2018), competitive (Leach et al 2016; Mpakairi et al. 2017) and plant-pollinator interactions (Araújo & Luoto 2007; Espíndola & Pliscoff 2019; Kass et al. 2020). However, the improvements made by including biotic interactions depend on the spatial scale of the biotic variable in the model (Heikkinen et al 2007), and correlative relationships observed in models do not always represent biotic interactions (Giannini et al. 2013).
The biotic variable can be included in the SDM as an explanatory variable (Araújo & Luoto 2007; Kass et al. 2020) and it can also be introduced as either the raw distribution (Giannini et al. 2013; Leach et al. 2016) or a modelled distribution of the species that the modelled species interacts with (Bateman et al 2012; Giannini et al. 2013) (hereafter “interacting species”). The biotic variable can also be introduced as an approximation of the intensity of the interaction, e.g. the genomic background of parasite hosts can help identify populations with resistance genes (Mathieu‐Bégné et al. 2021), distance towards sighting of a competitor for competition (Mpakairi et al. 2017) or the distribution of diet resources for resource-consumer interactions (Araújo et al. 2014).
Besides the methodological considerations, spatial resolution may strongly affect the contribution of biotic interactions to modelled distribution patterns (Pearson and Dawson, 2003; Soberón & Peterson, 2005; Wisz et al. 2013). Heikkinen et al. (2007) showed that the impacts of facilitation between owls and woodpeckers are more visible at a resolution of 10 km than 40 km. This is consistent with Pearson and Dawson (2003), who hypothesized that at broader scales and coarse resolutions, climate variables are more dominant and biotic interactions less apparent (Heikkinen et al. 2007). However, an obligate parasite with a strong interaction with its host may always be more dependent on its interacting species’ distribution at any resolution. There is insufficient evidence as to how the explanatory power of biotic factors changes with spatial resolution, which is crucial for improving SDMs of species with strong hypothesized biotic interactions.
Whether the interaction is essential for survival depends on factors that include body size, dietary breadth, the distribution of the interacting species and the dependency on one another, e.g. whether the flowers are similarly dependent on the pollinator as the pollinator on the plant. The bee body size shows a strong relation with the foraging distance of different bees (Greenleaf et al. 2007; Kendall et al. 2019 and references therein) and smaller bees with a smaller foraging distance would require their host plant closer to their nest. It has been hypothesized that the dietary breadth of the species could influence the importance of the biotic factor in the models (e.g. specialist vs generalist; Araújo et al. 2014). In the case of bees, it has been shown that the population trend of specialist bees is correlated to the plant that they are dependent on for their pollen (Scheper et al. 2014). We expect that a smaller distribution range of the interacting species would have a higher contribution to the models, as it more likely to be the limiting factor of the modelled species. Specialist bee species have a tendency to decline more than generalist bee species and their decline is correlated to their host plant (Biesmeijer et al. 2006) and this leads us to expect that the specialist species show a higher contribution of the interacting species to their models.
In this paper, we aim to use a priori knowledge to investigate the importance of biotic interactions in species distribution models of bees at different spatial resolutions. Wild bees are a group of well-studied organisms that include species with a great importance to ecosystem resilience and that play a key role in pollination services to wild plants and crops (Kleijn et al. 2015; Senapathi et al. 2015; Weekers et al. 2022). Bees depend on pollen and nectar provided by plants and diets range from narrow (oligolectic bees, using few plant species) to broad (polylectic bees, using many plant species) (Rasmussen et al. 2020). Other species, up to 30%, are cleptoparasitic, meaning they are brood parasites which lay eggs in nests of other bee species (Cardinal et al. 2010). They may have one or multiple host bee species. The Netherlands is a good case study for the effects of biotic interactions on the distribution of wild bees, as there are more than 300 species of wild bees (Reemer 2018) and there is extensive data on plant-pollinator interactions, hosts of cleptoparasitic bees and occurrence data. By integrating knowledge of plant visitation and cleptoparasitic interactions, we aim to (1) assess the importance of biotic factors in explaining distributions of polylectic, oligolectic and cleptoparasitic bees and (2) assess the relative importance of different factors, including flower and host specialization, spatial resolution, taxonomic level, distribution of the interacting species and bee body size (which relates to foraging range) in explaining the contribution of the interacting species to the models. By modelling a large number of bee species and using different input variables and methods, we identified important factors that are related to the implementation of the biotic interactions to the models.