DOI: 10.6086/d1rd69
title: Wild bee functional diversity and plant associations in native and conventional plant nurseries

description: An ongoing challenge in ecology is predicting how characteristics of
 communities correspond to habitat features. Examining variation in
 functional traits across species may reveal patterns not discernible from
 measures of mere abundance or richness. For beneficial insects like wild
 bees, functional trait-based approaches are often used to characterize
 communities in different agricultural habitats. However, no such approach
 has yet been applied in horticultural plant nurseries, which represent
 intensively managed artificial flowering plant assemblages. Certain
 nurseries mostly cultivate regionally native flowering plants, allowing
 one to test how differences between local plant assemblages may correlate
 with bee functional traits. We surveyed bee assemblages at native and
 conventional plant nurseries in southern California from spring through
 autumn over two years, while also documenting the native status of
 blooming plants in sampling plots. Bees were classified into different
 functional categories based on their diet breadth, nesting location, and
 social organization. At native plant nurseries, we netted proportionally
 more specialist bee species and fewer generalist species than at
 conventional nurseries. Nesting location and social organization of bee
 samples did not differ between nursery types. Meanwhile, landscape-level
 features were not associated with any observed functional trait metrics of
 bee collections. Furthermore, network-level specialization of bee-plant
 interactions was higher at conventional nurseries. Our results suggest
 that a horticultural cultivation practice is quantifiably correlated with
 the functional distribution of resident bee assemblages. These results are
 important and encouraging to pollinator conservation efforts in nursery
 systems and other human-modified landscapes dominated by ornamental
 plants.
description: <i>Study sites</i> We
 conducted this study at 13 plant nurseries in southern California, USA
 (Table S1). All nurseries were open-air yards with containerized flowering
 plants grown outside or under shade houses, thus being freely accessible
 to foraging bees (Fig. S1). From March 2016 to May 2018, we sampled each
 nursery every three months, except in winter: once each in spring
 (March–May), summer (June–August), and autumn (September–November). Due to
 the number of nurseries and distance between them, it was not feasible to
 sample all nurseries simultaneously. Within each season, we sampled
 nurseries in a semi-randomized order within a period of about 6.5 weeks.
 We sampled each nursery at least twice, with the exception of one,
 resulting in 58 sampling events (Table S2). We classified a nursery as a
 “native” nursery if it met two criteria: (1) it advertised itself as
 cultivating, mostly or exclusively, plants native to California, and (2)
 on average, greater than 50% of blooming plant taxa documented in surveys
 (see below) across sampling events were native to California. All five
 nurseries that met the second criterion also met the first, with mean (±
 SEM) percent native richness at these nurseries ranging from 72 ± 4% to 91
 ± 2%. Eight nurseries not meeting both criteria were classified as
 “conventional”, with mean percent native richness ranging from 1 ± 0.7% to
 45 ± 5%. <i>Bee sampling
 in nurseries</i> In each
 sampling event, we used both blue vane traps (SpringStar, Woodinville, WA)
 and sweep nets to collect bees. Vane traps detect a wide diversity of bees
 (Rhoades <i>et al.</i>, 2017) even near attractive floral
 resources. Moreover, conventional bowl traps for bees are vulnerable to
 flooding from overhead irrigation and disturbance from passersby. We
 deployed vane traps at a density of roughly one trap per hectare of
 nursery property (2–12 traps per nursery). Each vane trap was placed in
 the open, partially filled with soapy water, and suspended from a 1.2-m
 hook, then collected after 72 h. In every sampling event, we placed vane
 traps in roughly the same location at each nursery. Each vane trap
 represented the center of a circular sampling plot used for collecting
 data on blooming plants (see below).
 Each time we visited a nursery to deploy or collect
 vane traps (i.e., twice per sampling event), one observer sweep netted
 wild bees from flowering nursery stock for 30 minutes. Daily high
 temperatures on sweep netting days were at least 18 °C, wind velocities
 were low, and there was no precipitation. We recorded the associated plant
 species for each collected bee, and spent no more than five minutes
 (including specimen handling time) collecting from any one plant species.
 Non-native honey bees, <i>Apis mellifera</i>, were not sweep
 netted due to their high prevalence on flowering plants at
 nurseries. We identified all
 bees to species (or morphotaxon, if species-level identification was
 unfeasible) using published keys (Michener, 2007; Ascher &
 Pickering, 2017), reference collections, and consultation with bee
 specialists at the University of California, Riverside Entomology Research
 Museum. Using published literature, we categorized each species on the
 basis of three functional traits: pollen diet breadth, nesting location,
 and social organization (Table S3). If these data were not available for a
 taxon, we assigned traits based on those of its congeners (as in Wilson
 & Jamieson, 2019). For diet breadth, species were classified as
 generalists (collecting pollen of plants from many unrelated families to
 provision offspring, i.e., polylectic), specialists (oligolectic,
 mesolectic, or monolectic), or cleptoparasites that do not collect pollen
 (Cane & Sipes, 2006). For nesting location, species were
 categorized as above-ground nesters (including cavity-renting and
 -excavating species nesting in stems or wood), below-ground nesters (soil
 excavators and burrow renters), or cleptoparasites that reproduce using
 host nests. For social organization, species were classified as eusocial
 (including primitively eusocial), solitary (including subsocial and
 communal), or cleptoparasites reliant on a host for reproduction (Buchholz
 & Egerer, 2020).
 <i>Local and landscape features of
 nurseries</i> In each
 sampling event, we surveyed blooming plants inside a 15-m radius (700
 m<sup>2</sup>) circular plot around each vane trap. We
 recorded all plant species, including weeds, currently blooming in each
 plot (as in Wilson & Jamieson, 2019) regardless of floral
 abundance. We did not record graminoids, angiosperms not currently
 blooming, or non-angiosperms. Using the online CalFlora database
 (calflora.org), we classified each species as native if any part of its
 native range fell within California state boundaries (as in Avolio
 <i>et al.</i>, 2020). Hybrid cultivars were classified as
 native if all parent taxa were native. We also estimated percent floral
 cover within these same plots. We divided plots into quadrants, then
 assigned an ordinal value from 0 to 4 to each quadrant: ‘0’ = no flowers
 present, ‘1’ = few flowers present, ‘2’ = more than ‘1’ but covering less
 than 50% the area of the quadrant, 3 = flowers covering &gt; 50%
 quadrant, but with large patches lacking flowers, and ‘4’ = flowers
 occupying near 100% quadrant area. Values were averaged across quadrants
 within plots. Using QGIS (QGIS
 Development Team, 2020), we calculated nursery area and the proportion of
 natural area within buffers of varying radii (500 m and 1 km) around each
 nursery’s perimeter. These distances encompass the flight ranges of many
 wild bee species (Greenleaf <i>et al.</i>, 2007). We used the
 2016 National Land Cover Database (mrlc.gov, 30-m resolution) to classify
 land into cover types. We defined “natural” area as cover classes
 including forests, shrubland, grassland, and wetlands (classes 41–43, 52,
 71, 90, 95) and excluding water or land that is developed, managed for
 agriculture, or barren (11, 21–24, 31, 81, 92).
 <i>Statistical analyses: Bee functional
 traits</i> We conducted all analyses in R version
 3.3.3 (R Core Team, 2021). All means are reported ± SEM.<i>
 </i>In 7 of 58 sampling events, no bees were detected while sweep
 netting, so these events were excluded from sweep netting analyses
 (<i>N</i> = 51). We pooled captures from all vane traps per
 sampling event per nursery so as to render analyses comparable between
 sweep netting and vane trapping, as sweep netting was not conducted within
 plots. As honey bees were not sweep netted, we excluded from analyses the
 633 specimens incidentally collected in vane traps. We re-analyzed the
 vane trap dataset <i>including</i> honey bees to detect
 whether their removal affected model results (Table S4). For each
 collection method in each sampling event, we tallied the total number of
 bee species collected (richness) and the number of species in each
 category of the functional trait of interest (diet breadth, nesting
 location, or social organization). We converted these count data to
 proportions to account for differences in the numbers of species collected
 across samplings. To examine variation in
 the proportions of species across functional trait categories, we
 constructed generalized linear mixed models in <i>lme4</i>
 (Bates <i>et al.</i>, 2015) with binomial error distributions.
 We created six models, one for each functional trait and collection
 method. In models, the dependent variable was the proportion of bee
 species in each category of that functional trait; thus, each sampling
 event corresponded to three observations in a model. We weighted each set
 of three proportions by the total number of species collected in that
 sampling event. We included as fixed effects the corresponding functional
 trait (factor, <i>k</i> = 3), season, and nursery type, and
 all interactions between these factors, and nursery area and proportional
 surrounding natural area (1 km). We included nursery and year as random
 factors. Similarly, we constructed six extra models using proportions of
 <i>individuals</i> collected in sampling events as dependent
 variables (Table S4). We excluded average blooming plant richness and
 floral cover per sampling event from all mixed models, as their inclusion
 resulted in significantly inferior model fits (Table S5) and these
 variables had no effects (all P &gt; 0.10) on the proportion of
 species collected in any functional trait category or on functional
 dispersion of bee samples (Table S5).
 We calculated functional dispersion (FDis) of bees
 collected in each sampling event using function ‘dbFD’ with a Cailliez
 correction (Forrest <i>et al.</i>, 2015) in
 <i>FD</i> (Laliberte & Legendre, 2010). FDis
 quantifies the functional trait diversity of a sample, is weighted by
 abundances, and is not influenced by richness. We constructed two linear
 mixed models, one for each collection method, with FDis as the dependent
 variable and included the same fixed and random effects as
 above. We retained all
 non-significant (<i>P</i> &gt; 0.05) interaction terms in
 models to facilitate comparisons across models. We confirmed the absence
 of multicollinearity (VIF &lt; 3) among fixed effects using function
 ‘vif’ (<i>car</i>) (Fox & Weisberg, 2011). Lastly, we
 conducted post-hoc Tukey’s HSD tests using function ‘emmeans’
 (<i>emmeans</i>)<i> </i>(Lenth,
 2019)<i>.</i> Additionally, to determine if spatial scale
 influenced our model results, we also compared our original models to a
 different set of models which quantified the predictor variable of
 proportional natural area at a 500-m buffer (Table S5). Models using 500-m
 buffers fit the data as well as those using 1-km buffers (ΔAICc &lt;
 2.0; Table S5), and the proportions of natural area at the two buffer
 radii were highly correlated across the 13 nurseries (Pearson
 <i>r</i> = 0.956). For simplicity, we report only models using
 the 1-km buffer data.
 <i>Statistical analyses: Bee-plant
 networks</i> Using our
 sweep netting dataset, we subdivided our overall bee-plant interaction
 network matrix into sub-matrices representing each sampling event. We
 calculated the size of each sub-network as the product of the number of
 bee species collected and the number of plant species off which bees were
 netted (de Mendonça Santos <i>et al.</i>, 2010). We examined
 variation in sub-network size by constructing a linear mixed model with
 the aforementioned fixed and random effects.
 To examine generalization in bee-plant interactions
 in nurseries, we employed a null network approach using package
 <i>econullnetr</i> (Vaughan <i>et al.</i>, 2017).
 We compared our observed bipartite bee-plant networks to null networks
 created using function ‘generate_null_net’. Null networks were structured
 such that bees interact with plants in proportion to their relative
 abundances in observed networks, using seasons as subdivisions. Using
 function ‘bipartite_stats’, we compared observed values of network metrics
 to 95% confidence intervals from 1,000 iterations of null models;
 significantly different values fall outside this interval. Standard effect
 sizes for metrics were generated for comparisons between networks. To
 detect differences in generalization between networks from native and
 conventional nurseries, we divided our overall network into two
 sub-networks corresponding to nursery type, and compared their nestedness
 (NODF), specialization (H2), and linkages (connectance) to their
 respective null networks. To test for differences in generalization
 between native and non-native plant species and among bee diet breadth
 categories across nurseries, we used these groupings as factors in
 separate linear models to compare standardized effect sizes for centrality
 (closeness), specialization (d’), and linkages (effective
 partners). Finally, to assess
 if any bee species in our sweep netting dataset were significantly
 associated with either native or non-native ornamental plants, we
 conducted a multi-level pattern indicator species analysis using function
 ‘multipatt’ in package <i>indicspecies </i>(De Cáceres
 & Legendre, 2009) with 10,000 iterations. This analysis identifies
 species whose occurrence is associated with groups of sites organized into
 categorical habitat types. Here, we treated plant species analogously to
 sites, categorizing them based on native status.
description: Each Excel file contains a metadata tab (first tab) which
 explains the nature of each sheet in the file as well as each variable
 (column) in the file.
 "plant.nativity.richness.xlsx": richness,
 including proportional native richness, of plants at nurseries.
 "functional.traits.sweep.netting.xlsx": functional
 trait distributions of wild bees collected via sweep netting.
 "functional.traits.vane.trapping.xlsx": as above, but
 for vane traps.
 "bipartite.network.analyses.xlsx": bee-plant
 interaction network data.