Contrasting effects of exotic plant invasions and managed honeybees on plant–flower visitor interactions

To explore how a highly invasive plant species (Buddleja davidii Franch.), managed honeybees and flower diversity affected plant–flower visitor interactions over the whole elevational range distribution of the exotic plant.

In addition to the anthropogenic drivers described above, natural processes can also influence plant-pollinator interactions. In mountains, the degree of specialization of plant-pollinator networks usually decreases with elevation (Hoiss, Krauss, & Steffan-Dewenter, 2015) due to a more unpredictable and unfavourable environment that narrows resource accessibility at extreme cold temperatures (Adedoja, Kehinde, & Samways, 2018;Johnson & Steiner, 2000).
This reduction in network specialization is often linked to the replacement of specialized pollinator guilds (e.g., many wild bees) with groups with wider niche breadth, such as flies and wasps (Benadi et al., 2014;Chacoff et al., 2012;Ramos-Jiliberto et al., 2010).
However, elevation and plant diversity are often related in mountains, with native plants showing a hump-shaped relationship with elevation (Guo et al., 2013) and exotic plants typically decreasing with increasing elevation (Alexander et al., 2011;Marini et al., 2013).
Plant diversity is expected to play a key role in influencing plant-pollinator interactions as it is related to the availability and heterogeneity of floral resources. For instance, pollinator specialization has been showed to increase with increasing richness of flowering plant species (Ebeling, Klein, & Tscharntke, 2011).
Here, we explored how a highly invasive plant species (Buddleja davidii Franch.) impacted plant-flower visitor interactions over its whole elevational range distribution, and whether this effect was modified by managed honeybees and available alternative flower resources. We compared invaded and equivalent non-invaded riparian habitats along a gradient in honeybee abundance and species richness of flower resources by sampling Hymenoptera, Diptera and Lepidoptera pollinators. First, we tested for a potential interaction between B. davidii and honeybee in shaping the stability of plantflower visitor networks. We expected that, in the presence of the highly rewarding invasive plant species, honeybees would focus on B. davidii, thus decreasing network stability. In non-invaded sites instead, honeybees would forage on many plant species, consequently resulting in a higher network stability via increased network connectance. Second, we tested how B. davidii, honeybee abundance and species richness of flowering plants influenced resource overlap of wild flower visitors with honeybees and their species-level specialization, and whether these effects changed among insect orders.
We expected that, in sites invaded by the exotic plant species, a large proportion of flower visitors would concentrate their foraging on the novel resource, consequently increasing their specialization.
However, due to their aggressive foraging behaviour, honeybees are expected to influence how other species visit flowers (Hung et al., 2019;Santos et al., 2012), usually diverging wild flower visitors on less abundant resources (Goulson, 2003;Hung et al., 2019;Thomson, 2016) and thus increasing their apparent specialization.

K E Y W O R D S
Alien species, Apis mellifera, bipartite networks, competition, elevational gradient, network stability, specialization 2 | ME THODS

| Study system
Fieldwork was carried out in the mountain areas of Veneto and Trentino-Alto Adige regions (Northern Italy; Figure S1). Climate varies with elevation: it is continental in the lowlands, with warm summers and mild winters, whereas higher elevations are characterized by alpine climate, with cool summers, cold winters and frequent snowfalls. Precipitation is abundant year round, especially at higher elevations (c. 1,200 mm). The area is invaded by the perennial shrub B. davidii, which can be found up to 1,200 m a.s.l., especially in disturbed areas such as railways and river banks. It is native of China and Japan and was introduced in Europe in the late 1800s. Selfincompatibility is high (>95%), and the species requires insect crosspollination for reproduction (Ebeling, Schreiter, Hensen, Durka, & Auge, 2012). The species is highly attractive to pollinators due to its conspicuous inflorescences rich in nectar. Because of their floral morphology (narrow and long corolla tube), species of the genus Buddleja are mostly pollinated by butterflies, honeybees and other large bees (Gong et al., 2015). The species exhibits a long flowering period that usually starts in late spring and lasts until late summer.

| Sampling design
We selected nine pairs of 40 × 20 m sites across the whole elevational range distribution of B. davidii (c. 100-1,200 m a.s.l.; Table S1).
All the sampled sites were open riparian habitats located along valley bottoms. For each pair, we chose one site invaded and one noninvaded by B. davidii. Within each pair, the surrounding landscape of each site was similar in terms of habitat composition (Table S2).
In invaded sites, the abundance of B. davidii ranged from 2.3% to 38% of all flowering species (mean = 15.9, SD = 9.7%), whereas it was completely absent in control sites. We tested for differences in

| Surveys of flowering plants and flower visitors
During spring and summer 2018, we observed plant-flower visitor networks in the 18 study sites. Sites were visited approximately every three weeks, for a total of five surveys covering the full flowering season of B. davidii. The first survey was conducted in June, while the last survey occurred at the end of August, as the flowering started to cease. Plants were identified in the field when possible, or collected and prepared in a herbarium and identified later.
Insect sampling occurred between 09.00 hr and 17.00 hr in dry and sunny conditions with low wind and temperature above 18°C. At each round, each flowering species was observed for a 5-min period, during which all insects touching the reproductive parts of the flowers were counted or collected using a butterfly net (Ø 35 cm). Every 15 min, we paused for 10 min to reduce the level of disturbance and let the flower visitor community recover. Our survey included the most abundant and diverse groups of flower visitors, that is all the bees and sphecids (Hymenoptera), hoverflies, conopids and tachinid flies (Diptera) and butterflies (Lepidoptera). Coleopterans, less than 1% of all observed insects, were not collected and identified. Insects were identified in the field when possible or placed in vials filled with 70% ethanol and sorted and identified in the laboratory to species or morphospecies level. As females of Bombus terrestris (Linnaeus) and Bombus lucorum (Linnaeus) are often difficult to distinguish, we pooled both species (males and females) as B. terrestris/lucorum.

| Network stability
The network stability was quantified using the robustness index (Burgos et al., 2007;Memmott, Waser, & Price, 2004). Robustness was calculated by removing flower visitor species from the network according to their abundance, from the rarest to the most abundant. Robustness ranges from 0 (highly unstable network) to 1 (highly stable network). It is well known that number of species (i.e. network size) and sampling effort can affect the number of observed interactions in real ecological networks, and thus several measures of network structure (Delmas et al., 2019). This issue can seriously limit the interpretation of raw network measures and their use for multiple network comparison (Chagnon, 2015). In order to account for network size effect when comparing different networks, the robustness index was standardized using null models as ΔRobustness = observed robustness-robustness null , where robustness null represents the mean robustness value from 1,000 randomized networks obtained using the Patefield, the Vazquez and the swap.web algorithms (Dormann, Fründ, & Gruber, 2008;Patefield, 2012;Schleuning et al., 2012;Vázquez et al., 2007).
The three algorithms respectively constraint marginal totals of rows and columns, connectance and both marginal totals and connectance. As standardized robustness values obtained with the different algorithms were highly correlated with each other and with non-standardized ones ( Figure S3), we decided to present and discuss only the values obtained with the Vazquez algorithm.
Both raw and standardized robustness values were calculated using the bipartite package (Dormann et al., 2008).

| Resource overlap and flower visitor species specialization
For each wild flower visitor species within each network, the resource overlap with honeybees was calculated using the Morisita's index (Morisita, 1959) in the spaa package (Zhang, 2016). The metric ranges from 0 (no overlap) to 1 (complete overlap).
As species-level specialization can be defined in many different ways, we calculated three metrics that convey complementary information (Maglianesi, Blüthgen, Böhning-Gaese, & Schleuning, 2014): (1) normalized degree, (2) dʹ and (3) PDI (Paired Difference Index). The normalized degree is one of the simplest specialization indices, based on the presence/absence of flower visitors on plants, and calculated as the standardized ratio between the number of plant species on which a flower visitor species was observed and the total richness of flowering plant species in the network. It ranges from 0 (highly specialized) to 1 (highly generalized). The dʹ index (standardized Kullback-Leibler distance) additionally takes into account the abundance of flower visitors observed on a plant species, and how much the resource is shared by different flower visitor species (Benadi et al., 2014;Blüthgen et al., 2006). In PDI, the strongest link between a flower visitor species and its most visited plant is contrasted with those over all remaining plant species within the network (Poisot, Canard, Mouquet, & Hochberg, 2012).
Both dʹ and PDI range from 0 (no selectiveness) to 1 (extreme selectiveness). Normalized degree was transformed as the inverse of the index (1-normalized degree) in order to make it comparable with the other two specialization indices. Defining specialization for rare species may be critical because of the risk of overestimating their specialization (Dorado, Vázquez, Stevani, & Chacoff, 2011;Dormann, 2011). To calculate the species-level specialization indices, we thus removed all singletons and doubletons from the dataset before using the specieslevel function in the bipartite package (Dormann et al., 2008). Moreover, we computed the three species-level specialization indices excluding the honeybee from the network. The three specialization metrics were strongly correlated with those calculated with the honeybee in the network, so we present and discuss only the values calculated including the honeybee within the network. However, we run the analyses with both methods and found similar results (Table S3).

| Effect of B. davidii and honeybee abundance on network stability
A linear mixed-effects model, with the site pair as random factor to account for spatial dependence in the design, was used to test how network stability (i.e., robustness index) responded to the presence of B. davidii (invaded vs. control sites), honeybee abundance and flowering plant species richness. We also included in the model the interaction between B. davidii and honeybee abundance. Both abundance of honeybees and species richness of flowering plants were ln-transformed to improve linearity and model residual distribution.
The model was fitted using the function lme in the nlme package (Pinheiro, Bates, DebRoy, Sarkar, & Team, 2017). Interactions were removed using a backward stepwise model selection procedure (pvalue > .05). Model assumptions were checked using residual diagnostic plots. Effects were represented using partial residual plots using the visreg package (Breheny & Burchett, 2017). wild Hymenoptera, Diptera and Lepidoptera). In the case of resource overlap, one pair was removed from the dataset due to honeybee absence in one of the site. We included sampling site nested within the pair as random factor to account for spatial dependence in the design. Resource overlap, normalized degree, honeybee abundance and flowering plant species richness were ln-transformed in order to improve linearity and model residual distribution. Non-significant interactions were removed using a backward stepwise model selection procedure (p-value > .05). Model assumptions were checked using residual diagnostic plots. All models were fitted using the function lme in the nlme package (Pinheiro et al., 2017). Effects were represented using partial residual plots using the visreg package (Breheny & Burchett, 2017).

| Potential collinearity between predictors
The selection of the sites guaranteed statistical independence between honeybee abundance and elevation (Pearson's correlation = −0.207, p-value = .409) and between honeybee abundance and flowering plant species richness (Pearson's correlation = 0.368, p-value = .133), while the species richness of flowering plants significantly increased with elevation (Pearson's correlation = 0.549, p-value = .018). We also measured multicollinearity in our models by computing the variance inflation factors (VIF). The VIF values were all below 1.5, indicating no multicollinearity in our models.

| Elevation versus number of flowering plants
As elevation was correlated with flowering plant species richness and exhibited a narrow range of variation, we only used species richness of flowering plants as explanatory variable in all the models described above. However, in preliminary analyses, we ran all models using elevation instead of number of flowering plants as explanatory variable, and found no effect of elevation on network stability, resource overlap and species specialization. This result may be explained by the relatively low upper limit of B. davidii distribution (c. 1,200 m a.s.l.) that implied a narrow temperature gradient between our low and high elevation sites. In temperate mountains such as the European Alps, pollinator communities do not usually exhibit strong species turnover between lowlands and mid-altitudes (c. 1,500 m a.s.l.; Lefebvre, Villemant, Fontaine, & Daugeron, 2018), meaning that foraging preferences and species composition are quite stable along our elevational gradient.

| General results
Overall, we recorded 9,563 interactions between 370 flower visitor species and 150 flowering plant species (25 plant species were not visited by any insect; Tables S4 and S5), for a total of 1,737 unique plant-flower visitor interactions. Beside B. davidii, we found 10 more exotic plant species (Table S5). As often found in real mutualistic networks (Devoto, Bailey, Craze, & Memmott, 2012), we only observed a subset of the expected interactions based on rarefaction curves, but with no evident systematic differences in sampling completeness between sites invaded by B. davidii and control sites ( Figure S4).

F I G U R E 1 Partial residual plots
showing the effects of (a) honeybee abundance (ln-transformed) and (b) species richness of flowering plants (lntransformed) on standardized network robustness. The robustness index was standardized as ΔRobustness = observed robustness-robustness null , where robustness null represents the mean robustness value from 1,000 randomized networks obtained using the Vazquez algorithm (Vázquez et al., 2007)  Network robustness standardized using the Vazquez algorithm was positively influenced by both the abundance of honeybees ( Figure 1a) and the species richness of flowering plants (Figure 1b), while no effect of B. davidii presence was observed ( Table 1). The effects were consistent between models using robustness standardized with the Patefield and the swap.web algorithms, while effects were slightly different for the model using the raw robustness values (Table S8) with the effect changing among insect orders (Table 3); that is, dʹ increased with increasing honeybee abundance for dipterans and lepidopterans, whereas for hymenopterans, it was high irrespective of honeybee abundance (Figure 4). The PDI specialization index was only positively influenced by the number of flowering plant species (Table 3; Figure 5).

| D ISCUSS I ON
The effect of invasive plants on pollinators has been intensively studied, but often in isolation from other biotic and abiotic drivers. Here, for the first time, we tested for the interactive effect of a highly invasive plant species (B. davidii) and a competitive managed pollinator species (the honeybee) on plant-flower visitor interactions. We found no interaction between B. davidii and honeybee abundance on both network-and species-level metrics. However, the two explanatory variables showed contrasting effects on network stability and flower visitor specialization. While the presence of B. davidii decreased flower visitor specialization, honeybee abundance generally increased both network robustness against flower visitor extinctions and species specialization, but the latter with some differences among insect orders.

| Network stability increased with honeybee abundance and richness of flowering plant species
Super-generalist species such as the honeybee are able to exploit many plant species and thus to potentially increase network robustness against flower visitor extinctions. Although a previous study found that the removal of honeybees from networks did not affect overall stability (Santos et al., 2012), here on the contrary we showed that increasing honeybee abundance increased network stability. Highly stable networks, that is those with a high honeybee abundance, are more robust to the loss of rare flower visitor species, as only a limited number of plant species will not be visited anymore when rare species become extinct. However, as honeybees are sometimes less efficient pollinator than wild species, this increased robustness might not necessarily translate into higher reproductive success of those plants that are visited by honeybees (Magrach et al., 2017;Valido et al., 2019). The positive effect of species richness of flowering plants on network stability was probably related to the

TA B L E 2
Results from the linear mixed-effects model testing the effects of explanatory variables on resource overlap of wild flower visitors with honeybees, and how it changes among insect orders (i.e., Hymenoptera, Diptera and Lepidoptera) positive correlation between the number of flowering species and both pollinator richness and abundance; that is, sites rich in plant species hosted more species and individuals that in turn promoted network robustness (Thébault & Fontaine, 2010).

| Resource overlap between wild flower visitors and honeybees depended on insect order
The interactions between honeybees and wild flower visitors may be explored using the resource overlap index, where high levels of overlap indicate potential competition between wild and managed species (Franco, Aguiar, Ferreira, & De Oliveira-Rebouças, 2009;Paini, 2004). Resource overlap between honeybees and other pollinators has been recorded before (Goulson, 2003), and an increase in honeybee abundance is expected to intensify competition for floral resources between wild and managed bees (Thomson, 2016).
Here, we found that resource overlap was higher for lepidopterans than for hymenopterans and dipterans. Despite B. davidii being highly attractive for both butterflies and honeybees, the observed effect did not change in invaded vs. uninvaded sites, meaning that the high level of overlap between lepidopterans and honeybees is probably due to other reasons. One hypothesis might be related to the morphology and length of the mouthparts in the three insect orders: while lepidopterans have medium-long proboscises that allow to efficiently exploit flowers with different shapes and corolla depth, most hymenopteran (except for bumblebees) and dipteran species are short-tongued, and usually prefer open flowers (Kevan & Baker, 1983). Therefore, the honeybee, as a medium-tongued bee (c. 6 mm), mostly overlaps with lepidopterans. Consistent with previous studies (Gillespie & Elle, 2018;Steffan-Dewenter & Tscharntke, 2000), we found that honeybees did not appear to share flower resources with both hymenopterans and dipterans.

| Flower visitor specialization decreased with B. davidii presence and increased with increasing honeybee abundance
As invasive plants are often rich in pollen and nectar, they can potentially compete with native plants for pollinators (Dietzsch, Stanley, & Stout, 2011;Flanagan, Mitchell, & Karron, 2010), causing a reduction in the range of consumed resources and thus potentially resulting in increased specialization of pollinator species (Maruyama et al., 2016).
Conversely, we found that, in invaded sites, flower visitors used a broader range of plants, reflected by a decrease in their specialization. A possible explanation is that the invasive plant attracts flower visitors, which in turn visit also co-occurring species in the neighbourhood (Molina-Montenegro, Badano, & Cavieres, 2008). For example, some bumblebee species perceive adjacent plants (closer than 6 m) as a single resource, even when they belong to different species (Klinkhamer, de Jong, & Linnebank, 2001). Additionally, when foraging on the preferred species, other flowering species are regularly visited in order to compare their resource levels (Albrecht et al., 2016). B. davidii, thanks to its scented and showy inflorescences, is highly attractive and might act as a facilitator for co-occurring native plants. This "magnet species" effect was also recorded in areas invaded by other exotic plant species (Albrecht et al., 2016;Molina-Montenegro et al., 2008).
The attractive effect of exotic plants was previously found to be particularly evident for highly generalist pollinator species, such as the honeybee (Hung et al., 2019;Magrach et al., 2017). Because of its foraging behaviour, this species may force wild pollinators to shift on less abundant resources (Goulson, 2003;Hung et al., 2019;Thomson, 2016). Depending on how flexible other species are in modifying their diet breadth, the effect of honeybee on wild pollinators may be different. Here, we found that an increase in honeybee abundance did not modify the number of visited plant species or how flower visitors were partitioned among flower resources, but affected how the latter were shared by flower visitor species (dʹ). As honeybee abundance increased, lepidopterans and dipterans tended to switch to less visited plant species, possibly in order to avoid competition (Hung et al., 2019). On the other hand, the effect of honeybee abundance on wild hymenopterans was less evident. It did not lead to a switch on less visited species, but it decreased the evenness of wild hymenopteran visits on floral resources, suggesting that some resources became more preferred. Despite the common knowledge of lepidopterans having intermediate level of specialization compared with hymenopterans (more specialized) and dipterans (less specialized) (Benadi et al., 2014;Chacoff et al., 2012), in our study lepidopterans displayed the highest degree of specialization in term number of visited plant species. This is probably due the high attractiveness of some plant species to butterflies, including B. davidii (Gong et al., 2015). Indeed, in our study over the 76% of lepidopterans was collected in sites invaded by the exotic species, and the 36% was directly foraging on B. davidii. Among the remaining 24% lepidopterans collected in uninvaded sites, over the 33% was collected on Campanula persicifolia L. in a single sampling site. The observed high selectivity of lepidopterans is thus likely to depend on the identity of the plant species rather than the ecology of the order.

F I G U R E 2
In contrast to previous studies (Ebeling et al., 2011;Fründ, Linsenmair, & Blüthgen, 2010), we found that the number of flower-

| CON CLUS IONS
Contrary to our expectations, we found no interactive effect of B.
davidii and honeybee abundance on plant-flower visitor interactions.
However, we showed that the invasive plant and managed honeybees can jointly affect flower visitor communities. Although the honeybee is sometimes considered to be a less efficient pollinator

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13132.

DATA AVA I L A B I L I T Y S TAT E M E N T
Primary data are provided in the main text or in the Appendix S1.