In human-modified landscapes, such as agricultural areas, understanding how habitat characteristics influence the diversity and composition of beneficial organisms is critical both to conservation efforts and to modeling ecosystem services. Assessing turnover, or changes in species identity across sites or through time, is crucial to determine how shifts in community composition relate to changes in ecosystem services. For pollinators like wild bees, factors affecting their diversity have been well studied, but variables influencing temporal turnover, particularly across seasons within a year, remain relatively poorly understood. To investigate how local and landscape characteristics correlate with bee diversity and turnover across seasons, we recorded wild bee and flowering ornamental plant assemblages at 13 plant nurseries in southern California between spring and autumn over two years. Nurseries cultivate a broad diversity of flowering plant species that differ widely across sites and seasons, providing a field setting to test for correlations between turnover and diversity of plants and bees. As expected, we documented strong seasonal trends in wild bee diversity and composition. We found that local habitat factors, such as increased cultivation of native plants, were positively associated with bee diversity in sweep netting collections, whereas we detected moderate influences of landscape-level factors like proportion of surrounding natural area in passive trap collections. We also detected a moderate positive correlation between overall gains in plant species and gains in bee species at nurseries across consecutive seasons. Our results have important considerations for the conservation of wild bees in landscapes dominated by ornamental plants, and highlight the utility of plant nurseries for investigating hypotheses related to diversity and turnover in plant-pollinator systems.

Study sites and bee collection

We collected wild (non-Apis) bees at 13 wholesale or retail plant nurseries (each > 2 ha in area) in southern California, USA (Table S1, Fig. S1). We sampled each nursery from March 2016 – May 2018 once per season, except winter: once each in spring (March–May), summer (June–August), and autumn (September–November) (Table S2). We used two methods to collect bees. First, we deployed blue vane traps (SpringStar, Woodinville, WA) suspended from 1.2-m hooks in established sampling plots (see below) at a density of roughly one plot per hectare of nursery property (2-12 plots per nursery). We collected vane traps 72 h after deployment. Second, a single observer used a sweep net to collect bees visiting flowering plants in nurseries for 30-minute periods on each day we deployed and collected vane traps, spending no longer than 5 minutes collecting from a given plant species. Non-native honey bees, Apis mellifera, were abundant on many plants at these nurseries and were intentionally not sweep netted, but were sometimes collected in the vane traps. We identified all wild bees to species (or morphotaxon, if they could not be confidently assigned to species) using keys (Michener 2007, Ascher and Pickering 2020), reference collections, and assistance from bee specialists at the University of California Riverside Entomology Research Museum.

Local and landscape features

We surveyed flowering plants that were currently blooming, regardless of floral abundance, within a 15-m radius (~700 m²) plot around each vane trap during each nursery visit. We designated each plant species as native or non-native to California referencing the CalFlora database (calflora.org). We further classified five nurseries as “native nurseries” because they (1) advertised themselves as cultivating mostly plants native to California, and (2) on average, over two-thirds of plant taxa documented in plots were California-native. We classified the eight other nurseries as “conventional nurseries”. We estimated floral cover as percent cover within plots by dividing the plot into quadrants and assigning a rank to each quadrant using an ordinal scale from 0 to 4: ‘0’ = no flowers of any species, ‘1’ = few flowers with sparse cover, ‘2’ = more than ‘1’ but covering less than 50% of quadrant, 3 = > 50% cover but including large patches with no flowers, and ‘4’ = near 100% cover, then averaging these values for the plot. Using QGIS (QGIS Development Team 2020) and the 2016 National Land Cover Database (mrlc.gov), we calculated the proportion of natural area inside 500-m and 1-km buffers around each nursery’s perimeter. We treated land cover classes 41–43, 52, 71, 90, and 95 as “natural” landscape, excluding classes denoting water and developed, agricultural, or barren land (11, 21–24, 31, 81, and 92). We also recorded from local weather stations the daily high temperature for each of the four calendar days that vane traps were deployed, which we then averaged.

Statistical analyses

We conducted all statistical analyses in R version 3.3.3 (R Core Team 2021). All means are reported ± SEM. The dataset represented 58 sampling events for blue vane traps (of which 21 were at native nurseries) and 80 sampling events for sweep netting (32 at native nurseries). For seven of these sweep netting sampling events, no wild bees were detected. We excluded all Apis mellifera (633 specimens) from the vane trap dataset so as to focus on patterns in wild bee assemblages. Comparisons of analyses with the full dataset and the dataset excluding honey bees revealed that their removal had minimal to no influence on the outcomes of our statistical models (Table S3).

We performed individual-based rarefaction using vegan (Oksanen et al. 2018) and iNEXT (Hsieh et al. 2016) to examine sampling completeness of wild bees across sites, nursery types, and collection methods. For each sweep netting and vane trap sampling event, we calculated bee Shannon diversity and evenness. For each survey plot, we calculated the richness of blooming plants and the proportion of richness represented by native plant species, then averaged across plots within a given sampling event. In statistical models, we substituted proportional native plant richness with the residuals of this variable obtained from a regression against nursery type, as these predictors were collinear (variance inflation factor > 3) (Graham 2003). The same was also done with average daily high temperature in blue vane trap models, which was collinear with season. We constructed four separate linear mixed models (LMMs) using lme4 (Bates et al. 2015) for each of bee Shannon diversity and evenness and each of two collection methods. We included as fixed effects season and nursery type (and their interaction), blooming plant richness, proportional native plant richness, percent floral cover, daily high temperature, nursery area, and proportional surrounding natural area at 1 km. Models in which natural area was quantified at 500 m were not significantly better than those at 1 km (all ΔAIC < 2.0), thus we only report models using the 1-km buffer data. Nursery, study year, and date of sampling served as random effects in LMMs. For LMMs, we confirmed lack of multicollinearity (VIF < 3) using ‘vif’ (car) (Fox and Weisberg 2011). Using the function ‘emmeans’ in the package emmeans (Lenth 2019), we conducted post-hoc Tukey’s tests to compare treatment means for statistically significant factors and corrected for multiple comparisons.

We next examined variation in the composition of plant and wild bee assemblages in nurseries. For bees, we conducted these analyses only on vane trap samples, as sweep net samples were not large enough to sufficiently represent assemblage composition. We tested for differences in the multivariate centroids of samples using permutational multivariate analyses of variance (permANOVA) on dissimilarity matrices using ‘adonis’ (vegan). For bees, we used Bray-Curtis dissimilarity as it accounts for the relative abundances of species, whereas for plants, we used Jaccard dissimilarity as plant species were recorded in plots as present or absent. We included the aforementioned fixed effects and stratified permutations (10,000) by nursery. We also examined beta diversity by analyzing multivariate homogeneity in dispersion of plant and wild bee assemblages around their centroids using ‘betadisper’ (vegan) between nursery types, seasons, and years. We performed Mantel tests using ‘mantel’ (vegan) to explore correlations between dissimilarity matrices of plant and bee (vane trap samples) assemblages across nurseries separately in spring, summer, and autumn.

We used function ‘RAC_change’ in the package codyn (Hallett et al. 2016) to calculate relative changes in plant and bee richness through time (Kazenel et al. 2020), pooling bees from both collection methods. Temporal species turnover can be partitioned into gains and losses (Faleiro et al. 2018). Function ‘RAC_change’ quantifies changes in richness between a pair of time periods as the number of species gained or lost divided by the total unique species present in both time periods. We calculated these changes for each nursery between consecutive seasons in which samplings occurred (spring to summer, summer to autumn, and autumn to the following spring, over two years), totaling six inter-season comparisons. To determine if turnover in bee species was related to that of plants, we constructed LMMs with changes in bee richness (as gains or losses) as dependent variables, and the corresponding plant turnover metric, nursery type, category of inter-season comparison, and their interaction as fixed effects.

Finally, to determine if any bee species in our sweep netting or vane trap collections were significantly associated with local or landscape factors, we conducted a multi-level pattern indicator species analysis using function ‘multipatt’ (indicspecies) (De Cáceres and Legendre 2009). This analysis identifies species whose occurrence is associated with groups of sites categorized by a factor. We conducted analyses for season, nursery type, proportional native plant richness (which we dichotomized as > or ≤ 50% native), floral cover (> or ≤ the average percent cover of all plots), and proportional natural area within 1 km of the nursery (> or ≤ 50% natural area).

Descriptions of all sheets and variables in file are given in the first sheet ("metadata").