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.

Study sites

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%.

Bee sampling in nurseries

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 et al., 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, Apis mellifera, 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).

Local and landscape features of nurseries

In each sampling event, we surveyed blooming plants inside a 15-m radius (700 m²) 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 et al., 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 > 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 et al., 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).

Statistical analyses: Bee functional traits

We conducted all analyses in R version 3.3.3 (R Core Team, 2021). All means are reported ± SEM. In 7 of 58 sampling events, no bees were detected while sweep netting, so these events were excluded from sweep netting analyses (N = 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 including 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 lme4 (Bates et al., 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, k = 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 individuals 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 > 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 et al., 2015) in FD (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 (P > 0.05) interaction terms in models to facilitate comparisons across models. We confirmed the absence of multicollinearity (VIF < 3) among fixed effects using function ‘vif’ (car) (Fox & Weisberg, 2011). Lastly, we conducted post-hoc Tukey’s HSD tests using function ‘emmeans’ (emmeans) (Lenth, 2019). 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 < 2.0; Table S5), and the proportions of natural area at the two buffer radii were highly correlated across the 13 nurseries (Pearson r = 0.956). For simplicity, we report only models using the 1-km buffer data.

Statistical analyses: Bee-plant networks

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 et al., 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 econullnetr (Vaughan et al., 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 indicspecies (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.

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.