Honey bees (Apis mellifera) modify plant-pollinator network structure, but do not alter wild species’ interactions

Abstract

Honey bees (Apis mellifera) are widely used for honey production and crop pollination, raising concern for wild pollinators, as honey bees may compete with wild pollinators for floral resources. The first sign of competition, before changes appear in wild pollinator abundance or diversity, may be changes to wild pollinator interactions with plants. Such changes for a community can be measured by looking at changes to metrics of resource use overlap in plant-pollinator interaction networks. Studies of honey bee effects on plant-pollinator networks have usually not distinguished whether honey bees alter wild pollinator interactions, or if they merely alter total network structure by adding their own interactions. To test this question, we experimentally introduced honey bees to a Canadian grassland and measured plant-pollinator interactions at varying distances from the introduced hives. We found that honey bees increased the network metrics of pollinator and plant functional complementarity and decreased interaction evenness. However, in networks constructed from just wild pollinator interactions, honey bee abundance did not affect any of the metrics calculated. Thus, all network structural changes to the full network (including honey bee interactions) were due only to honey bee-plant interactions, and not to honey bees causing changes in wild pollinator-plant interactions. Given widespread and increasing use of honey bees, it is important to establish whether they affect wild pollinator communities. Our results suggest that honey bees did not alter wild pollinator foraging patterns in this system, even in a year that was drier than the 20-year average.

Introduction

With the widespread use of managed honey bees (Apis mellifera, Linnaeus, 1758) for crop pollination and honey production in parts of the world outside their native range in Europe, Africa, and the Middle East [1], many studies have suggested that high honey bee densities negatively affect wild pollinators through interference or exploitative competition (reviewed in [2, 3]. Wild insect pollinators include many species, including other bees, wasps, and many other orders, including flies, butterflies, moths and beetles, which are diverse in size, sociality, and floral trait preferences [4]. Honey bees visit a wide variety of flower species [5], and high honey bee densities have been linked to declines in wild pollinator diversity [6], abundance [7], floral resource use [8], and fitness [9], all via competition [10].

However, the first sign of competition between honey bees and wild pollinators may be changes to which flower species wild pollinators visit in the presence of honey bees, versus when honey bees are not present [7, 8]. If a wild pollinator is displaced by honey bees from one or more of its usual flower species (as reported for other competitive relationships between pollinators; e.g., [11, 12], and for honey bees, [9], the total set of flowers that the wild pollinator visits will either decrease (as in [13]), or will shift to include new flower species not previously visited by that pollinator species [14]. Additionally, if honey bee use of a flower depletes nectar or pollen [15], the frequency of visits by wild pollinator species to that flower species may decrease, even if visitation does not cease altogether [16]. The ecological effects of these changes may cascade across the plant-pollinator community, in that a wild pollinator species displaced by honey bees onto a new flower species may in turn displace the pollinators that had previously focused on that flower species [17]. Furthermore, when honey bees dominate floral resources [9], this may also affect a plant’s ability to attract specific pollinators and ultimately affect plant fitness and persistence [18]. Thus, shifts in resource overlap patterns in the plant-pollinator community may signal whether honey bees could cause eventual wild pollinator (or plant) diversity losses. To examine this, an interaction network approach, considering interaction patterns across the whole community, can be useful [19].

Plant-pollinator interactions can be visualized as bipartite (two trophic-level) networks (Fig 1A), which comprise the set of interactions (“links”) observed between species (“nodes”). These networks represent the pattern of interactions between all species within the plant-pollinator community that are included in the network model [20]. Honey bees may change the structure and function of a plant-pollinator network in two different ways: (i) by adding new interactions (their own interactions with plants), thereby changing network structure without necessarily indicating competitive displacement [21]; and/or by (ii) altering pre-existing plant-pollinator interactions, signaling competitive displacement. Changes to network structure caused by honey bees can be detected by comparing the values for network structural metrics under different honey bee densities (e.g., [19]). To detect changes to the entire network caused by honey bees (i), total networks including honey bee interactions must be examined (e.g., [22, 23]; Fig 1A.i–1a.iii). To detect changes in just the set of wild species interactions (ii), networks not including honey bee interactions must be examined at high versus low honey bee densities [22, 23; Fig 1A.iv]. This is because network metrics are heavily influenced by the most abundant species [24] which honey bees often are [19]. Therefore, changes to the pattern of wild species interactions cannot be detected when honey bee interactions are included in the network model. Here we assess both types of networks, while previous studies have typically tested only total networks (i) [19, 25, 26], but see [27–29]. Importantly, excluding honey bees from a network model (ii) does not exclude honey bees’ ecological effects, because the interactions being visualized in the network occurred when honey bees were present (ii). This is simply a method of building a network model to clearly see wild pollinator interactions, as in [27–29] (S1 File).

[Figure omitted. See PDF.]

Visual representation of the plant-pollinator networks tested (a) and the hypothesized response of each network metric against honey bee abundance (b). a) Plant-pollinator interaction data from each of 19 sampling transects (left), four datasets were created (right): i) a full season all taxa dataset; ii) a mid-season all taxa dataset; iii) a full season bees-only dataset, where only bee interactions, including honey bee interactions, were kept from the full season network; and iv) a without-honey bees dataset, where only honey bee interactions were removed from the full season network. b) Predictions of the effects of increasing honey bee abundance on each network metric, for six metrics of resource overlap and six metrics associated with network stability. Yellow trend lines show predictions for metrics of total network structure, whereas blue trend lines show predictions for metrics calculated from just wild pollinator-plant interactions, excluding honey bee interactions. The rationale and references to support each prediction are in the Supplementary Information. Insect images are from SVG SILH [30].

We used six commonly employed metrics to quantify patterns of resource overlap in plant-pollinator communities [31], from both the pollinator perspective (how pollinator species overlap in what plant species they visit), and the plant perspective (how plant species overlap in what pollinator species visit them). The first metrics provide information on resource overlap from the pollinator perspective. The simplest is “generality”, which is the mean number of plant species visited by each pollinator species [31]. We also calculated pollinator “niche overlap” [31], which is Horn’s information-theory-based index of overlap, ranging from 0 when all species’ interactions within a trophic level are distinct, to 1 when they are identical [32]. Finally, we calculated pollinator “functional complementarity”, a multivariate measure of the dissimilarity of resource use among pollinators. It calculates Euclidean distances between species based on the identity and frequency of their interactions, then clusters the distance matrix to produce a dendrogram, and calculates total branch length of the dendrogram, which provides a univariate expression of this multivariate difference in the set of plant species that each pollinator species interacts with [33]. From the plant perspective, we calculated the equivalent three metrics: “vulnerability” (the mean number of pollinator species visiting each plant species; the name of this metric is from food webs and does not actually refer to vulnerability in plant-pollinator networks), plant “niche overlap” (the Horn’s information theoretic measure of overlap in pollinator species visiting each plant species), and plant “functional complementarity” (the index of multivariate difference in the set of pollinator species visiting each plant species).

We hypothesized that honey bees would compete with wild pollinators [10], excluding some wild pollinator species from some of their normal flower species [9], and thus narrowing their floral niches [7]. We therefore predicted that in networks restricted to wild plant-pollinator interactions, we would see decreases in generality and pollinator niche overlap, and an increase in pollinator functional complementarity, with increasing honey bee abundance (Fig 1B). We also hypothesized that honey bees would dominate interactions with their preferred plant species [8], causing these plant species to be visited by fewer pollinator species [27]. We thus predicted that in networks of just wild plant-pollinator interactions, we would see decreases in vulnerability and plant niche overlap and an increase in plant functional complementarity, with increasing honey bee abundance (Fig 1B).

We also calculated six other commonly reported network metrics to maximize comparability of our results with other studies reporting effects of honey bees on plant-pollinator networks [19, 26–28]. These metrics were: nestedness, modularity, connectance, link density, interaction strength asymmetry, and interaction evenness. These additional metrics have been theoretically linked to network stability [34, 35], but given that it is harder to interpret how changes in these metrics relate specifically to competition, we use them mostly as indicators of network structural change, and restrict our discussion of these metrics to the Supplementary Information. Our predictions for the effects of honey bees on the set of stability-related network metrics are shown in Fig 1B, and explained in detail in the Supplementary Information.

We expected to see the strongest effects of honey bees on native plant-pollinator interactions in midsummer (Fig 1A.ii) and therefore analyzed a midsummer dataset separately from the full season dataset. Honey bee colonies comprise between 10,000 to 50,000 adult individuals, attaining their maximum size [36] and foraging distances [37] in midsummer. Larval numbers peak just before adults, demanding greater pollen collection [38].

We also examined how honey bee abundance affected network metrics for a bees-only version of the full season network (Fig 1A.iii). As most studies reporting effects of honey bees on wild pollinators have only assessed wild bees [3], we wondered if including all pollinating taxa (Hymenoptera, Diptera, Lepidoptera, Hemiptera, Coleoptera) might mask effects on just the wild bee-plant interactions.

Plant-pollinator interaction network studies are complicated by the effects of flower abundance and diversity, which vary spatially (in that pollinators will travel greatly varying distances for plant resources) and temporally (in that different plant species flower at different times of the year), and are not always controlled for. Network metrics also depend on sampling completeness (the proportion of the true interactions that were actually detected; [39], which is not always estimated, making comparability between studies low. Our study controls for effects of all of flower diversity, flower abundance, and sampling effort, and tests experimentally whether changes to plant-pollinator network structure are caused by honey bees changing wild plant-pollinator interactions. Our objectives were to: i) determine if honey bees affect network metrics related to resource overlap and network stability; ii) determine if the effect of honey bees on network metrics is more detectable with higher honey bee colony size; and iii) determine if any effect of honey bees on plant-pollinator network metrics is due to changes in wild plant-pollinator interactions, or due to honey bee-plant interactions.

Materials and methods

Experimental design

In 2019, we arranged three clusters of honey bee hives, set at least 3 km apart, on the Mattheis Research Ranch, which is in mixedgrass prairie rangeland near the town of Brooks, Alberta, Canada (Fig 2; GPS locations in S1 Table). Livestock have grazed these rangelands for a century, and this ranch has used fluctuating but light cattle stocking rates (approximately 0.7 animal unit months ha⁻¹) since 1970 [40]. This region has average high and low temperatures of 18.3°C in July and -12.9°C in January, respectively, and 300 mm of rainfall annually. Honey bees have been present in southern Alberta for over a century, with hive numbers in Alberta increasing by about 4000 hives per year over the past 40 years [41]. Though sporadic bee keeping has occurred in and around the study area in the past, the Eastern Irrigation District, owner of the land surrounding the Mattheis Research Ranch, reported no known apiaries within 19 km of the ranch’s boundary line within the study year or the previous year. Apiaries are required by law to register their hive(s) annually to the Provincial Apiculturist. The hive clusters contained 48, 32, and 16 hives respectively, emulating the range of honey bee densities around commercial apiaries. We established eighteen 30 x 2 m transects at 100 m, 500 m, and 5000 m distances on either side of each cluster of honey bee hives for a total of six replicates (Fig 2). The 100 m and 500 m distances were chosen to achieve high and medium honey bee densities. We chose 5000 m as the minimum distance for sites intended to be without honey bees, because on average they do not tend to travel farther than this distance [42].

[Figure omitted. See PDF.]

Sampling transects are shown as circles and honey bee hives as diamonds. Each apiary contained a cluster of hives (northernmost: 48 hives, central: 32 hives, southernmost: 16 hives). The transect G5000 indicates the new location established when sampling was discontinued at F5000 after July 9^th, 2019. Land outside the Mattheis Ranch is managed by the Eastern Irrigation District. Map shapefiles from Natural Resources Canada under an Open Government license [43]. Map created in QGIS [44].

Between May 28 and August 28, 2019, observers visited each transect and surveyed insect flower visitation almost once per week (weather permitting), for a total of 10 collection rounds. The first transect sampled at the beginning of each day was chosen on a rotating schedule, and subsequent transects were sampled along a route of highest efficiency within that day. During some collection rounds, some transects could not be sampled due to cattle presence or a lack of flowers, resulting in different amounts of sampling for each transect (S1 Table). Transects were visited only on warm (≥ 15°C), sunny days with winds under 50 km/h [45, 46]. We measured wind speed with a Brunton Sherpa (Riverton, Wy, USA). Sampling occurred between 9:30 am and 5:00 pm, when insect flower visits are most frequent [47].

Transects were observed for 30 minutes by two people for a total of 60 person-minutes per transect per collection round (4200 total collection minutes). When large volumes of insects were captured, sampling was paused via a timer for processing and then resumed. All insects that visibly contacted the anthers or stigma of open flowers were collected with a hand net and placed in labelled individual vials, frozen, and identified to species in the lab. Some of these specimens (5%) were identified to morphospecies. When insects observed landing on flowers were missed they were recorded, but these data were excluded because the species-level identifications for these observations were heavily biased toward larger-bodied pollinators. All identifications were completed using a dissecting microscope (Bausch & Lomb Stereozoom, Wild M8, or Leica MZ 125), dichotomous keys, expert taxonomist help, and comparison with reference specimens from the University of Alberta Strickland Museum and University of Calgary Museum of Zoology entomological collections. Following flower visitor sampling, all flowering plants on the transect were identified and their flowers were counted. A list of all references used for the insect and plant species identifications is included in S2 Table. Occasionally, when flowers were not present on the marked transect, we moved the transect ≤ 10 m from its original location to reach flowers near the original transect demarcation. Moving the transect did not change the distance to the hives.

After July 9, 2019, sampling at one 5000 m transect “F5000” was halted after many honey bees were observed (presumably from a feral or unregistered hive, see S1 Fig), since this replicate was intended to be a low-honey bee density transect. A new 5000 m transect, “G5000”, was selected approximately 8000 m away from the northern 48-hive cluster (Fig 2). We treated these transects as distinct site replicates but controlled statistically for the lower number of collections at these by including number of collections as a covariate in our multiple regression models (see below and S1 Table).

This study involved the collection of plant and animal materials from the Mattheis Research Ranch, a research area owned by the University of Alberta, and the Eastern Irrigation District (EID) in Southern Alberta. Permission was granted by both and did not require a field permit number.

Network metric calculation

We pooled the observed insect-flower visitor interactions across the season for each transect, to construct one bipartite network for each transect, for a total of 19 networks (6 x 3 distances from hives, plus the new 5000 m transect “G5000” that was established mid-season).

We calculated the six network metrics described above to analyze honey bee effects on resource use overlap: generality, vulnerability, pollinator niche overlap, plant niche overlap, pollinator functional complementarity, and plant functional complementarity [31]. “Weighted” versions were used, meaning that interactions were weighted by their frequency, rather than just considering interaction presence or absence in calculating the metrics [48].

Additionally, we calculated the weighted versions of six network metrics to analyze the effects of honey bee abundance on network structural metrics related to stability: WNODF, a nestedness metric based on overlap and decreasing fill (“nestedness”) [49], modularity, interaction strength asymmetry, interaction evenness, connectance, and link density. Calculation of these metrics and their interpretation is described in detail in the Supplementary Information. Each metric was calculated using the bipartite package [50]. All statistical analyses were completed using R version 3.2.4. [51]. Additionally, for comparability to [26], in a separate analysis, we pooled the network interactions across the season for the 100 m, 500 m, and 5000 m transects into meta-networks, and calculated each aforementioned network metric for each meta-network and qualitatively compared them in the results.

Statistical analysis

We used linear regression to test the effect of honey bee abundance on each of the six resource-overlap-related and the six stability-related network metrics. We assessed Pearson correlation coefficients between these metrics and found pollinator and plant functional complementarity were highly correlated, as well as link density and vulnerability (S3 Table). We kept all metrics in the analysis, because these were ecological results rather than consequences of the metric definitions and might differ for different datasets. We pooled honey bee abundance (total number of honey bees observed visiting flowers) across the full season for each transect, and divided by number of collection rounds at that transect, to calculate “honey bee abundance” (S1 Fig), which we used as the predictor variable in all full season statistical models. We used honey bee abundance rather than treatment (distance from honey bee hives) as the predictor variable because distance from hives was not as closely associated with honey bee abundance as might have been expected (S1 Fig). Honey bee abundance seemed like a better representation of honey bee “effect” at that transect, as in [25, 28].

We tested the effect of honey bee abundance on each network metric for several datasets: i) each transect’s network from the full season (“Full season all taxa”; Fig 1A.i); ii) each transect’s network for only the mid-season sampling rounds, using honey bee abundance pooled across the mid-season divided by number of mid-season collection rounds at that transect as the predictor variable in this case (“Mid-season all taxa”; Fig 1A.ii described below); and iii) only bee interactions from the full season network (“Full season bees-only”; Fig 1A.iii). In these analyses each network included the observed honey bee interactions, as well as the wild pollinator interactions.

We first ran a simple linear regression (SLR) with honey bee abundance as the only predictor variable for each response variable (network metric) in each dataset. Some of the transects (replicates) were close together, so we plotted the standardized residuals of each SLR to visualize whether similar residuals were close together in space, and to examine whether a special correlation structure was necessary to account for any spatial autocorrelation (S2 Fig).

We then ran generalized least squares (GLS) models with different correlation structures (no correlation, corEcp, corGaus, corSpher, corLin, corRatio) for each response variable (nlme package, [52]). After all models with correlation structures (including no special correlation structure) were run, the model with the lowest AICc value was selected (MuMin package, [53]). The best models for each predictor variable did not include special correlation structures, so simple linear regression models (SLR) were run. The exception was for modularity in the without-honey bees network, so for this metric in this dataset, we ran a GLS model with the corRatio correlation structure.

We then ran an additional multiple regression model (MR) for each response variable in each dataset for which honey bee abundance was significant in the SLR. This was because flower abundance, number of available flower species, and collection effort were positively correlated with honey bee abundance (S4 Table), so we included them as additional predictor variables in this second set of models, to test for an effect of honey bees above and beyond any correlation of honey bee abundance with these other variables [54]. The predictor variables included in each full MR model were honey bee abundance, flower abundance, flower species richness, total number of collection rounds, the interaction between honey bee abundance and flower abundance, as well as the interaction between honey bee abundance and flower species richness. All continuous predictor variables were standardized by subtracting each transect’s value from the mean and dividing by the standard deviation. Then, the full model and all possible simpler models were run, for a total of 21 regression models per response variable (nlme package, [52]), and the best model for each response variable was selected using AICc. In this second set of models, if honey bee abundance was retained after model selection and was significant, it would suggest that there was an effect of honey bees on that response variable that could not plausibly be attributed to other correlated predictor variables [54]; Supplementary Information). If honey bee abundance was significant in the SLR model, but was either not retained or not significant for the MR model, it would suggest that response variable is related to honey bee abundance, but any effects of honey bees cannot be unambiguously attributed to honey bees. They may be due instead to flower abundance or diversity or collection effort, which honey bees were responding to themselves.

The assumptions of normality, homogeneity of variance, and linearity for each response variable were tested. If any assumptions were not met, that response variable was log transformed and model selection was repeated, after which, if assumptions were still not met, a Box-Cox transformation was applied (MASS package, [55]). Transformation did not improve normality or homogeneity of variances for the model for bees-only interaction strength asymmetry, so the results of this model were not interpreted as we could not find a robust modeling method.

To test whether the effects of honey bee abundance on the network metrics differed mid-season, we created a mid-season all taxa dataset (Fig 1A.ii). To determine which sampling dates the mid-season dataset should include, we divided the total honey bee abundance across the season, from all transects pooled, into three roughly equal-length periods: “early”, “mid” and “late”, based on natural breaks in the abundances of honey bees (S3 Fig). We then tested the effects of honey bee abundance on all the same response variables as above, using the procedure outlined above, for this mid-season all taxa dataset. Because some transects’ networks were too species-poor in the mid-season to obtain accurate values (A100, B100, B500, E100, F5000), they were removed from this analysis.

Finally, we created a reduced network from the full season all taxa dataset where all honey bees were removed from the network, to look at the network metrics for wild pollinators only (Fig 1A.iv; [28]). In this “without-honey bees” dataset we ran SLR and MR models that tested the effect of honey bee abundance and all predictor variables, respectively, on all the same response variables. If honey bee abundance was a significant predictor in these without-honey bees models, it would indicate that honey bees altered how wild pollinators interacted; if it was not, it would indicate that the changes to that network metric caused by honey bee abundance in models for our other datasets were just due to the honey bee node and links contributing to the network structure.

This analysis involved interpreting 65 separate models, so to maintain Type I Error at 0.05, a Bonferroni-Holm correction was used in assessing significance of the final models (at P = 0.00077). However, all P-values are presented to allow a less conservative consideration of uncorrected P-values. To estimate the interaction sampling completeness, we divided the raw interaction richness across all transects by the Chao1 estimated total interaction richness across all transects (SpadeR, [56]). This calculated the proportion of the estimated total number of unique interactions that our sampling detected (as in [39]). We repeated this at the transect level for the full season all taxa, mid-season all taxa, and full season bees-only datasets.

Results

Diversity of pollinators, plants, and interactions in the grassland plant-pollinator community

A total of 281 pollinator species and 37 plant species were identified in 1,814 interactions, with 654 unique species pair interactions. Of these interactions, 425 were only observed once. The Chao1 estimated true total number of unique interactions was 1,779 interactions (95% confidence interval: 1500–2148 interactions), meaning that we observed 31–44% of the estimated total number of plant-pollinator interactions occurring in this study region (S4 Fig). We also calculated the Chao1 estimated true number of unique interactions at each transect, for the full season, mid-season, and full season bees-only datasets (S5 Table). At 100 m transects, honey bees accounted for 29.4% of the total interactions, at 500 m transects, 14.9% of the total interactions, and at 5000 m transects, 1% of the total interactions. Bees (including honey bees) accounted for 872 (48.1%) of the total interactions. Of the remaining interactions, flies accounted for 494 (27.2%), butterflies for 130 (7.2%), ants for 104 (5.7%), beetles for 100 (5.5%), wasps for 62 (3.4%), true bugs for 27 (1.5%), and day-flying moths for 25 (1.4%). In the mid-season all taxa dataset, at 100 m transects, honey bees accounted for 49.6% of the interactions, at 500 m transects, 25.2% of the interactions, and at 5000 m transects, 2.4% of the interactions. The meta-networks of pooled interactions for each transect are shown in Fig 3. All species identifications are listed in S6 (pollinators) and S7 (plants) Tables.

[Figure omitted. See PDF.]

Data from transects 100 m (top), 500 m (middle) and 5000 m (bottom) from honey bee hives are pooled across the season. The upper row depicts the pollinator species, coloured by their given taxonomic “group” or Order (see S6 Table for the list of pollinator species). “Moths” are Lepidoptera excluding butterflies (Papilionoidea). “Wasps” are Aculeata excluding ants (Formicoidea) and bees (Anthophila). The bottom row (in black) depicts the plant species (see S7 Table for the list of plant species). The width of each upper and lower bar represents the relative frequency of interactions observed for that species. The interactions between plants and wild pollinators are represented by blue lines; width indicates the frequency of the interaction. Honey bees and their interactions are indicated in red.

Effect of honey bees on the structure of networks including honey bee interactions

In the full season all taxa dataset, in the SLR models (testing honey bee abundance as the only predictor against each network metric), increasing honey bee abundance was associated with significant changes in two network metrics related to resource overlap. Pollinator functional complementarity (t = 5.92, P < 0.0001) and plant functional complementarity (t = 6.39, P < 0.0001) significantly increased with honey bee abundance, indicating decreasing similarity in resource use between species within a trophic level with increasing honey bee abundance (Fig 4A, Table 1A). Increasing honey bee abundance was also associated with a significant decrease in interaction evenness (t = -4.33, P = 0.0005; Fig 4A, Table 1A). For the mid-season all taxa dataset (Fig 4A, Table 1B), and the full season bees-only dataset (Fig 4A, Table 1C), the results were very similar. All reported P-values are significant after the Bonferroni-Holm correction.

[Figure omitted. See PDF.]

(a) and without-honey bees networks (b). a) Relationships between honey bee abundance and the network metrics: (i) pollinator functional complementarity, (ii) plant functional complementarity, and (iii) interaction evenness that were significantly related to honey bee abundance across three datasets (full season all taxa, mid-season all taxa, and full season bees-only). Solid regression lines indicate significant relationships with Bonferroni-Holm correction, whereas dotted lines indicate insignificant relationships. b) Relationships between honey bee abundance and the network metrics: (i) pollinator niche overlap, (ii) pollinator functional complementarity, (iii) plant functional complementarity, (iv) weighted connectance, and (v) interaction evenness in the without-honey bees dataset. Dashed regression lines indicate non-significant relationships with Bonferroni-Holm correction. Although (i) and (iv) appear to contain trends, these cannot be unambiguously attributed to honey bee abundance, as they can be explained by flower community variables and collection effort (Table 2D) that were somewhat correlated with honey bee abundance (S4 Table).

[Figure omitted. See PDF.]

[Figure omitted. See PDF.]

In the second set of MR models (including all predictor variables), flower species richness and total number of collections best explained most response variables (Table 2A). However, increasing honey bee abundance was still associated with significant increases in pollinator functional complementarity (t = 5.51, P < 0.0001) and plant functional complementarity (t = 5.75, P < 0.0001), and a significant decrease in interaction evenness (t = -4.33, P = 0.0005) (Table 2A), indicating that these effects on plant-pollinator network structure could be attributed to honey bees, rather than to flower community variables that were correlated with honey bee abundance. For the mid-season all taxa dataset (Table 2B) and the full season bees-only dataset (Table 2C), the results were very similar, except that in the full season bees-only dataset, pollinator functional complementarity was not significantly increased by honey bees after the Bonferroni-Holm correction, even though it still had a very low P-value (t = 3.57, P = 0.0026).

Effect of honey bees on the without-honey bees network structure

In contrast to all of the above results, the full season without-honey bees dataset with only honey bee abundance as a predictor variable (SLR models) showed that honey bees did not significantly affect pollinator functional complementarity, plant complementarity, or interaction evenness, or any other network metric (Fig 4B, Table 2D). Although some P-values were low in the SLR models (Table 1D), in the MR models, honey bee abundance was removed during model selection in every case (Table 2D), suggesting that the low P-values in the SLR models could be explained by correlated effects of the flower variables or collection effort better than they could be explained by honey bee abundance. This disappearance of any significant effect of honey bee abundance on any network metrics when honey bees were not considered as part of the network was not an artifact of having decreased the network size by removing the honey bee interactions (see comparison of all network sizes in S8 Table).

Discussion

As predicted, total plant-pollinator network structure (including honey bee interactions), changed with increasing honey bee abundance, as reported in other studies [19, 25, 26, 28]. Across three datasets (full season all taxa, mid-season all taxa, and full season bees-only), the changes due to honey bees (rather than potentially due to flower community variables correlated with honey bee abundance) were that pollinator functional complementarity and plant functional complementarity significantly increased, and interaction evenness significantly decreased as honey bee abundance increased. The effect sizes were highest in the mid-season all taxa dataset, when honey bee abundances were highest (compare regression coefficients in Table 2A to partial regression coefficients in Table 2B). The exception was that pollinator functional complementarity did not change significantly in the bees-only dataset (Table 2C). When looking at the meta-networks (Fig 3), which had interactions pooled across transects by distance from hives, similar to [26], rather than pooled by transect, the results were qualitatively similar (S9 Table), suggesting that our results are robust across different analysis methods.

However, when considering only wild plant-pollinator interactions in the without-honey bees dataset, honey bee abundance (which varied from high to almost zero between these sites) did not affect any of the network metrics (Table 2D), suggesting that all changes to network structure associated with honey bee abundance in this study were caused by the honey bees’ interactions contributing to the network structure, rather than by honey bees altering the interactions of wild pollinators and plants. Thus, adding honey bees had no detectable effect on the network structure of the wild plant-pollinator interactions in this region, whether considering metrics of resource overlap, or metrics associated with network stability.

Magrach et al. [19], who reported changes in plant-pollinator network structure between orange groves in Spain in which honey bees accounted for 72% vs. 38% of interactions, found similar results to our with-honey bee networks: honey bees decreased interaction evenness and increased functional complementarity in their system. The increase in functional complementarity indicates increasing dissimilarity in the interaction composition of pollinators and plants in the network, which can be thought of as a reduction of their resource use overlap [33]. Although this metric accounts for the interactions of all species in the community, its value can change in response to the addition of an abundant generalist species’ interactions, such as the honey bee, because an abundant generalist species’ interaction identities and frequencies are so different from those of less abundant and generalist species (corresponding to a long dendrogram branch length). However, Magrach et al. [19] did not test with- versus without-honey bee networks for honey bee abundance effects on these metrics. Although their with-honey bee network results are similar to ours, we found that in our without-honey bees network, the effect of honey bee abundance disappeared, suggesting no effect of honey bees on wild plant-pollinator interactions.

Similarly, Valido et al. [26], who tested changes in network structure before vs. after honey bees were added to one site in Tenerife, Canary Islands in three consecutive years, did not remove the honey bee node from the networks and re-test effects of honey bee abundance on just the wild pollinator interactions. Additionally, they did not control for floral abundance or flower species richness as potentially correlated predictor variables. Furthermore, we had similar sampling effort, range of honey bee abundances tested, and estimated sampling completeness to both [19] and [26]. Taken together, this suggests that if more potentially correlated variables had been controlled for, these two studies may also have found a non-significant effect of honey bees on wild plant-pollinator network structure. Alternatively, it may be that honey bees shift wild pollinator interactions in other systems more than in the grassland system where our study took place. Lázaro et al. [28], similar to our study, tested honey bee effects on total network structure, as well as on the network structure of just wild pollinator interactions, and found larger effects on total network structure. They reported that the effects of honey bee abundance on wild bee network metrics that they observed were caused by changes in wild bee species richness, whereas in our study, wild pollinator species richness did not change with increasing honey bee abundance.

The most unexpected result from our study was that we found no effect of honey bees on nestedness, whether honey bee interactions were included or not (Table 1A–1D). This contrasts with several studies that have reported increases in network nestedness when honey bees were present [21–23, 28, 57]. However, all these studies except [28] tested honey bee effects on “binary” nestedness, calculated without considering interaction frequencies, whereas we used “weighted” nestedness, which incorporates interaction frequency, because [58] showed that changes to weighted nestedness are a better indicator of shifts in resource partitioning. Our finding of no change to weighted nestedness associated with increasing honey bee abundance is further evidence that honey bees did not change patterns of resource partitioning by wild pollinators or plants [58]. Our findings for the other stability-related metrics for the with-honey bees networks are similar to the findings of other studies (see Supplementary Information).

Our failure to find evidence that honey bees shift resource use by wild pollinators, or shift visitation of plant species, may indicate that honey bees do not compete with wild pollinators in this system, due to abundant floral resources that are not limiting. Additionally, honey bees interacted frequently with non-native plants. Alberta’s grasslands are diverse and dominated by native plant species, but several exotic plants exist in the region; sweet clover (Melilotus spp.), sow-thistle (Sonchus spp.), and cicer milk-vetch (Astragalus cicer), in particular, exist in scattered patches throughout the study region. We observed that these exotic species were highly attractive to pollinators of all types. In their network study, Bendel et al. [59] determined that honey bees preferred exotic plants over native plants. We also observed that small-flowered and small-stemmed plants were less appealing to larger-bodied pollinators. This preference may have prevented honey bees from outcompeting wild pollinators for resources, preserving native plant-pollinator interactions.

The plant-pollinator network in this system appears to accommodate honey bee interactions, integrating them into the system without major changes, despite honey bees taking a central role (Fig 3). Honey bees may spread pathogens and parasites to some wild bee species via shared flowers [60], or cause other changes to wild plant or pollinator fitness undetectable from looking at changes to flower visitation [18]. Furthermore, we did not detect all occurring plant-pollinator interactions at our transects. Our estimated sampling completeness was 37% overall (see S4 Table), suggesting that we missed many rare interactions. Rare interactions typically involve rare species, which may be particularly sensitive to competition from honey bees [5]. However, understanding honey bee effects on the subset of detected interactions is useful because the most detectable interactions are the most frequent, and are those that contribute most to ecosystem function [61].

Additionally, our study was conducted in only one season, and plant-pollinator interactions are known to be highly variable year to year [62]. This variability is largely because of species and interaction turnover, however, and network structure is expected to be more consistent between years [62], so our conclusions about network structure might be expected to hold in other years. The sampled year 2019 was drier than the 20-year average (S9 Table), meaning floral resources should have been more limiting than in a normal year [63], and therefore we should have been more likely to see competitive effects of honey bees on native pollinators in 2019. We also sampled in 2018, another very dry year (S10 Table), and though the data were not of high enough quality to publish, they indicated the same results.

Conclusions

Honey bees have been frequently associated with wild pollinator decline, and because regulating honey bee densities might assist wild pollinator conservation, it is worth establishing whether honey bees are a significant threat to wild pollinators. Our study took place in a part of Canada with high bee diversity, in a dry year when floral resources for pollinators should have been more limiting than usual. In this setting we found no effects of honey bees on wild pollinator-plant network metrics, suggesting that even at high densities, honey bees do not detectably cause shifts in wild pollinator foraging habits in this region. Future studies should control for the effects of flower abundance and diversity on plant-pollinator network structure, and distinguish the effects of honey bee interactions on total network structure from effects of honey bees on just the wild pollinator interactions with plants.

Supporting information

S1 File.

https://doi.org/10.1371/journal.pone.0287332.s001

(DOCX)

S1 Fig. Abundance of honey bees caught visiting flowers (full season all taxa dataset), pooled across the full season per transect, with transects ordered by increasing honey bee abundance, and coloured by distance from bee hives.

In the transect names, 100 indicates 100 m, 500 indicates 500 m, and 5000 indicates 5000 m distances from hives.

https://doi.org/10.1371/journal.pone.0287332.s002

(TIFF)

S2 Fig. Visual representation of the spatial autocorrelation of residuals from the model for each response variable for the full season all taxa dataset, for which honey bee abundance was a significant predictor (before Bonferroni-Holm correction) when it was the only predictor in the model, for full season dataset.

Each panel is a map of the transect locations in space (compare to Fig 2). Circles indicate the residual from each transect, with circle size proportional to residual size (smaller circle = better model fit for that transect). Colour indicates the sign of the residual; blue shows values lower than 0 and red values higher than 0. Circles close together in space having the same colour and size would indicate that spatial autocorrelation might be a problem, in which case a special correlation structure would be likely to be selected during model selection.

https://doi.org/10.1371/journal.pone.0287332.s003

(TIFF)

S3 Fig. Total abundance of honey bees caught visiting flowers, with data pooled across all transects.

The entire season was split into three: collection rounds 1–4 represented “early” season (May 28^th-July 7^th), 5–7 represented “mid” season (July 8^th-July 31^st), and 8–10 represented “late” season (August 1^st-August 28^th).

https://doi.org/10.1371/journal.pone.0287332.s004

(TIFF)

S4 Fig. Rarefaction curve for interaction richness.

This figure shows, for a given number of re-sampled interactions from our full season all taxa dataset (x-axis), the mean number of unique interactions detected (y-axis). There was a total observed interaction richness of 654, and Chao1 estimated 1,779 interactions (95% confidence interval: 1500–2148 interactions), meaning 31–44% of the estimated interactions were observed.

https://doi.org/10.1371/journal.pone.0287332.s005

(TIFF)

S1 Table. Longitude and latitude for each hive location and sampling transect, and collection effort at each transect.

The northernmost, central, and southernmost hive locations are listed as Bee48, Bee32, and Bee16 respectively, indicating their number of hives. Each transect is indicated by its treatment (100 m, 500 m, or 5000 m distance from a hive location). Letters indicate each replicated set of distances from hives (See Fig 1). G5000 indicates the new location sampled at once sampling at F5000 was discontinued mid-season.

https://doi.org/10.1371/journal.pone.0287332.s006

(DOCX)

S2 Table. List of references and resources used in species identifications.

Full citations are listed in the References section of the Supporting Information.

https://doi.org/10.1371/journal.pone.0287332.s007

(DOCX)

S3 Table. Pearson correlations between the response variables (network metrics).

Each row and column header refers to the following metrics, respectively: G, generality; V, vulnerability; Poll. NO, pollinator niche overlap; Pl. NO, plant niche overlap; Poll. FC, pollinator functional complementarity; Pl. FC, plant functional complementarity; WNODF, weighted nestedness based on overlap and decreasing fill; M, modularity; WC, weighted connectance; LD, link density; ISA, interaction strength asymmetry; IE, interaction evenness.

https://doi.org/10.1371/journal.pone.0287332.s008

(DOCX)

S4 Table. Pearson correlations between predictor variables.

“Collections” refers to the number of collection rounds at each transect, “Honey bee abundance” refers to honey bee abundance at each transect, “Flower species” refers to the number of flowering species at each transect, and “Flower abundance” refers to the number of individual flowers at each transect.

https://doi.org/10.1371/journal.pone.0287332.s009

(DOCX)

S5 Table. Percent sampling completeness and 95% confidence intervals, calculated as the number of unique interactions observed, divided by the Chao1 estimated total number of unique interactions, and multiplied by 100, for the full season all taxa, mid-season all taxa, and full season bees-only datasets.

The “Lower 95%” column indicates the observed interactions divided by the lowest Chao1 estimated total number of true interactions, and the “Upper 95%” column indicates the observed interactions divided by the highest Chao1 estimated total number of true interactions (which is why the “Lower” values are higher). Shaded cells denote networks that had too few species for some network metrics to be meaningful, so these replicates were excluded from analysis for the metrics: generality, vulnerability, plant and pollinator niche overlap, and nestedness.

https://doi.org/10.1371/journal.pone.0287332.s010

(DOCX)

S6 Table. Identifications of insect pollinators to species-level or morphospecies level.

Morphospecies identifications are listed by “[Genus] sp. #”. Some species could not be differentiated between genera, and so both genera are listed along with the epithet “sp”. Specimens listed beside “cf” (confer, meaning compare with) are specimens that were damaged or for which taxonomic keys are insufficient, and these were compared to other specimens to determine identification. Numbers of each (morpho)species are given for each distance from hives, despite the fact that we used honey bee abundance, rather than distance from hive as the predictor variable in our analyses.

https://doi.org/10.1371/journal.pone.0287332.s011

(DOCX)

S7 Table. Identifications of flowering species from each distance from honey bee hives to species level.

https://doi.org/10.1371/journal.pone.0287332.s012

(DOCX)

S8 Table. Network size (number of interactions) for each network that was analyzed from each transect.

Shaded cells denote networks that had too few species for some network metrics to be meaningful, so these replicates were excluded from analysis for the metrics: generality, vulnerability, plant and pollinator niche overlap, and nestedness. Mean network size and range in network size are shown for each dataset.

https://doi.org/10.1371/journal.pone.0287332.s013

(DOCX)

S9 Table. Comparison of the full season meta-network metrics calculated for networks created by pooling data from all transects at each distance from hives, over the whole season.

There was a general increase in each network metric with distance to hives, with the exceptions of plant niche overlap, functional complementarity, and generality. Bolded values indicate results that were unexpected based on the literature and/or expectations if competition between honey bees and wild pollinators is occurring. Positive (+) indicates a positive correlation between honey bee abundance and the metric, while negative (-) indicates a negative correlation, with distance from hives as a proxy for honey bee abundance (even though this was not a perfect proxy, as can be seen in S1 Fig).

https://doi.org/10.1371/journal.pone.0287332.s014

(DOCX)

S10 Table. Historical data obtained from the Alberta Climate Information Service (ACIS) Verger AGCM weather station in Southern Alberta, located on the University of Alberta’s Mattheis Research Ranch.

https://doi.org/10.1371/journal.pone.0287332.s015

(DOCX)

Acknowledgments

We thank Brittany Wingert, Irene Jimenez Roncancio, James Glasier, Greg Pohl, Lincoln Best, and the University of Calgary for their assistance in the identification of specific insect taxa. We also thank assistants Alexandra Johnson and Janelle Goodine, and volunteers Connor Nelson, Olivia DeBourcier, Zachary Roote, Ferf Brownoff, Sara Peterson, Olivia Hrehoruk, Rykkar Jackson, and Victoria Dubord for all their support on this project. We thank Edwin and Ruth Mattheis for donating the Mattheis Research Ranch to the University of Alberta, and Marcel Busz and Lisa Raatz for assistance with logistics at the ranch. We also thank the Eastern Irrigation District for allowing the use of their land. Honey bees were generously supplied by Pankratz Beekeeping. We thank Cameron Carlyle and Mark Poesch for helpful discussion.

Citation: Worthy SH, Acorn JH, Frost CM (2023) Honey bees (Apis mellifera) modify plant-pollinator network structure, but do not alter wild species’ interactions. PLoS ONE 18(7): e0287332. https://doi.org/10.1371/journal.pone.0287332

About the Authors:

Sydney H. Worthy

Roles: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Visualization, Writing – original draft, Writing – review & editing

E-mail: [email protected]

Current address: Parks Department, City of Saskatoon, Saskatoon, Canada

Affiliation: Department of Renewable Resources, University of Alberta, Edmonton, Alberta, Canada

ORICD: https://orcid.org/0000-0002-6696-2283

John H. Acorn

Roles: Conceptualization, Data curation, Investigation, Methodology, Supervision, Writing – review & editing

Affiliation: Department of Renewable Resources, University of Alberta, Edmonton, Alberta, Canada

Carol M. Frost

Roles: Conceptualization, Funding acquisition, Investigation, Methodology, Resources, Supervision, Writing – review & editing

Affiliation: Department of Renewable Resources, University of Alberta, Edmonton, Alberta, Canada

ORICD: https://orcid.org/0000-0003-3805-5528

References

1. Grixti JC, Wong LT, Cameron SA, Favret C. Decline of bumble bees (Bombus) in the North American Midwest. Biological Conservation. 2009 Jan;142(1):75–84.

2. Iwasaki JM, Hogendoorn K. Mounting evidence that managed and introduced bees have negative impacts on wild bees: an updated review. Current Research in Insect Science. 2022;2:100043. pmid:36003276

3. Mallinger RE, Gaines-Day HR, Gratton C. Do managed bees have negative effects on wild bees?: A systematic review of the literature. Raine NE, editor. PLoS ONE. 2017 Dec 8;12(12):e0189268.

4. Potts SG. The assessment report on pollinators, pollination and food production: summary for policymakers. Bonn, Germany: Secretariat of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services; 2016.

5. Waser NM, Ollerton J, editors. Plant-pollinator interactions: from specialization to generalization. Chicago: University of Chicago Press; 2006. 445 p.

6. Badano EI, Vergara CH. Potential negative effects of exotic honey bees on the diversity of native pollinators and yield of highland coffee plantations. Agricultural and Forest Entomology. 2011 Nov;13(4):365–72.

7. Forup ML, Memmott J. The relationship between the abundances of bumblebees and honeybees in a native habitat. Ecol Entomol. 2005 Feb;30(1):47–57.

8. Nielsen A, Reitan T, Rinvoll AW, Brysting AK. Effects of competition and climate on a crop pollinator community. Agriculture, Ecosystems & Environment. 2017 Aug;246:253–60.

9. Hudewenz A, Klein A. Red mason bees cannot compete with honey bees for floral resources in a cage experiment. Ecol Evol. 2015 Nov;5(21):5049–56. pmid:26640681

10. Cane JH, Tepedino VJ. Gauging the Effect of Honey Bee Pollen Collection on Native Bee Communities. Conservation Letters. 2017;10(2):205–10.

11. Morse DH. Resource Partitioning in Bumble Bees: The Role of Behavioral Factors. Science. 1977 Aug 12;197(4304):678–80. pmid:17776271

12. Taggar AK, McGrath E, Despland E. Competition between a native and introduced pollinator in unmanaged urban meadows. Biol Invasions. 2021 Jun;23(6):1697–705.

13. Heinrich B. Resource Partitioning Among Some Eusocial Insects: Bumblebees. Ecology. 1976 Aug;57(5):874–89.

14. Fontaine C, Collin CL, Dajoz I. Generalist foraging of pollinators: diet expansion at high density: Flexibility of pollinator generalization. Journal of Ecology. 2008 Sep;96(5):1002–10.

15. Wignall VR, Brolly M, Uthoff C, Norton KE, Chipperfield HM, Balfour NJ, et al. Exploitative competition and displacement mediated by eusocial bees: experimental evidence in a wild pollinator community. Behav Ecol Sociobiol. 2020 Dec;74(12):152.

16. Cervantes‐Loreto A, Ayers CA, Dobbs EK, Brosi BJ, Stouffer DB. The context dependency of pollinator interference: How environmental conditions and co‐foraging species impact floral visitation. Irwin R, editor. Ecology Letters. 2021 Jul;24(7):1443–54.

17. Montero‐Castaño A, Vilà M. Influence of the honeybee and trait similarity on the effect of a non‐native plant on pollination and network rewiring. Ayasse M, editor. Funct Ecol. 2017 Jan;31(1):142–52.

18. Aizen MA, Morales CL, Vázquez DP, Garibaldi LA, Sáez A, Harder LD. When mutualism goes bad: density-dependent impacts of introduced bees on plant reproduction. New Phytol. 2014 Oct;204(2):322–8.

19. Magrach A, González-Varo JP, Boiffier M, Vilà M, Bartomeus I. Honeybee spillover reshuffles pollinator diets and affects plant reproductive success. Nat Ecol Evol. 2017 Sep;1(9):1299–307. pmid:29046536

20. Delmas E, Besson M, Brice MH, Burkle LA, Dalla Riva GV, Fortin MJ, et al. Analysing ecological networks of species interactions. Biological Reviews. 2019;94(1):16–36. pmid:29923657

21. Geslin B, Gauzens B, Baude M, Dajoz I, Fontaine C, Henry M, et al. Massively Introduced Managed Species and Their Consequences for Plant–Pollinator Interactions. Advances in Ecological Research. 2017;57:147–99.

22. Santos GM, Aguiar CML, Genini J, Martins CF, Zanella FCV, Mello MAR. Invasive Africanized honeybees change the structure of native pollination networks in Brazil. Biol Invasions. 2012 Nov;14(11):2369–78.

23. Giannini TC, Cordeiro GD, Freitas BM, Saraiva AM, Imperatriz-Fonseca VL. The Dependence of Crops for Pollinators and the Economic Value of Pollination in Brazil. Journal of Economic Entomology. 2015 Jun 1;108(3):849–57. pmid:26470203

24. Vázquez DP, Melián CJ, Williams NM, Blüthgen N, Krasnov BR, Poulin R. Species abundance and asymmetric interaction strength in ecological networks. Oikos. 2007 Jul;116(7):1120–7.

25. Prendergast KS, Ollerton J. Spatial and temporal scale of analysis alter conclusions about the effects of urbanisation on plant–pollinator networks. Arthropod-Plant Interactions. 2022 Dec;16(6):553–65.

26. Valido A, Rodríguez-Rodríguez MC, Jordano P. Honeybees disrupt the structure and functionality of plant-pollinator networks. Sci Rep. 2019 Dec;9(1):4711. pmid:30886227

27. Ropars L, Dajoz I, Fontaine C, Muratet A, Geslin B. Wild pollinator activity negatively related to honey bee colony densities in urban context. Blenau W, editor. PLoS ONE [Internet]. 2019 Sep 12 [cited 2021 Jul 1];14(9). Available from: pmid:31513663

28. Lázaro A, Müller A, Ebmer AW, Dathe HH, Scheuchl E, Schwarz M, et al. Impacts of beekeeping on wild bee diversity and pollination networks in the Aegean Archipelago. Ecography. 2021 Sep;44(9):1353–65.

29. Ropars L, Affre L, Thébault É, Geslin B. Seasonal dynamics of competition between honey bees and wild bees in a protected Mediterranean scrubland. Oikos [Internet]. 2022 Apr [cited 2023 Feb 3];2022(4). Available from: https://onlinelibrary.wiley.com/doi/10.1111/oik.08915

30. SVG SILH. SVG SILH [Internet]. SVG SILH. 2022 [cited 2022 Apr 29]. Available from: https://svgsilh.com/

31. Dormann CF, Gruber B, Frueund J. Introducing the bipartite package: Analysing Ecological Networks. R News. 2008;8(2):8–11.

32. Horn HS. Measurement of “Overlap” in Comparative Ecological Studies. The American Naturalist. 1966 Sep;100(914):419–24.

33. Devoto M, Bailey S, Craze P, Memmott J. Understanding and planning ecological restoration of plant-pollinator networks. Ecology Letters. 2012;15(4):319–28. pmid:22251948

34. Okuyama T, Holland JN. Network structural properties mediate the stability of mutualistic communities. Ecol Letters. 2008 Mar;11(3):208–16. pmid:18070101

35. Thébault E, Fontaine C. Stability of Ecological Communities and the Architecture of Mutualistic and Trophic Networks. Science. 2010 Aug 13;329(5993):853–6. pmid:20705861

36. Canadian Honey Council. Bee Facts [Internet]. Canadian Honey Council. 2018 [cited 2021 Jun 18]. Available from: https://honeycouncil.ca/bee-facts/

37. Couvillon MJ, Schürch R, Ratnieks FLW. Waggle Dance Distances as Integrative Indicators of Seasonal Foraging Challenges. Smagghe G, editor. PLoS ONE [Internet]. 2014 Apr 2 [cited 2022 Apr 28];9(4). Available from: pmid:24695678

38. Torres DJ, Ricoy UM, Roybal S. Modeling Honey Bee Populations. Rueppell O, editor. PLoS ONE [Internet]. 2015 Jul 6 [cited 2023 May 6];10(7). Available from: pmid:26148010

39. Chacoff NP, Vázquez DP, Lomáscolo SB, Stevani EL, Dorado J, Padrón B. Evaluating sampling completeness in a desert plant-pollinator network: Sampling a plant-pollinator network. Journal of Animal Ecology. 2012 Jan;81(1):190–200.

40. Thomas BW, Gao X, Zhao M, Bork EW, Hao X. Grazing altered carbon exchange in a dry mixedgrass prairie as a function of soil texture. Can J Soil Sci. 2017 Dec 18;CJSS-2017-0052.

41. Laate EA. Economics of beekeeping in Alberta 2016. Economics Section, Economics and Competitiveness Branch, Alberta Agriculture and Forestry; 2017.

42. Hagler JR, Mueller S, Teuber LR, Machtley SA, Van Deynze A. Foraging Range of Honey Bees, Apis mellifera, in Alfalfa Seed Production Fields. Journal of Insect Science. 2011 Nov;11(144):1–12.

43. Government of Canada. Map Projects Inventory [Internet]. Natural Resources Canada; 2023 [cited 2023 Apr 10]. Available from: https://open.canada.ca/data/en/dataset/f5f2db55-31e4-42fb-8c73-23e1c44de9b2

44. QGIS Development Team. QGIS Geographic Information System [Internet]. QGIS Association; 2023. Available from: https://www.qgis.org

45. Pasek JE. Influence of wind and windbreaks on local dispersal of insects. Agriculture, Ecosystems & Environment. 1988;22–23:539–54.

46. Sanderson RA, Goffe LA, Leifert C. Time-series models to quantify short-term effects of meteorological conditions on bumblebee forager activity in agricultural landscapes. Agricultural and Forest Entomology. 2015;17(3):270–6.

47. Lefebvre V, Fontaine C, Villemant C, Daugeron C. Are empidine dance flies major flower visitors in alpine environments? A case study in the Alps, France. Biology Letters. 2014;10(11). pmid:25376804

48. Bersier LF, Banašek-Richter C, Cattin MF. Quantitative descriptors of food-web matrices. Ecology. 2002 Sep;83(9):2394–407.

49. Almeida-Neto M, Guimarães P, Guimarães PR, Loyola RD, Ulrich W. A consistent metric for nestedness analysis in ecological systems: reconciling concept and measurement. Oikos. 2008 Aug;117(8):1227–39.

50. Dormann CF, Frueund J, Bluethgen N, Gruber B. Indices, graphs and null models: analyzing bipartite ecological networks. The Open Ecology Journal. 2009;2:7–24.

51. R Core Team. R: A language and environment for statistical computing. [Internet]. R Foundation for Statistical Computing; 2020 [cited 2021 Aug 2]. Available from: https://www.r-project.org/

52. Jose P, Douglas B, Saikat D, Sarkar D, R Core Team. nlme: Linear and Nonlinear Mixed Effects Models. 2020.

53. Bartoń K. MuMIn: Multi-model inference. [Internet]. [cited 2021 Aug 2]. Available from: https://r-forge.rproject.org/R/?group_id=346

54. Morrissey MB, Ruxton GD. Multiple Regression Is Not Multiple Regressions: The Meaning of Multiple Regression and the Non-Problem of Collinearity. Philosophy, Theory, and Practice in Biology [Internet]. 2018 Jun [cited 2022 Apr 28];10(20220112). Available from: http://hdl.handle.net/2027/spo.16039257.0010.003

55. Venables WN, Ripley BD. Modern Applied Statistics with S [Internet]. 4th ed. Springer; 2002 [cited 2021 Aug 2]. Available from: http://www.stats.ox.ac.uk/pub/MASS4/

56. Chao A, Ma KH, Hsieh TC, Chiu CH. SpadeR: Species-Richness Prediction and Diversity Estimation with R [Internet]. 2016 [cited 2021 Aug 2]. Available from: https://cran.r-project.org/package=SpadeR

57. Traveset A, Heleno R, Chamorro S, Vargas P, McMullen CK, Castro-Urgal R, et al. Invaders of pollination networks in the Galápagos Islands: emergence of novel communities. Proc R Soc B. 2013 May 7;280(1758):20123040.

58. Staniczenko PPA, Kopp JC, Allesina S. The ghost of nestedness in ecological networks. Nat Commun. 2013 Jun;4(1):1391. pmid:23340431

59. Bendel CR, Kral‐O’Brien KC, Hovick TJ, Limb RF, Harmon JP. Plant–pollinator networks in grassland working landscapes reveal seasonal shifts in network structure and composition. Ecosphere [Internet]. 2019 Jan [cited 2021 Aug 2];10(1). Available from: https://onlinelibrary.wiley.com/doi/10.1002/ecs2.2569

60. Graystock P, Goulson D, Hughes WOH. Parasites in bloom: flowers aid dispersal and transmission of pollinator parasites within and between bee species. Proc R Soc B. 2015 Aug 22;282(1813):20151371. pmid:26246556

61. Winfree R W. Fox J, Williams NM, Reilly JR, Cariveau DP. Abundance of common species, not species richness, drives delivery of a real-world ecosystem service. Shea K, editor. Ecol Lett. 2015 Jul;18(7):626–35. pmid:25959973

62. Zografou K, Swartz MT, Tilden VP, McKinney EN, Eckenrode JA, Sewall BJ. Stable generalist species anchor a dynamic pollination network. Ecosphere [Internet]. 2020 Aug [cited 2022 Apr 30];11(8). Available from: https://onlinelibrary.wiley.com/doi/10.1002/ecs2.3225

63. Phillips BB, Shaw RF, Holland MJ, Fry EL, Bardgett RD, Bullock JM, et al. Drought reduces floral resources for pollinators. Glob Change Biol. 2018 Jul;24(7):3226–35. pmid:29652102