Pollinator declines can leave communities less diverse and potentially at increased risk to infectious diseases. Species-rich plant and bee communities have high species turnover, making the study of disease dynamics challenging. To address how temporal dynamics shape parasite prevalence in plant and bee communities, we screened >5,000 bees and flowers through an entire growing season for five common bee microparasites (Nosema ceranae, N. bombi, Crithidia bombi, C. expoeki and neogregarines). Over 110 bee species and 89 flower species were screened, revealing 42% of bee species (12.2% individual bees) and 70% of flower species (8.7% individual flowers) had at least one parasite in or on them, respectively. Some common flowers (e.g., Lychnis flos-cuculi) harboured multiple parasite species, whilst others (e.g., Lythrum salicaria) had few. Significant temporal variation of parasite prevalence in bees was linked to bee diversity, bee and flower abundance, and community composition. Specifically, we found that bee communities had the highest prevalence late in the season, when social bees (Bombus spp. and Apis mellifera) were dominant and bee diversity was lowest. Conversely, prevalence on flowers was lowest late in the season when floral abundance was the highest. Thus, turnover in the bee community impacted community-wide prevalence, and turnover in the plant community impacted when parasite transmission was likely to occur at flowers. These results imply that efforts to improve bee health will benefit from promoting high floral numbers to reduce transmission risk, maintaining bee diversity to dilute parasites, and monitoring the abundance of dominant competent hosts.

Bee and flower collections

To determine parasite prevalence within bees and upon the surfaces of flowers, we collected samples weekly across three old-field sites in upstate New York, between 18^th April 2017 and 22^nd September 2017. The sites were named Lansing (Lat: 42° 32' 24.4932'' N, Long: 76° 29' 47.9076'' W), McDaniels (Lat: 42° 32' 11.5872'' N, Long: 76° 25' 3.7668'' W), and Whipple (Lat: 42° 29' 23.6328'' N, Long: 76° 25' 49.818'' W). The McDaniels and Whipple sites are managed by Cornell University and required no permits for their use. The Lansing site is privately owned and we obtained permission to use the site. The distance between sites ranges from 5.3km to 7.5 km. Weather depending, we dedicated a day of collecting (~7 hours) per site, every week. On a typical collection day 3 members would be collecting samples from the same site (21 person-hours per site per week). In total we spent 58 days collecting samples over the five-month period (equivalent of approx. 1,200 person-hours). We collected bees in sterile tubes either directly or with the use of a net. Variety of floral form meant that ‘flowers’ were either single flowers or inflorescence depending on the plant species. The number of florets used per inflorescence was the same number a bee would forage upon in a single interaction with that species. We placed flowers that the bees were foraging upon directly into separate sterile tubes using sterilized forceps. Bee and flower samples were placed on dry ice immediately after collection, then transported to the lab and stored at -80C until we commenced DNA extraction.

 

Floral abundance and diversity

Each week throughout the growing season, we randomly selected three quadrats (each 10 m ´ 10 m) at each site, within which we identified the floral species in bloom and counted the number of stems for each (two out of 26 weeks were not surveyed due to weather). Per species, we estimated the number of floral units per stem and averaged this across the sites and weeks the species was in bloom (minimum five measurements per species). We define a floral unit as being the typical unit (single flower or inflorescence) that a bee typically foraged from ¹. This definition of a floral unit was also consistent with the amount of plant material we used for each sample during parasite screening. We therefore define floral abundance of each plant species within the quadrat as the product of stem counts with the estimated number of floral units per stem. The floral survey data are found in Supplementary Table 2.

 

Plant and bee identification

We identified plants either in the field or via specimens and/or photos brought back to the lab and keyed out ^2–5. We confirmed the identity of each individual plant screened for parasites before DNA extraction occurred. Overall, we screened 2,624 flowers from 89 species for parasites (Supplementary Table 5). We identified bees after the gut dissections occurred for each specimen using reference materials located in the Cornell University Insect Collection (CUIC: http://cuic.entomology.cornell.edu/) and published keys ^6–9. All identifications were conducted by P. Muñiz, taxa verifications were conducted by M. Arduser, and all voucher specimens are housed in the McArt lab or the CUIC. Overall, we screened 2,672 bees from at least 110 species for parasites (Supplementary Table 6).

 

Bee parasites 

We selected to screen samples for parasites with known effects on bee health and links to population declines ^10,11. Historically, research has focused on parasites of honeybees and a small number of bumblebee species. This research bias has contributed to a fundamental lack in our understanding of the full host ranges of bee parasites.  Because we were interested in screening across bee taxa, we specifically sought to screen for parasites that have been identified in multiple bee species. We considered common microparasites of bees as three groups; microsporidians, trypansomatids, and neogregarines. The microsporidians are comprised mostly of Nosema species whose spores are transmitted via the faecal-oral route ¹². The effects of microsporidia infections in honeybees and bumblebees includes wing deformity, reduced foraging efficiency, reduced colony fitness and increased mortality, with N. bombi being found in multiple species of bumblebees and N. ceranae found in bumblebees and honeybees  ^13–19 In addition N. ceranae may be able to infect the mason bee, Osmia bicornis²⁰. While this study did not find impacts on survival, a similar study that assessed impacts of N. ceranae on larval O. bicornis did find negative impacts on survival²¹. The trypanosomatids are mostly Crithidia species whose cells are transmitted via the fecal-oral route. The effects of trypanosomatid infection in bumblebees include reduced foraging efficiency, reduced queen fitness, and increased mortality of infected bees. Crithidia spp. can infect bumblebees, honeybees and several solitary bee species ^22–30. The neogregarines are an understudied group, with one described bee parasite, Apicystis bombi. This neogregarine has been detected in honeybees, a range of bumblebee species and solitary bees. Infected bumblebees have reduced fat bodies, increased mortality, and queen bumblebees are less likely to survive hibernation ^31–37.  

 

Parasite DNA extraction from individual bees and inflorescences

Parasite detection does not confirm an active infection, and although we endeavoured to reduce the likelihood of reporting uninfected bees by eliminating parasites on the outer cuticle via surface sterilization, and only report parasites in tissues known to harbour the selected parasites, some of the detections may be due to transient parasite material and not an active infection ³⁸. We performed dissections in bleached, UV sterilized hoods, and sterilised instruments between samples using a dry bead steriliser set to 250°C. We carefully removed the alimentary canal from the mid-gut to the rectum using standard techniques ³⁹. If the gut broke apart inside the bee, in addition to dissecting out the torn gut, we washed 10ml of Phosphate-buffered saline in and out of the bee cavity with a pipette to recapture any spilled gut contents before adding it to the gut tissue for DNA extraction. We performed dissections with minimal destruction to the cuticle to allow accurate species identification.

 

Similarly with parasite detections in bees, we do not know if parasites molecularly detected on flowers are viable, or pathogenic to all foraging bees. We term flowers as potential ‘transmission hubs’ because of their role in providing a physical platform for microbes to be deposited by bees and then acquired by subsequent foraging bees^40–42. Using the 96-well plate Qiagen DNeasy blood and tissue Extraction Kit protocol, we washed single flowers/inflorescences in 600ml ATL buffer by pulse vortexing for 30 seconds. We then transferred 450ml of the wash to a 2ml screwcap tube with ~100 µL of 0.10mm zirconia beads and one 5mm steel bead. The wash was then lysed for 30 seconds at 6.5M/s on a Omni Bead Ruptor 24 homogenizer. For bee guts, we added 180 µl buffer ATL, ~ 100 µL 0.1 mm zirconia beads, and one 5 mm steel bead to each sample. We then homogenized the bee guts at 30 hz for 3 minutes on a Qiagen Tissue Lyser II. For both flower washes and bee gut homogenate, we next added 50 µL of Proteinase K before allowing the samples to incubate at 56°C overnight. Following incubation, we followed the standard “Quick Start" protocol provided by Qiagen, until the DNA was eluted in 100 µL of AE buffer.

 

Parasite screening

We took a 2-step approach to determine the parasite prevalence on flowers and in bee guts; (1) we screened for the presence of two broad taxonomic groups (trypanosomatids and microsporidians) known to contain multiple parasite genera, before (2) screening these positively identified groups for likely species of bee parasites within those groups. In addition, the multiplex PCR also identified samples containing neogregarine parasites. The broad multiplex panel was designed to efficiently screen flowers and insect guts for the most common bee-infecting parasites within a single reaction ⁴³. Samples positive for Crithidia or Nosema were then diagnosed with species-specific multiplex panels (C. bombi, C. expoeki, and N. bombi, N. ceranae, respectively). Primers were either newly designed (see supplementary methods) or chosen from existing literature ^43–45. The concentrations of each reagent and the thermocycling conditions are given in Supplementary Table 1. PCR products were run alongside a size standard on a 2%  agarose gel, stained with GelRed® to visualise and confirm amplicon size. Each assay included a negative and a positive control. The raw data for broad range and species-specific parasite screens can be found in the Supplementary Table 3 and 4. Furthermore, a subset of positive samples had their amplicons Sanger sequenced to confirm correct amplification (Sequence details can be found on Dryad and via NCBI database with accession numbers MT212154-MT212159, MT296581-MT296586, MT302779-MT302784, MT359894-MT359896, MT366919, MT387450 and MT387451).  

 

Statistical analysis

1. Parasite prevalence of the bee community through time

For each of the broad taxonomic groups and parasite species that we screened for, temporal trends in the overall parasite prevalence of the bee communities were evaluated using generalized linear mixed models (GLMM), with parasite status of the samples as binomial response, week number within the field season as the explanatory variable, and site as a random factor. All subsequent GLMMs also used site as the random factor. Many issues can potentially arise from the low parasite prevalence in the data (e.g., prevalence less than 1% for some parasite species); for instance, maximum likelihood estimates (MLE) are only asymptotically unbiased, meaning that they can be significantly biased if the number of positives is small ⁴⁶. Therefore, we excluded any parasite groups or species with less than 20 positives in the analysis, which in this case meant leaving out N. bombi. The same criterion was used in all subsequent analysis involving parasite status. In addition, we repeated the analysis for the four species and neogregarines combined as a single group of bee parasites. The broad groups microsporidia and trypanosomatids were not included in this combination since they also contain species that may not infect bees. To account for multiple testing given the number of parasite species/groups being considered, throughout this manuscript we use a Bonferroni-corrected significance level of a = 0.05/N_p, where N_p is the number of parasite species/groups considered in each analysis. To justify the use of GLMMs instead of autoregressive time-series models, we checked for temporal auto-correlation using Durbin-Watson tests on scaled residuals.

 

Next we investigated bee community composition and diversity as possible drivers of any observed temporal trends. These are discussed in details below.

a. Bee community composition

Should there be significant heterogeneity in host competence among taxa, community turnover may drive changes in parasite prevalence. As the four most common bee genera Apis, Bombus, Ceratina and Lasioglossum comprised as much as 76% of the bee samples, we explored how the individual contributions of these genera (as well as all other genera as a single group) to the overall parasite prevalence of the community varied over time. First, temporal trends in the relative abundance of these genera were visualized with smoothing splines fitted using multinomial generalized additive models (GAM) with site included as a random factor. Relative abundances were further quantified by partitioning the season into three equal windows of 8 weeks each, and performing for each site and window an exact multinomial goodness-of-fit test (against the hypothesis of equal multinomial proportions), followed by post-hoc tests (with Bonferroni-corrected a given the number of genera) to identify the genera that fall above or below the expected proportions.

Second, to evaluate heterogeneity in parasite prevalence across genera, binomial GLMMs were again fitted for each parasite, but this time with genus and its interaction with week number as additional explanatory variables. Main effects of genus were evaluated regardless of whether the interactions were significant. Testing for main effects in the presence of significant interactions is known to violate the principle of marginality ^47,48. More specifically, the main effects of genera are given by the differences in fitted logit links at week 0 between genera (i.e. differences in y-intercepts); since significant interactions imply different temporal trends between genera (different slopes), the main effects of genera become ambiguous since they now depend on what date we choose to assign as week 0. To address this, we defined week 0 of each genus to be the median week number among all samples of that genus; week numbers of the samples were then shifted accordingly (the shifts were not shown in Supplementary Fig. 2 to avoid confusion).  By doing so, the main effects could be unambiguously interpreted as the difference in parasite prevalence between genera, evaluated for each genus at a characteristic date of that genus. Post-hoc tests of pairwise contrasts (with multiple testing corrections for the number of contrasts) were performed whenever the main effects of genus turned out to be significant. Significance of temporal trends in the parasite prevalence of each genus were also assessed (with multiple testing corrections for the number of genera).

Finally, temporal trends in the individual contributions of each genus to the overall parasite prevalence of the community (relative abundance ´ parasite prevalence of genus) were visualized using smoothing splines fitted using GAMs with site as a random factor.

 

b. Bee diversity

Should bee diversity affect parasite prevalence, temporal trends in the former may drive that of the latter. To assess how the bee diversity varied across the season, we calculated the weekly Shannon index at each site, based on the collected bee samples, and assessed the temporal trend using a linear mixed model. Before fitting the model, rarefaction was implemented for the bee Shannon indices to address non-uniformity in bee collection efforts. Note that rarefaction tended to reduce larger indices more than smaller ones, hence potentially affecting the strength of any temporal trends. Therefore, we evaluated a range of subsample sizes, rather than simply have the size be determined by the smallest site/week sample. Each time, samples smaller than the size being considered were discarded. This allowed us to identify an optimal size large enough for most indices to remain close to their non-rarefied values, but yet small enough to minimize the number of discarded samples. All analyses were conducted using rarefied indices at this optimal size; robustness of any results to the choice of subsample size were also assessed. We also assessed whether there were differences in bee diversity between sites using one-way ANOVA.

Next, we evaluated the relation between parasite prevalence and diversity using a GLMM with parasite status of the weekly bee samples as binomial response, and rarefied Shannon indices calculated using the samples as the explanatory variable. To explore whether any observed associations were entirely driven by temporal trends in the abundance of Apis and Bombus, we also repeated the analysis, but this time without Apis and Bombus samples when calculating parasite prevalence.

 

2. Parasite prevalence of the floral community through time

As was done for bees, temporal trends in the overall parasite prevalence of the floral community were evaluated using GLMMs.  Flowers were not obtained at random, instead their collections were directed by the co-collection of foraging bees upon them. Hence the floral community being screened should be thought of as being biased to some extent by bee foraging preferences – unlike the floral surveys used for abundance/diversity analyses. N. ceranae was excluded from the analysis since less than 20 flower samples tested positive.

Temporal trends in the total floral abundance within quadrats were also evaluated using GLMMs, with counts as negative binomial response. Negative binomial was chosen to allow for overdispersion due to clustering. To investigate whether floral abundance can dilute parasite prevalence on flowers, we evaluated the relation between parasite prevalence and floral abundance using a binomial GLMM, with parasite status on flower samples as binomial response and log₁₀(mean total floral abundance per quadrat) as explanatory variables. We chose to take the logarithm of floral abundance, first because of model convergence issues, and second because it was the more appropriate linear predictor for the logit link if we assumed prevalence to be inversely proportional to floral abundance.

Finally, temporal trends in the floral diversity (Shannon index based on the quadrat floral abundance surveys) were also evaluated using linear models, and the relation between parasite prevalence on flowers and floral diversity evaluated using a binomial GLMM.

 

Software and packages used

All analyses were performed in R v3.5.1 ⁴⁹. Packages used were lme4 ⁵⁰ and glmmTMB ⁵¹ for fitting GLMMs, DHARMa ⁵² for performing Durbin-Watson tests using scaled residuals, DescTools ⁵³ for generating binomial and multinomial confidence intervals in the plots, mgcv ⁵⁴ for fitting smoothing splines using GAM, XNomial ⁵⁵ for exact multinomial tests, multcomp ⁵⁶ for pairwise contrasts, and Vegan 2.54 for calculating Shannon indices ⁵⁷.  References attached as text file.