MetaComNet: A random forest‐based framework for making spatial predictions of plant–pollinator interactions

Predicting plant–pollinator interaction networks over space and time will improve our understanding of how environmental change is likely to impact the functioning of ecosystems. Here we propose a framework for producing spatially explicit predictions of the occurrence and number of pairwise plant–pollinator interactions and of the species richness, diversity and abundance of pollinators visiting flowers. We call the framework ‘MetaComNet’ because it aims to link metacommunity dynamics to the assembly of ecological networks. To illustrate the MetaComNet functionality, we used a dataset on bee–flower networks sampled at 16 sites in southeast Norway along with random forest models to predict bee–flower interactions. We included variables associated with climatic conditions (elevation) and habitat availability within a 250 m radius of each site. Regional commonness, site‐specific distance to conspecifics, social guild and floral preference were included as bee traits. Each plant species was assigned a score reflecting its site‐specific abundance, and four scores reflecting the bee species that the plant family is known to attract. We used leave‐one‐out cross‐validations to assess the models' ability to predict pairwise plant–bee interactions across the landscape. The relationship between observed occurrence or absence of interactions and the predicted probability of interactions was nearly proportional (GLMlogistic regression slope = 1.09), matching the data well (AUC = 0.88), and explained 30% of the variation. Predicted probability of interactions was also correlated with the number of observed pairwise interactions (r = 0.32). The sum of predicted probabilities of bee–flower interactions were positively correlated with observed species richness (r = 0.50), diversity (r = 0.48) and abundance (r = 0.42) of wild bees interacting with plant species within sites. Our findings show that the MetaComNet framework can be a useful approach for making spatially explicit predictions and mapping plant–pollinator interactions. Such predictions have the potential to identify areas where the pollination potential for wild plants is particularly high, and where conservation action should be directed to preserve this ecosystem function.


| INTRODUC TI ON
Nearly nine of every 10 species of flowering plants rely on interactions, mainly with insects, for cross-pollination (Ollerton et al., 2011). However, wild plants are experiencing a shortage of pollinators in both natural and managed terrestrial ecosystems (Bennett et al., 2020), which can considerably affect plant reproduction and persistence (e.g. Gomez et al., 2010;Ollerton, 2017;Thomann et al., 2013). Our ability to predict how plant communities will respond to environmental change partly relies on our ability to predict how plant-pollinator interactions vary in heterogeneous landscapes (Tylianakis & Morris, 2017).
Most approaches to modelling plant-pollinator interactions in heterogeneous landscapes have focused on network properties, such as modularity (reviewed in Pellissier et al., 2018), but these properties may not fully encapsulate the underlying community assembly processes (Olito & Fox, 2015). An alternative approach is to model pairwise plant-pollinator interactions directly, and to derive network properties by aggregating model predictions (Graham & Weinstein, 2018), for instance by including plant and pollinator abundances, traits and phylogenies (e.g. Benadi et al., in press;Pichler et al., 2020;Stock et al., 2021). However, Benadi et al. (in press) found a drop in model performance when predicting into novel habitats, which may arise from differences in habitat environmental conditions influencing pollinator distributions (e.g. Hoiss et al., 2012). Hence, accounting for processes behind community assembly (as per Vellend, 2016) is central for making spatial predictions of plant-pollinator interactions. We therfore hypothesise that the number of pairwise interactions between plants and pollinators can be effectively modelled as a function of plant and pollinator affiliations and of variables underlying pollinator community assembly at different spatial scales. We propose a conceptual framework to make spatial predictions of pairwise interactions between plants and pollinators, which we refer to as 'MetaComNet' because it aims at linking meta-community structuring factors to the structure of ecological networks.
MetaComNet is pollinator-oriented in that it focuses on predicting the occurrence, or number, of plant-pollinator interactions by modelling the distribution of wild bees across plant species while considering regional and landscape level factors. The response variables considered are the number, or presence/absence, of observed interactions between pollinator and plant species observed in specific localities. We focus on pollinators because they tend to display more pronounced spatial turnover in interaction networks than plants do (Trøjelsgaard et al., 2015). MetaComNet builds on Tylianakis and Morris (2017) who define the occurrence, or number, of interactions between plants and pollinators as the endproduct of processes and conditions that determine species composition at different spatial scales (Figure 1). We pose the following hypotheses relevant at different spatial scales: Regional scale: If community, or network, assembly is ecologically neutral, the abundance of species within communities (Vellend, 2016) will be proportional to species commonness in the regional species pool, in turn dependent on climatic requirements and the biogeography of species (Cornell & Harrison, 2014). Accounting for differences in regional commonness is also important because species abundance is related to the level of random interactions (reviewed in Krishna et al., 2008;Tylianakis & Morris, 2017;Vázquez et al., 2007).
In MetaComNet, the regional commonness of pollinators is included as a regional level predictor of pairwise interactions to account for neutral network assembly processes.
Landscape scale: Pollinator communities are assembled through dispersal processes (Hagen et al., 2012) and through mechanisms of species sorting determined by the suitability of an area as habitat for the species (environmental filtering). Dispersal rates depend in part on geographic distance and barriers (Carstensen et al., 2014;Trøjelsgaard et al., 2015); hence the likelihood of a species occurring in a habitat patch, and of it interacting with the plants therein, is expected to decrease with the distance to the nearest population in the region. In MetaComNet, we accounted for the influence of immigration rates by including site-and species-specific estimates of the distance to the nearest population. To support viable pollinator populations, a location must contain enough nesting and foraging resources within enough proximity (Westrich, 1996). The amount of semi-natural habitat or degree of landscape diversity within a 250 m radius can be used as a proxies for habitat amount since solitary bee of predicted probabilities of bee-flower interactions were positively correlated with observed species richness (r = 0.50), diversity (r = 0.48) and abundance (r = 0.42) of wild bees interacting with plant species within sites. 4. Our findings show that the MetaComNet framework can be a useful approach for making spatially explicit predictions and mapping plant-pollinator interactions. Such predictions have the potential to identify areas where the pollination potential for wild plants is particularly high, and where conservation action should be directed to preserve this ecosystem function.

K E Y W O R D S
interactions, network, plants, pollinators, predict, random forest species richness increases with habitat area at this scale (Steffan-Dewenter et al., 2002). The distance to soil deposits with high sand concentrations can be used as a proxy for distance to high quality nesting substrates, because many species prefer such substrates for nest sites (Antoine & Forrest, 2021;Heneberg et al., 2013).
Local scale: At the scale of individual flowers within habitats, the occurrence and number of bee-plant interactions will depend on the attractiveness of the flower. Flower attractiveness depends in part on their relative abundance (Fowler et al., 2016;Stavert et al., 2019), even though visitation rates to flowers may saturate (Totland, 1994) or show unimodal responses (Benadi & Pauw, 2018). Other factors determining flower attractiveness to pollinators are, for example, particular morphological characteristics (e.g. Benadi et al., in press;Pichler et al., 2020;Stock et al., 2021), non-visual cues such as floral scent (Larue et al., 2016) and pollen toxicity that require adaptations to overcome (reviewed in Rivest & Forest, 2020). Trait-matching reduces the frequency or even excludes some combinations of partners in plant-pollinator networks (Olesen et al., 2011). However, which floral traits select for specific bees is not always easy to predict and species may respond to different traits on the same plant (Rowe et al., 2020). Also, in bumble bees, pollen preferences can be more directly related to phylogenetic relationships than to proboscis length (Wood et al., 2021). In MetaComNet, trait-matching between plants and pollinators can be accounted for by assigning a set of functional and/or floral preference traits reflecting the expected plant-pollinator associations.
The MetaComNet model integrates data from the three geographic levels indicated above into a data frame illustrated in Table 1. The model parameters include response variables (number, or the presence or absence of interactions), grouping variables (pollinator species, plant species, site identity) and predictor variables such as pollinator traits, plant traits and site-specific environmental F I G U R E 1 Conceptualisation of the MetaComNet framework and the hierarchical assembly of plant-bee interaction networks. The likelihood of a bee interacting with a plant species growing in a habitat patch is a product of a series of scale-dependent conditions, that is, regional commonness of pollinators, landscape habitat suitability (nesting and foraging resources) and distance to neighbouring populations, and local composition of flower plants, their abundance and level of attractiveness conditions. Using empirical data from Norway, we illustrate the ability of MetaComNet to predict the occurrence and number of interactions between wild bees and plants.
MetaComNet uses a random forest modelling framework (Breiman, 2001) to model pollinator-plant associations because of its ability to define nonlinear interactions between predictor variables. Random forest models also offer an intuitive way to view the hierarchical assembly of pairwise interaction networks ( Figure 1). Moreover, random forest models have been shown to outperform other machine learning and glm-based techniques when predicting pairwise plant-pollinator interactions in both simulated and empirical data (Pichler et al., 2020) and to perform as well as mechanistic likelihood-based models (Benadi et al., in press). We used three random-forest approaches to test if predictor variables (see Table 2) could predict spatial variation in: • Occurrence and number of pairwise interactions between wild bee species and plants.
• Species richness, Shannon diversity and abundance of wild bee floral visitors.

| Bee-flower network sampling
We sampled bee-flower networks along 16 roadsides (sites) in Southeast Norway in 2017. Eight of the study sites were located on sandy sediments and the remaining eight were located on clay dominated sediments (Skoog, 2018). At each study site, flowervisiting bees were collected during 1 hr by two observers along a 50 m transect, once during early to mid-July and once during early The site-specific plant abundance was used as a proxy for the number of flowers from that species within a site, and was based on species which were flowering at the time of the inventory. Plants that did not occur within any of the six 1 m 2 vegetation quadrats within a site, but with bee visitation records were given a plant cover abundance value of 0.05 to indicate that the plant was locally present.
We assigned variables to the recorded plants and bees that we expected would influence the number of plant-pollinator interactions. For plants we used the site abundance since we expected a high correlation of bee visitations with plant abundance. Because our aim was not to identify functional traits that determine beeflower interactions, but to account for floral associations of bees in our models, we assembled a binary network of 207 bee species and 61 plant family interactions based on existing information from interaction records sampled at a greater temporal and spatial extent than our study area (Rasmussen et al., 2021;Wood et al., 2021). This approach provided a more inclusive measure of host plants of bees, closer to the fundamental niche, than what would be achieved from a single survey, such as ours. We used a detrended correspondence analysis (DCA) in the Vegan package in r (Oksanen et al., 2018) to establish the four main axes of correspondence between plant families and non-parasitic bee species associations ( Figure S1). We used the plant and bee DCA scores to account for plant-bee associations. TA B L E 1 Variables included in the MetaComNet network model. The data frame contains columns with response variables including: (i) number, or presence or absence, of observed interactions between a pollinator species and a plant species in a particular study site. Grouping variables include: pollinator species (Pol.), plant species (Plnt.) and site identity (Site). Predictor variables that can be linked to pollinator species, plant species or the site. Pollinator variables include the regional commonness (RC); the distance to the nearest known population of that species (Dst); and traits (Tr8s) related to environmental requirements or the floral preference of the species. Plant variables include the local commonness (LC), or abundance, and traits (TR8s) related to pollinator affiliations. Site-specific variables include geographic coordinates that are used for calculating distances to potential source populations and for extracting georeferenced environmental variables (Env) such as m a.s.l., or area of semi-natural habitat within 250 m radius in the surrounding landscape In our model we distinguished between solitary and social wild bees (Bombus spp.) because solitary bee diversity responds more strongly to landscape conditions at local spatial scales than bumble bees (Steffan-Dewenter et al., 2002). We used records from GBIF.org (GBIF, 2021) to estimate the regional commonness of bee species and the distance to the nearest known occurrence of each species. We downloaded occurrence records covering our region ( Figure 1, x min = 9.99, x max = 12.19, y min = 59.09, y max = 61.31, projection: WGS84) and excluded records older than 20 years and with a coordinate uncertainty >100 m. To estimate regional commonness, we tallied the number of 10 km grid cells within which a species had been observed. For each site we calculated the geographic distance to the nearest GBIF record of each species.
For each site we obtained information on environmental conditions known to affect wild bee distributions. As a proxy for climatic conditions, we used elevation above sea level, obtained from a digital elevation model with a 50 m resolution (Norwegian Mapping Authority, 2016). As proxies for landscape conditions, we used the European ELC10 land cover map (Venter & Sydenham, 2021) to calculate: the proportion of grassland area and Shannon landscape diversity within a 250 m radius. As a proxy for nesting conditions, we calculated the distance to sand-dominated geological deposits, that is, soils with a high permeability (Geological Survey of Norway, 2011).

| Modelling and predicting empirical bee-plant interactions
We assembled a data frame where each row was defined by a study site, a plant species found within the study site, and one of the bee species occurring in the 16 study sites (Tables 1 and 2 We used random forest (Breiman, 2001) with the ranger package (Wright & Ziegler, 2017) via the Caret package (Kuhn, 2018) in r (R Core Team, 2020). We fitted three models to the data depending on the type of response variable: the presence versus absence of interactions, using classification trees; the presence (one) versus absence (zero) of interactions using regression trees; and number of interactions using regression trees. The three resulting model outputs were: predicted probability of interactions, that is, class probability; predicted frequency of interaction, that is, predicted proportion of presences; and predicted number of interactions. We used leaveone-out cross-validation, by iteratively training models on data from 15 sites and predicting onto the remaining site to allow assessing model performances across all 16 sites as well as their variability in terms of predicting pairwise interactions within sites. We also con- For all three models we used NagelKerkes log likelihood-based R 2 from r package MuMIn (Barton, 2018), from the GLM models, and the area under the curve (AUC) from r package prOC (Robin et al., 2011) to assess model performance. We calculated the regression slopes for the first two models (with presence/absence as response variables), and the average R 2 and AUC. To assess the models' power to predict the number of interactions, we calculated the Pearson correlations between the observed number of interactions and predicted probability, frequency and number of interactions. All validation metrics were calculated by (a) including predictions and observations from all 16 sites, (b) and by calculating the mean, and standard deviation from predictions for each site individually. We calculated the scaled importance of predictor variables from each of the 16 models in order to assess if the variation in variable selection across models differed between the three modelling strategies.

| Predicting flower-visitor richness, diversity and abundance
We tested the level of correspondence between the sum of predicted pairwise bee-plant interactions and observed flower-visitor species richness, diversity and abundance within sites. For each of the three models, we calculated predicted flower-visitor species richness from the sum of predicted probabilities of interactions and from the sum of predicted frequencies of interactions for each plant species per site. We also calculated the predicted abundance of flower visitors as the sum of predicted number of interactions across bee species.
We then calculated the Pearson correlation between predicted species richness or abundance and observed: flower-visitor species richness; species diversity; and species abundance.

| Mapping flower-visitor species richness
To illustrate how the predicted flower-visitor species richness can be mapped and thus used to identify areas where plants are most or least likely to be pollen limited, we re-fitted the random forest regression to the occurrence of interactions using all the data from the 16 sites and used this model to produce prediction maps. We created predic-

| Predicting pairwise interactions between bee and plant species
There was a positive relationship between observed and predicted occurrence and number of bee-plant interactions irrespective of the modelling strategy (Figure 2a Bee and plant species-based leave-one-out validations showed that classification trees and regression trees on occurrence or absence of interactions, performed equally well and outperformed regression trees based on number of interactions when attempting to predict pairwise interactions for bees and plant species (Figures S2 and S3).
Comparing predictor variable importance (Figure 2c,f,i) across the 16 models for each of the three modelling strategies also showed that models for interaction frequencies (Figure 2f)

| Predicting flower-visitor species richness, diversity and abundance
Predicted bee species richness and abundance were positively correlated with observed flower-visitor species richness, diversity and abundance. The Pearson correlation coefficient between observed flower-visitor species richness, diversity or abundance, and predicted species richness was similar for classification (Figure 3a-c) and regression random forest models (Figure 3d-e). In both cases median values for observed species richness, diversity and abundance, increased with predicted species richness but the rate of increase seemed to saturate at a predicted species richness at and above two. In comparison, correlations between predicted flowervisitor abundance and observed flower-visitor species richness, diversity and abundance were weaker (Figure 3g-i).

| Mapping flower-visitor species richness
The random forest regression models produced the most accurate For Compositae, areas with predicted values at or above the third quantile were mainly concentrated on areas with sand dominated geological deposits (around the Gardermoen airport, Figure 4d; Figure S2). For Leguminosae, the area predicted with highest diversity of wild bees was found southwest of the Airport (Figure 4g), where the soil substrate is dominated by marine, clayish deposits.

| DISCUSS ION
The aim of this study was to develop and test a framework for producing spatially explicit predictions of plant-pollinator networks.
Despite the relatively low predictive importance of landscape level variables, there was a considerable spatial difference in the predicted species richness of bees that visit plants belonging to Compositae and Leguminosae (Figure 4).
Approaches to modelling plant-pollinator networks can be classified according to the main strategies of species distribution modelling identified by D' Amen et al. (2017). In the 'assemble first and predict later' approach (sensu D' Amen et al., 2017), network indices are modelled as functions of environmental conditions (reviewed in Pellissier et al., 2018). The advantages of this strategy are that structural properties of entire networks are captured; and network indices have hypotheses affiliated to their drivers and relationships to ecosystem functioning. A problem with this strategy is that different species compositions can result in similar network properties (Olito & Fox, 2015), so that processes such as, for example, competition, trait-matching and neutrality can all theoretically give rise to similar degrees of modularity within ecological networks (reviewed in Dormann et al., 2017). An alternative strategy is the 'predict first assemble later' (sensu D' Amen et al., 2017), where distributions of flower visitors are modelled individually, and the resulting predictions are then aggregated to network properties. This approach is typically adopted in single plant species systems and, for example, used to predict pollinator abundances in crops (e.g. Gardner et al., 2020;Lonsdorf et al., 2009). A drawback of this approach is that because species are modelled individually, it is difficult to link predictions to community assembly processes (sensu Vellend, 2016) and meta-community ecological theory (sensu Leibold et al., 2004).
In the final strategy, the 'assemble and predict together' strategy (sensu D'Amen et al., 2017), species interactions are modelled simultaneously for all pollinator species across plant species and communities. This strategy has been adopted in recent frameworks for predicting plant-pollinator interactions and allows linking network assembly to metacommunity ecological theory such as trait-based species sorting among communities (Leibold et al., 2004). For instance, Graham and Weinstein (2018) predicted plant hummingbird interactions simultaneously for hummingbird species along a gradient of elevation and outlined a strategy for how functional traits could be integrated into their modelling framework. In addition to trait-based pollinator species sorting, the role of pollinator abundances has previously been included into models of plant-pollinator interactions that identify linkage rules within networks (Bartomeus et al., 2016). Furthermore, in addition to trait-based filtering, and random (i.e. abundance-based) encounters, plant-pollinator interactions, for example, between mustards and wild bees, have been shown to depend on habitat isolation (Steffan-Dewenter & Tscharntke, 1999). By including the influence of dispersal limitation, in addition to those of trait-based species sorting and random encounters, MetaComNet extends existing 'assemble and predict together' frameworks for predicting plant-pollinator interactions.
A direct comparison of the predictive power of MetaComNet and other frameworks for predicting pairwise interactions (e.g. Pichler et al., 2020;Stock et al., 2021;Benadi et al., in press) is partly hindered by differences in study designs, and particularly that we in our models included spatial predictors in order to predict across networks of pairwise interactions (sites). For instance, Stock et al. (2021) devised a cross-validation strategy to assess withinnetwork model prediction of: pairwise interactions; interactions for particular bee or plant species; or particular combinations of bees and plants. The most similar strategy to that we followed would be that of the 'pairwise interactions', however, in some instances particular plant species occurring in the validation data, would not have been found in the combination of sites used to train our models.
Predicting across networks will therefore produce a mix of the crossvalidation strategies proposed by Stock et al. (2021) and we only found marginal differences in model performances when predicting pairwise interactions within and across sites ( Figure 2) and when predicting interactions for individual species of bees ( Figure S2) and plants ( Figure S3). However, despite including spatial variables and cross-validating predictions into new sites which we would expect would reduce model performances, our models yielded AUC statistics for predictions of interaction occurrences that were comparable to those obtained using simulated and empirical data for within network predictions (Pichler et al., 2020;Stock et al., 2021). This may suggest that plant-pollinator trait-matching and neutral, abundancebased, processes, are likely to be the main structuring processes F I G U R E 3 Relationships between observed flower-visitor species richness, diversity or abundance and predictions flower-visitor species richness and abundance modelled with three random tree approaches: (i) predicted probability of interactions based on classification trees (a-c); (ii) predicted frequency of interactions from regression trees on presence/absence data (b-h); (iii) and predicted number of interactions from regression trees (g-i). The Pearson correlation coefficient between observed and predicted values is shown for each of the modelling approaches. Data points are plotted as red points together with their boxplot summary statistics behind network assembly, although our models did identify spatial signals in plant-pollinator interactions (Figure 4).
We used detrended correspondance scores to reflect the floral associations of bees, and bee associations of plant families ( Figure   S2). Another option would be to use morphological and phenological traits instead (Pichler et al., 2020). A benefit of the latter approach is that novel interactions can be predicted, for example, for invasive species, if one has information on their traits. Alternatively, or additionally, to traits one can use information on phylogenetic relatedness when modelling interactions (e.g. Benadi et al., in press;Stock et al., 2021) because plant-pollinator associations are often somewhat phylogenetically conserved (e.g. Wood et al., 2021). If information on host plants of some pollinator species is not known, and one wishes to apply the approach using DCA scores adopted in F I G U R E 4 Spatial predictions of flower-visitor species richness in a representative landscape in the study area, the area surrounding Oslo airport Gardemoen (a). Predicted flower-visitor species richness to plant species in the Compositae (b-d) and Leguminosae (e-g) families depended on plant abundance, shown by predicting with plant abundances held constant at low (b, e) and high (c, f) levels. The predicted flower-visitor species richness differs spatially between the two plant taxa (d, g), illustrated by masking out areas with predicted values less than the 75th quantile of predicted values from (c) and (f) respectively. Satellite imagery from Copernicus Sentinel-2 data (2019)/processed by the Norwegian Mapping Authority this study, a potential would be to use the average DCA score values from the closest relatives for which one has information, or assume that the species will be restricted to the same host-plant families as in its native range (Vaudo et al., 2020). Using floral preference traits, inferred from, for example, DCA scores, should therefore not nescecarily hinder the prediction of novel interactions. However, a limitation of using floral association (e.g. DCA) scores is that they do not provide tests of how specific plant-pollinator trait-combinations influence pairwise interaction probabilities. Thus if an aim is to identify linkage rules in ecological networks then, hypothesised, functional traits should be used in the MetaComNet framework instead of floral preference scores.
In our models, we treated absences of interactions with the same degree of confidence as presences. However, because one is unlikely to detect all interactions when sampling plant-pollinator interaction networks (Chacoff et al., 2012), data on multiple interaction net-  (Liu et al., 2015). However, this approach requires some degree of subjectivity in terms of which traits are included when estimating and setting thresholds for interaction credibility. We are also unsure how one would weight the credibility of an interaction in a spatial setting, where some pairwise interactions might be considered credible, but were unobserved because of environmental conditions at larger spatial scales than within the network. Removing proportional to the predicted probabilities of occurrences (Elkan & Noto, 2008). We therefore suggest that absence values, and the information they contain, are included in models aimed at predicting plant-pollinator interactions.
Despite a reasonable fit to the validation data, a considerable amount of variation in pairwise bee-plant interactions was left unaccounted for by our models (Figure 2a,d,g). While some of this unexplained variation is likely attributable to random error, it seems fair to assume that a large fraction of it was due to unmeasured predictor variables. Negative biotic interactions, such as interspecific competition which can suppress bee-flower visitations (Wignall et al., 2020), were not accounted for in our models.
However, we would expect that negative biotic interactions would result in our models overestimating local interaction probabilities, which was not the case. By contrast, our models tended to underestimate the occurrence of pairwise interactions (i.e. Figure 2a,d, regression slopes >1) which may suggest that there were elements related to habitat conditions, such as habitat continuity (Morandin & Kremen, 2013), that could result in higher-than-expected occurrences of bees. A potential solution for future implementations of MetaComNet would be to use annually updatable land cover maps (e.g. Venter & Sydenham, 2021) to estimate the continuity of habitat patches. It could also be that bee species were found more frequently than predicted because the distance to source population variable was too conservative, that is, because species occurrences in the GBIF record did not reflect the actual distribution of species. Ideally one would have more information on the location of potential source populations, than from potentially scattered, and potentially biased, species observation records. Another important contributor to the unexplained variation in our models was the small number of study sites available (15 for building models, 1 for validating). In order to accurately estimate the contributions of landscape level variables, we therefore suspect that increasing the sample size of study sites will enable stronger predictions of bee-flower interactions than what was possible given the limited data available to us. Despite these potential limitations, our models did produce prediction maps that correspond well with how we would expect floral visitation patterns to be distributed (Figure 4), that is, with Compositae-visiting solitary bees being concentrated near sandy soils, and Leguminosae-visiting bumble bees being less concentrated on sandy soils ( Figure S2).

| CON CLUS IONS
Spatial models of pollinator diversity should theoretically allow for identifying areas (a) where plant populations would benefit from pollinator enhancement schemes, and (b) where plant populations are likely to have a higher genetic diversity because of high levels of pollination. We believe that the modelling framework presented here provides a promising avenue for producing spatially explicit predictions of plant-pollinator interactions. If possible, future studies adopting the framework should attempt to assess the degree of pollen limitation within plant populations to empirically test if predicting interaction partner diversity allows to also predict the degree of pollen limitation experienced by plant species.

ACK N OWLED G EM ENTS
We thank the Norwegian public roads administration for facilitating our field work along road sides. The associate editor, and Dr Gita Benadi and one anonymous reviewer provided thoughtful and constructive reviews of previous versions of the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/2041-210X.13762.

DATA AVA I L A B I L I T Y S TAT E M E N T
R code for reproducing analyses can be found at https://doi. org/10.5281/zenodo.5644742 (Sydenham et al., 2021a). Data for running the R code are available through the Dryad Digital