Wild bees encounter and collect environmental microbes whilst foraging. While environmental context affects bee diversity, little is known about it how affects the wild bee microbiome. We used field surveys in 17 urban gardens to examine whether and how variation in local and landscape habitat features shapes the microbiome of the solitary Blue Orchard Bee, Osmia lignaria. We installed O. lignaria cocoons at each site, allowed bees to emerge and forage, then collected them. We measured local features of gardens using vegetation transects and landscape features with GIS. We found that in microbiome composition between bee individuals varied by environmental features such as natural habitat, floral resources, and bee species richness. We also found that environmental features were associated with the abundance of bacterial groups important for bee health, such as Lactobacillus. Our study highlights complex interactions between environment context, bee species diversity, and the bee-associated microbes.

Characterization of Study Sites

            We examined local and landscape characteristics related to floral and nesting resource availability for bees in 17 urban gardens, ranging in size from 444 m² to 15,525 m², each separated by 2 km, across three counties (Monterey, Santa Clara, and Santa Cruz) in the California central coast. In two sampling periods, one month apart (early March and early April 2016), we measured local habitat characteristics within a 20 x 20 m plot placed at the center of each garden. We counted and identified all perennial trees and shrubs within the 20 x 20 plot. Then, in each plot we randomly selected four 1 x 1 m plots within which we counted all flowers (from annual crops, weeds, and ornamentals), and assessed percent ground cover of bare soil, grass, herbaceous plants, leaf litter, rocks, and mulch. We also estimated the total garden size. For analysis, values were averaged across the two sampling periods. In all, we measured 10 variables: % rock cover, % mulch cover, % leaf litter, % bare soil, % herbaceous plant cover, species richness of flowers, abundance of flowers, species richness of trees and shrubs, abundance of trees and shrubs, and garden size.

At the landscape scale, we classified the land cover types within a 500 m buffer surrounding each garden with data from the 2011 National Land Cover Database (NLCD, 30 m resolution) [41]. We selected 500 m buffers because while Osmia lignaria females have a maximum foraging distance of up to 1,200 m [42], they tend to collect more pollen and nectar at flowers near to their nests within 500 m. [43]. We created four land-use categories: 1) natural habitat (composed of deciduous (NLCD number 41), evergreen (42), and mixed forests (43), dwarf scrub (51), shrub/scrub (52), and grassland/herbaceous (71)), 2) open habitat (composed of lawn grass, parks, and golf courses (21)), 3) urban area (composed of low (22), medium (23), and high intensity developed land (24)), and 4) agriculture area (pasture/hay (81) and cultivated crop (82)). Other land cover types that covered <5% of the total area at each site were not included. We assessed land cover with spatial statistics tools in ArcGIS v. 10.1.

Bee Community Diversity Assessment  

            We used bee richness data collected for two prior experiments conducted by our research team at these sites [44,45]. Bee community diversity at each site was measured across six sampling periods between June and September 2013 and between June and September 2015. We used aerial nets for 30 minutes (not including handling time) each at site and three pan traps for 8 hours, and netted and placed traps within the 20 x 20 m vegetation plots. We identified bees using dichotomous keys to genus, and when possible, to species (see Quistberg et al. [44] for details on bee sampling and identification methods). For analysis, values were averaged across sampling periods and then across years for each site. In June 2016, we conducted one visual survey for bees for 30 min at each site to confirm that the relative ranking of species richness and abundance between sites was similar across years.

Bee Installation

            Over the course of three days in mid-March 2016, we installed Osmia lignaria bees at each community garden. We first placed one UV-sterilized binderblock laminate nest (Pollinator Paradise, Parma, ID) at or near the center of each site. Each binderblock was stocked with overwintering cocoons (100 female and 150 males). We applied three sprays of orchard bee attractant on each nest (Crown Bees, Woodinville, WA). Bees were allowed to emerge and forage for 16 days. We then collected adult female bees. Each bee was placed into a sterile 2 ml vial and immediately stored in dry ice. We also collected a blank, no-template control air sample at each site. We sterilized gloves, forceps, nets between collecting each bee and between each sampling site with bleach then ethanol. Bees were transported to the lab and into -80 °C cold storage.

Control Treatment

            To confirm that the environment confers unique bacterial communities to foraging bees, we allowed six female bees to emerge from their pupal cocoon casing in a sterile, indoor lab environment in petri dishes. Upon emergence, each female was immediately collected in a sterile 2 ml vial and placed into -80 °C cold storage. We also removed and collected an additional six females from their pupal casing by cutting individual pupal casings with a blade and removing the female with forceps. We sterilized gloves, blades, and forceps between samples.

Illumina 16s Sequencing

            We collected a total of 344 O. lignaria (an average of 19 bees per site). We extracted DNA from each sample and 1 control blank per site with the Qiagen ‘DNeasy Blood and Tissue’ extraction kit (Qiagen, Valencia, CA), but with the addition of tissue lysing step using sterile 5mm stainless steel beads and 0.1 mm glass beads in a Qiagen Tissue Lyser II to ensure extraction of gram positive bacteria [46]. We used whole-insect samples without surface sterilization because previously the bacterial community structure between surface sterilized bees and unsterilized bees has not been significantly different [47]. Library preparation and sequencing (Illumina MiSeq) was performed using previously published protocols [48]. To amplify the 16s rRNA gene, we used the 799F (5′-GAGT TTGATCNTGGCTCAG-3′) and 1115R (5′-GTNTTACNGCGGCKGCTG-3′) primer pair and included negative controls (control blank samples). 

Parasite Detection

            We screened all O. lignaria bees with genus-specific primers for the presence of the neogregarine protozoan Apicystis spp., the trypanosomatid protozoean Crithidia spp, and the fungal Aspergillus spp., a group that includes both common non-pathogenic environmental fungi and strains that manifest as Stonebrood disease in some bee species. We used parasite specific primers and conditions for genus-level identification (Table S1). Products were run alongside a standard ladder on a 1% agarose gel stained with GelRed to confirm amplicon size. Each assay included a negative and positive control.

Data Analysis

We used QIIME 2 (2019.1) to demultiplex and filter 16s sequence libraries [49].  First, we visualized and trimmed the low-quality ends of the reads with QIIME2, then used DADA2 [50] to remove chimeras, remove reads with more than two expected errors, and assign sequences into amplicon sequence variants (ASVs). We used the qiime2 feature-classifier [51] and trained to the 799-1115 region of the 16S gene with the SILVA database [52]. We filtered out features from the resulting ASV table that corresponded to contaminants as identified in our blanks (such as Propionibacterium) [53] or were present at only one read (singletons). To account for variable sequencing depth, we subsampled to 1,746 reads per sample. This allowed us to retain most samples (14 bee specimens excluded) and capture the majority of sequence diversity found in our samples. We used the MAFFT aligner [54] and FastTree v2.1.3 in QIIME 2 to generate a phylogenetic tree of our sequences [55]. We used the resulting tree and ASV to tabulate weighted UniFrac distance matrices (phylogenetic distances, weighted by abundance) [56] for beta diversity comparisons.

We used QIIME2 and the R statistical environment [57] to conduct and visualize analyses on alpha and beta diversity.  To avoid co-linearity between independent variables, we performed a variable selection process. We divided variables into four groups (floral resources, nesting resources, landscape composition, and bee community diversity) and ran Pearson’s correlations within groups to identify correlated (P<0.05) variables within groups. We selected variables that were correlated with the largest number of other variables in that group for subsequent analysis and that were previously found to be ecologically meaningful in describing pollinator diversity in the same field system [44]. We selected two variables reflective of floral resources at a site (abundance of annual flowers and the abundance of perennial trees and shrubs), one variable describing nesting materials (% bare soil, as soil is used by Osmia lignaria construct nests), one variable describing the landscape cover (% natural cover within 500 m), and one variable describing the bee community (bee richness). We natural log transformed variables that did not meet conditions of normality. To test for multicollinearity, we calculated a variance inflation factor (VIF) using the car package [58] and found each predictor had a VIF score below 2. 

            To compare the microbiome communities of bees allowed to forage and bees from our control treatment, we used QIIME2 to first test for homogeneity of dispersion using ‘permdisp’ and then performed a pairwise permutation-based, nonparametric test (PERMANOVA, permutations = 999) [59] on the weighted Unifrac distance matrices (which accounts for phylogenetic relatedness between ASVs). We visualized differences by plotting 95% confidence ellipses around each centroid on nonmetric multidimensional scaling (NDMS) plots using ‘ordiellipse’ function in vegan package [60] and ggplot2 [61]. To analyze how floral abundance, tree and shrub abundance, natural cover within 500 m, bare soil, and bee richness are related to the variance between microbiome communities of experimental foraging bees at different sites, we first performed non-metric multidimensional scaling weighted by abundance using metaMDS in the ecodist package in R [62]. We partitioned the weighted Unifrace distance matrices between sources of vegetative variation with the vegan package by applying ‘adonis’ on the dissimilarity distance matrix and used ‘envfit’ to fit floral abundance, tree and shrub abundance, natural cover, bare soil, and bee richness onto the ordination. We then compared the dissimilarity between vegetation communities and the bee microbiome communities at each site using a Mantel test, with 999 permutations. We first calculated dissimilarity between vegetation communities at each site using the Bray-Curtis method with the ‘vegdist’ function in the vegan package. The vegetation community matrix consisted of the following variables: % natural cover (500m), abundance of trees and shrubs, bare soil groundover (%), and flower abundance. We then calculated the dissimilarity between microbiome samples with the Bray-Curtis method. We obtained a Mantel statistic describing the correlation between matrices based on the Pearson method and plotted the correlogram with the mgram function in ecodist.

To analyze how site factors are associated with microbe diversity and parasite prevalence, for each bee we calculated relative ASV abundance for bacterial sequences found in the following groups: class Gammaproteobacteria, order Betaproteobacteriales, genus Lactobacillus, and genus Wolbachia. We used linear models with the lm function in R to examine relationships between bacterial groups and the site variables (floral abundance, tree/shrub abundance, % natural cover within 500 m, % bare soil, and bee richness). We tested combinations of these variables using the glmulti package [63], with bacterial abundance of Gammaproteobacteria, Lactobacillus, and Wolbachia as our response variables. For models where the AICc for top models was within 2 points of the next best model, we ran model averages with the MuMIn package [64]. When the best models shared the same significant predictors as model averages, and we reported output from best models.

We indicated parasitism for each bee specimen with binary (0/1) prevalence data for each parasite. To examine the relationship between the prevalence of each parasite and the (scaled) abundance of each bacterial group, we used a generalized linear mixed effect model with a random effect of site and a binomial distribution using the glmer function in R.