DOI: 10.6086/d1h094
title: Environment shapes the microbiome of the Blue Orchard Bee, Osmia lignaria

description: 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.
description: <b><i>Characterization of
 Study Sites </i></b>             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<sup>2</sup> to 15,525 m<sup>2</sup>, 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 <i>Osmia lignaria</i> 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 &lt;5% of the total area at
 each site were not included. We assessed land cover with spatial
 statistics tools in ArcGIS v. 10.1. <b> </b> <b><i>Bee Community Diversity
 Assessment  </i></b>             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.   <b><i>Bee
 Installation</i></b>             Over the course of three days
 in mid-March 2016, we installed <i>Osmia lignaria</i> 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 <b>°</b>C cold storage.   <b><i>Control
 Treatment</i></b>             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 <b>°</b>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.   <b><i>Illumina 16s
 Sequencing</i></b>             We collected
 a total of 344 <i>O. lignaria</i> (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).    <b><i>Parasite
 Detection</i></b>             We screened all <i>O.
 lignaria</i> bees with genus-specific primers for the presence of
 the neogregarine protozoan <i>Apicystis </i>spp., the
 trypanosomatid protozoean <i>Crithidia</i> spp, and the fungal
 <i>Aspergillus</i> 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.   <b><i>Data
 Analysis</i></b> 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 <i>Propionibacterium</i>) [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&lt;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
 <i>Osmia lignaria</i> 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 <i>Gammaproteobacteria</i>, order
 <i>Betaproteobacteriales,</i> genus
 <i>Lactobacillus</i>, and genus <i>Wolbachia</i>.
 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,
 <i>Lactobacillus</i>, and <i>Wolbachia</i> 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.