Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This research was supported by the National Honey Board. In addition, the Flenniken Laboratory is supported by the National Sciences Foundation CAREER Program, the United States Department of Agriculture National Institute of Food and Agriculture, Agriculture and Food Research Initiative (USDA-NIFA-AFRI) Program, Montana Department of Agriculture Specialty Crop Block Grant Program, the National Institutes of Health IDeA Program COBRE grant GM110732, National Science Foundation EPSCoR NSF-IIA-1443108, Hatch Multistate Funding (NC-1173), Project Apis m., the Montana State Beekeepers Association, Montana State University, and the Montana State University Agricultural Experiment Station. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
1,2]. California almond production is the most renowned example of the role of honey bee pollination services in the US. Every year, over 60% of the approximately 2.5 million commercially managed honey bee colonies in the United States are transported to the Central Valley of California to pollinate the almond crop. California almond growers produce ~ 80% of the global almond supply, which was valued at $1.7 billion in 2014 . Honey bee colonies are important for the pollination of plants in both agricultural and non-agricultural landscapes [4,5]. Thus, high annual honey bee colony losses in the US, which have averaged ~ 33% annually since 2006 and increased from historic levels of approximately 12%, are concerning to stakeholders, including beekeepers, growers, and scientists [6,7]. To compensate for increased colony losses and meet pollination service demands, beekeeping operations have been splitting their colonies (i.e., producing two colonies from one colony) more frequently. Therefore, the total number of commercially managed honey bee colonies in the US has been maintained at approximately 2.5 million since the late 1990s  despite unsustainable annual colony losses.
7,9–18]. Honey bee colony health is typically evaluated by estimating colony population size [19–23]. Colonies are comprised of ~35,000 sterile female workers, hundreds of males (drones), and a single reproductive queen that can lay approximately 1,000 eggs per day . While no single factor has been deemed responsible for high annual honey bee colony losses, pathogens are major contributing factors [6,7,10,11,21,22,25–33].
33,34]. Honey bee infecting viruses include Acute bee paralysis virus (ABPV), Black queen cell virus (BQCV), Deformed wing virus (DWV), Israeli acute paralysis virus (IAPV), Kashmir bee virus (KBV), Sacbrood virus (SBV), Chronic bee paralysis virus (CBPV) [34–38] and the Lake Sinai viruses . In addition to viruses, pathogens of honey bees include eukaryotes such as the trypanosomatid Lotmaria passim (formerly Crithidia mellificae strain sf [39,40]), the microsporidial pathogen Nosema ceranae , and bacterial pathogens such as Paenibacillus larvae  and Melissococcus plutonius , the causative agents of American and European foulbrood diseases, respectively. In addition, the ectoparasitic mite, Varroa destructor, contributes to decreased colony health by feeding on developing bees (brood) and facilitating virus transmission  (i.e., DWV [14,45], KBV [46,47], and IAPV [34,44,48]). Mite parasitization of developing honey bees can result in physical deformities, reduced body weight, greater DWV virus levels, and/or death [44,49].
21,22,25,26,31]. Furthermore, there have been relatively few studies that have monitored individual commercially managed colonies [7,10,15–17,21,22]. The majority of studies in the US have monitored honey bee health at the apiary level with the main goals of detecting exotic pathogens and establishing a baseline understanding of honey bee colony health [6,7,15,17,30]. Additional studies that monitor colonies located in varying geographic locations are needed to better understand the impacts of multiple factors on honey bee colony health. To address this knowledge gap, we monitored Minnesota-based commercially managed honey bee colonies involved in the 2014 California almond pollination event. The goal of this study was to investigate the association between multiple factors, including pathogen seasonality, pathogen abundance, beekeeping operation, colony population size, and level of mite infestation, on honey bee colony health and colony losses. We determined that sampling date strongly influenced pathogen prevalence and abundance and confirmed a strong association between DWV abundance and Varroa destructor mite infestation. However, the abundance of other viruses (i.e., BQCV, LSV1, and LSV2) was not associated with mite infestation. Together, data from this and other temporal observational cohort studies that precisely account for sampling date (i.e., day of the year) and monitor individual honey bee colonies will lead to a better understanding of the association between factors affecting pathogen abundance and honey bee colony health.
S3 Table) [20,21]. Samples of honey bees were obtained from monitored colonies, up to four times each. Specific 2014 sample dates varied and were categorized into discrete sampling events in order to better compare data within and between colonies managed by two beekeeping operations. Specifically, “before almond pollination” samples were obtained from colonies located in California in January 2014 (Beekeeping Operation 1 on Jan. 27 and Beekeeping Operation 2 on Jan. 29). Samples that represent the pathogens associated with honey bee colonies “during almond pollination” were obtained immediately post-almond pollination in March 2014 (Beekeeping Operation 1 on March 3 and Beekeeping Operation 2 on March 14). “After almond pollination” samples were obtained from colonies that were located in Minnesota in early (June 9) and late (Aug. 25 and Sept. 2) summer (Fig 1 and S3 Table). Additionally, core samples containing pollen/bee bread, wax, and honey were obtained from a representative frame in each colony and analyzed for pesticide residues, as described in Kegley et al. 2017 . The sample cohort analyzed herein includes pathogen prevalence, abundance, and survivorship data obtained in 2014, including variable numbers of samples obtained from dead (n = 2), weak (n = 8), average (n = 34), and strong (n = 48) rated colonies (total = 92) spanning the duration of the 2014 California almond pollination season, and colony longevity data through January 2015 (Fig 1 and S3 Table). In total, 21 colonies in this sample cohort died by the end of January 2015.
Fig 1. Commercially managed honey bee colonies were longitudinally monitored before, during, and after the 2014 almond pollination season.
22]. The following equation from Pirk et al. 2013, N = ln(1-D) / ln(1-P) (N = sample size, ln = natural logarithm, D = probability of detection, P = proportion of infected bees) predicts that with a sample size of five bees, pathogenic infections affecting 45% or more of the individuals within a colony would be detected with 95% probability ; this sample size has proven sufficient for the pathogen-specific PCR detection of highly prevalent pathogens [21–23].
52]. At each sampling event, 150 adult bees were obtained from a brood-containing frame and rinsed in ethanol to dislodge attached mites into a collection jar. Mites infestation percentage was reported as the number of mites per 100 bees . In total 92, samples were included in the analysis (Fig 1 and S3 Table).
S1 Table. Specifically, PCR was used to determine the prevalence of DWV, BQCV, IAPV, KBV, LSV1, LSV2, N. ceranae, and L. passim (S1 Table). In brief, each PCR included 1 μl cDNA template, combined with 10 pmol each of forward and reverse pathogen specific primers, and amplified with ChoiceTaq polymerase (Denville) according to the manufacturer’s instructions using the following cycling conditions: 95°C for 5 minutes; 35 cycles of 95°C for 30 seconds, 55°C for 30 seconds, and 72°C for 30 seconds, followed by final elongation at 72°C for 4 minutes. The PCR products were visualized using 2% agarose gel electrophoresis followed by fluorescence imaging. In addition, the PCR assays utilized in this study were verified by sequencing . Positive and negative control reactions were included for all pathogen-specific PCR analyses and exhibited the expected results. IAPV was detected by PCR in only one of the samples, thus IAPV abundance was not examined by quantitative PCR.
S1 Table). All qPCR reactions were performed in triplicate with a CFX Connect Real Time instrument (BioRad) and the following reaction conditions: 2 μL of cDNA template in 20 μL reactions containing 1X ChoiceTaq Mastermix (Denville), 0.4 μM each forward and reverse primer, 1X SYBR Green (Life Technologies), and 3 mM MgCl2. The qPCR thermo-profile consisted of a single pre-incubation 95°C (1 minute), 40 cycles of 95°C (10 seconds), 58°C (20 seconds), and 72°C (15 seconds). Plasmid standards, containing from 109 to 103 copies per reaction, were used as qPCR templates to assess primer efficiency and quantify the relative abundance of each pathogen. The linear standard equation for primer efficiency of each primer set, which was generated by plotting the crossing point (Cp) versus the log10 of the initial plasmid copy number, is as follows: BQCV: y = -3.7336x+ 42.849, R2 = 0.996; DWV: y = -3.443x + 41.277, R2 = 0.99958; LSV1: y = -3.1994x + 38.71, R2 = 0.982, LSV2: y = -3.8147x + 44.805, R2 = 0.980; KBV: y = -3.4505x + 40.099, R^2 = 0.99927, L. passim: y = -3.4825x + 41.025, R2 = 0.98588 and N. ceranae: y = -3.9656 + 42.124, R2 = 0.9961; the minimum qPCR detection levels were ≤ 1,000 copies per sample. In addition, qPCR of a host encoded gene, Apis m. Rpl8, was performed using 2 μL cDNA template on each qPCR plate to ensure consistency and cDNA quality. qPCR products were analyzed by melting point analysis and 2% agarose gel electrophoresis. In addition, the qPCR assays utilized in this study were verified by sequencing . An estimate of the number of viral RNA copies per bee can be obtained by multiplying the reported qPCR copy number values by 25. This estimate is based on the following: typical RNA yield was approximately 50 μg per bee, each qPCR reaction was performed on cDNA generated from 2 μg RNA, therefore each well represents 1/25th of an individual bee .
22]. The natural logarithm (ln) was used to transform pathogen abundance (as determined by qPCR) to normalize the response variable; 1 was added to each observation since some observations had 0 total abundance. The confounding effects of seasonality on pathogen abundance was removed by using the residuals of a linear regression evaluating the relationship between pathogen abundance and the “day of year” of sampling as the response variable in the multiple linear regression models. As a result, models reported associations between the residuals of pathogen abundance after accounting for the effects of pathogen seasonality (“residual abundance”), and colony strength rating, beekeeping operation, and mite count. Most colonies were measured multiple times in the 92 samples; therefore, a random effect for individual colony was included in the models to account for repeated measures. Estimated models for each pathogen are described as:53]. The main effects as well as an interaction between sampling event and mite count were included in the reported PERMANOVA model. Colony health rating was included as a main effect in a separate PERMANOVA model, but was not considered a significant predictor. Additionally, we described any differences in the variability of pathogen composition by sampling event using an analysis of multivariate homogeneity of group dispersions (betadisper) . Furthermore, we used a similarity percentage (SIMPER) analysis to describe the percent dissimilarity explained by each pathogen in comparisons between consecutive sampling events .
54]. Generalized linear models were conducted using the nlme package . NMDS plots were generated using the labdsv package while all other multivariate analyses were performed using the vegan and labdsv packages [53,56].
Fig 1) . At each sampling event, colony health, using colony population size as a proxy for health, was monitored and samples of live honey bees were obtained for pathogen diagnostics, including Varroa destructor mite infestation levels and the prevalence and abundance of five commonly occurring viruses (i.e., BQCV, DWV, LSV1, LSV2, and KBV), and two eukaryotic pathogens (the trypanosomatid Lotmaria passim and the microsporidial pathogen Nosema ceranae) using pathogen specific PCR and qPCR, respectively [22,23,26,50]. Israeli acute paralysis virus (IAPV) was only detected in one sample using PCR, thus it was not included in further analysis. In this work, we use “pathogen prevalence” to indicate the number of different pathogens detected in a sample and “pathogen abundance” as the estimated number of relative RNA transcripts for each pathogen using RT-qPCR (S3 Table). Furthermore, we use “pathogen composition” to refer to the abundance of all co-infecting pathogens in a sample relative to all other samples. Therefore, both pathogen prevalence and abundance affect the pathogen composition of each colony throughout the course of the study.
S3 Table) [7,52,57,58]. Mean pathogen prevalence varied by colony health rating (ANOVA, F3,88 = 3.05, p-value = 0.03). Specifically, samples from dead colonies averaged 6.50 pathogens (n = 2, SE +/- 1.50), weak colonies 5.25 pathogens (n = 8, SE +/- 0.38), average colonies 5.38 pathogens (n = 34 SE +/- 0.18), and strong colonies 5.95 pathogens (n = 48, SE +/- 0.13) (S2 Table and S1 Fig). Strong rated colonies had higher pathogen prevalence than average rated colonies (S1 Fig). No other pairwise comparisons of mean pathogen prevalence between colony health ratings were statistically different.S3 and S4 Tables and S2 Fig). Specifically, colonies sampled early in the year before almond pollination averaged 5.00 pathogens (n = 28, SE +/- 0.17), and colonies sampled in March, immediately after almond pollination and thus representing the pathogen profile of colonies during the almond bloom, averaged 5.39 pathogens (n = 28, SE +/- 0.09). Colonies sampled in early summer when the colonies were in Minnesota (i.e., after 1, June) averaged 7.09 pathogens (n = 11, SE +/- 0.28), in late summer (i.e., after 2, Aug.) colonies averaged 6.43 pathogens (n = 14, SE +/- 0.13), and at the end of summer (i.e., after 3, Sept.) colonies averaged 6.00 pathogens (n = 11, SE +/- 0.23) (S4 Tableand S2 Fig). Results from Tukey’s HSD post-hoc tests indicate that the mean pathogen prevalence was greater in sampling events later in the year, when the colonies were located in Minnesota (after 1, 2, and 3) compared to sampling events in January and March, before and during almond pollination, respectively, when the colonies were located in California (S2 Fig).
S3 Table). Quantitative PCR values ranged from 0 to 6.7x108 per sample (S3 Table). Overall, longitudinal assessment of the abundance of common honey bee pathogens at the colony level from January to September 2014 revealed that individual colonies, both within a single beekeeping operation and between two typically managed commercial beekeeping operations, varied from the beginning of the monitoring period in January, until the end of the monitoring period in September (Fig 2 and S3–S9 Figs).
Fig 2. Longitudinal assessment of honey bee pathogen abundance at the colony level.
S3 Table). The natural log transformed data for each pathogen, BQCV (pink), DWV (orange), Lp (light pink), Nc (gray), LSV1 (light blue), LSV2 (dark blue), and KBV (black), is represented by the height of each segment of the vertical bar representing a sample from an individual colony (y-axis); the height of the bar represents the total pathogen abundance for each colony over the course of the study (January to September 2014, x-axis).Fig 3). Results from a PERMANOVA analysis indicated that sampling date explained the most variation in pathogen compositions (Fig 3 and Table 1, R2 = 0.61, p < 0.01), followed by levels of Varroa destructor mite infestation (Table 1, R2 = 0.01, p = 0.05), indicating pathogen composition varies seasonally and with changing levels of mite abundance . Unexpectedly, pathogen composition did not differ between colony health ratings (PERMANOVA, F3,83 = 1.467, p = 0.14, R2 = 0.018). To identify the specific pathogens that contributed the most to the changes in pathogen composition between sampling events, a SIMPER analysis was used to calculate the contribution of each pathogen to the dissimilarity of pathogen compositions between consecutive sampling events (Table 2) . Differences in pathogen composition of honey bee colonies between consecutive sampling events was not consistently explained by a single pathogen, but varied for each comparison, highlighting the dynamic nature of honey bee colony pathogen composition. DWV contributed the most to the differences between pathogen compositions of honey bee colonies when the mean difference in V. destructor mite infestation was greatest (Table 3), likely due to the role of mites in DWV transmission. Results from a multivariate homogeneity of group dispersions analysis (betadisper) indicated that the dispersion of pathogen composition among sampling events did not change (F4,87 = 0.27, p-value = 0.90) . Together, this analysis indicated that the composition of honey bee associated pathogens is most influenced by sampling date and the level of mite infestation.
Fig 3. Relative pathogen composition of honey bee colonies visualized by sample event.
Table 1 indicated that the sampling event explained the most amount of variance in the pathogen composition.
Table 1. Analysis of variance of pathogen composition.
Table 2. Similarity percentage analysis of the relative pathogen composition between consecutive sampling events.
Table 3. Pathogen abundance model summaries.S3 Table). Whereas, honey bee samples obtained from colonies later in the year, when the colonies were located in Minnesota after almond pollination (i.e. after1, after2, and after3), harbored a 7.17% mean mite infestation (mean = 7.17, SE+/- 1.21, n = 36) (S3 Table). To evaluate the risk Varroa mites pose to colonies throughout the sampling period, results from a generalized linear mixed effects model with a binomial family distribution demonstrate that the odds of a colony exceeding the recommended treatment threshold level of 3% Varroa mite infestation is e0.26 (1.3) times higher with each increase in the day of the year (log-odds of surpassing threshold, est = 0.23) (Fig 4and Table 3). These results indicate that honey bee colonies within this sample cohort frequently accumulated higher levels of mites toward the end of the pathogen monitoring period (Fig 4), consistent with the typical mite population growth pattern that tracks the seasonal colony buildup and decline over the course of a year.
Fig 4. The probability of honey bee colonies accumulating Varroa destructor mite densities above the threshold recommended for treatment is greater later in the year.S3 Table). It is notable that both beekeeping operations had numerous colonies with high levels of mite infestation, despite application of several anti-mite treatments including amitraz, oxalic acid, and formic acid over the course of the sampling period (see Methods).
S3–S9 Figs), therefore we further examined the relationship between pathogen abundance and other factors (e.g., beekeeping operation, colony health, and mite infestation) after removing the confounding effects of sample date. In order to remove the confounding effects of sample date on pathogen abundance, we examined the relationship between the residuals of a linear regression of pathogen abundance in response to “day of year” (i.e., residual abundance), and used residual abundance as the response variable in multiple linear regression to evaluate associations with residual abundance and beekeeping operation, mite infestation percentage, and colony health, which was represented in the models as frame count (Table 3).2 and 5). All the honey bee samples in this sample cohort (n = 92) were tested for the presence of DWV using PCR (S3 Table). Interestingly, none of the samples collected in January and March when the colonies were located in California (n = 56) were DWV positive, whereas 96.4% of samples obtained from colonies located in Minnesota (n = 36) were positive for DWV. Quantitative PCR was utilized to determine the relative DWV abundance in each sample (n = 92), which in general increased over the course of the monitoring period (Fig 5 and S3 Table). Multiple linear regression was used to investigate associations between residual DWV abundance with the predictor of interest, frame count (i.e., colony strength), while also accounting for beekeeping operation and level of mite infestation, indicating that median residual DWV abundance increased by e0.23 (126%) for each unit increase in mite infestation after accounting for beekeeping operation and frame count (Table 3). These data support previous results describing the relationship between mite infestation and DWV abundance [7,14,45].
Fig 5. DWV abundance increased from January to September 2014.Table 3). While not significant, a negative association was observed between residual KBV abundance and levels of mite infestation (mites, est = -0.09, SE = 0.06, p-value = 0.06) (Table 3). Mites transmit both DWV and KBV, therefore the positive correlation between DWV and mite infestation was expected, whereas further investigation is required to determine the relationship between KBV and mite infestation [44,46,47].21–23,31,60,61]. Consistent with other studies, LSV2 abundance was higher in weak colonies, as compared to healthy colonies with greater bee populations in this sample cohort [23,28], but this observation was only statistically supported when colony health was categorized into health ratings based on frame count. The median residual LSV2 abundance was e7.39 (161970.60%) greater in “dead” colonies compared to “weak” rated colonies (t-test, t-value56 = 3.02, p = 0.004), when categorical colony health data, as opposed continuous frame count data, was used in the models, while also accounting for the effect of beekeeping operation and mite infestation. However, samples were only collected from 2 of the 22 colonies that died. Despite the close phylogenetic relationship of LSV1 and LSV2 , changes in LSV1 abundance were associated with monitored factors that were different from those that correlated with LSV2 abundance. Although LSV1 was similarly widespread, testing positive in 86.9% of samples (80/92) (S3 Table), median residual LSV1 abundance was e2.27 (967.94%) times greater in Operation 2 compared to Operation 1 when also accounting for frame count and mite infestation (Operation 2, est = 2.27, SE ± 0.95, p = 0.02) (Table 3) (S6 Fig).Table 3). The median residual abundance of L. passim increased by e0.23 (125%) with each unit increase in frames counted, and the median residual abundance of N. ceranae increased by e0.18 (119.72%) in response to each unit increase in frames counted when also accounting for the effects of beekeeping operation, and mites (Table 3). These results differ from previous studies demonstrating that N. ceranaeinfluences co-infections of other pathogens important to reduced colony health in the US [10,28,49], or as the causative agent of colony loss in Spain [62,63].S10 Fig). For all pathogens except KBV, there was no correlation between abundance and Δ-Frame Count independent of day of year (S10 Fig). There was a negative correlation between KBV abundance and Δ-Frame Count independent of day of year (S10 Fig).
S3–S9 Figs). Together these data illustrate the similarities and differences in pathogen loads over the course of this study, and between colonies managed by different beekeeping operations, which encompasses numerous confounding variables including precise geographic location, intensities of exposures to pathogens, honey bee colony genetics, the availability and quality of bee forage, and specific management practices, all of which may impact colony health and pathogen prevalence and abundance. In this sample cohort, the abundances of BQCV, LSV1, KBV, and L. passim, and the total pathogen abundance differed between beekeeping operations (Table 3). However, there was a large gap in data between March and August for each beekeeping operation, thus this study cannot attribute differences in pathogen abundance to apiary management. In previous US based honey bee monitoring studies, the specific beekeeping operation affected pathogen abundances ; however, this was not observed in colonies monitored as part of the German Bee Monitoring Project . Additionally, in other studies, apiaries affected by CCD had more weak colonies with higher pathogen prevalence within closer proximity to other weak colonies with high pathogen prevalence . Future investigations aimed at better understanding the mechanisms and dynamics of pathogen transmission between colonies within the same apiary (e.g., via mites, foraging/pollen, robbing, beekeeper mediated) are important to developing strategies to reduce inter-colony transmission.
10,22,28,31,64]. To investigate pairwise interactions between the pathogens monitored in this study, we generated a correlation matrix illustrating the correlation coefficients (R-values) that quantify the strength and direction of the relationship between pathogen abundance data (Fig 6). The strongest positive correlation observed in this sample cohort was between V. destructor mites and DWV (R = 0.70, p-value < 0.001). DWV and L. passim (R = -0.65, p-value < 0.001), followed by V. destructor mites and L. passim (R = -0.49, p-value < 0.001), demonstrated the strongest negative correlations (Fig 6).
Fig 6. Pathogen co-infection correlation matrix.
7,10,21,22] and honey colonies located throughout the world [31,61,65–68] is influenced by multiple biotic and abiotic factors. In addition to pathogens, other factors including seasonal variation in colony size, foraging requirements, nutrition, and pesticide exposure all impact colony health [10,35,38,69,70]. Therefore, temporal monitoring projects that emphasize robust sampling designs (e.g., consistent monthly samples obtained on a specific date) and include large sample sizes and multiple monitored factors, are required to fully describe both general and acute threats to honey bee colonies throughout the year and in the context of pollination of specific crops. Additionally, the synergistic effects of co-infecting honey bee pathogens [35,64], which can impact the prevalence and abundance of other pathogens [7,28,31,49], must be investigated from a multidimensional approach to study the effects of pathogens on honey bee colony health. The combination of sampling individual honey bee colonies on specific days of the year and across multiple sites, and employment of statistical methods aimed at integrating and interpreting multifactorial data sets, will further our understanding of the factors with the greatest influence on colony losses and may result in the development of strategies to limit future losses.
S1 Table. Primers used in this study.
S2 Table. Pathogen prevalence by colony health rating.
S3 Table. Complete data set.
S4 Table. Pathogen prevalence by sampling event.
S1 Fig. Analysis of mean pathogen prevalence by colony health rating.
S2 Fig. Analysis of mean pathogen prevalence by sampling event.
S3 Fig. Changes in total pathogen abundance in honey bee colonies throughout the sampling period.
S4 Fig. Changes in KBV abundance in honey bee colonies throughout the sampling period.
S5 Fig. Changes in LSV2 abundance in honey bee colonies throughout the sampling period.
S6 Fig. Changes in LSV1 abundance in honey bee colonies throughout the sampling period.
S7 Fig. Changes in N. ceranae abundance in honey bee colonies throughout the sampling period.
S8 Fig. Changes in L. passim abundance in honey bee colonies throughout the sampling period.
S9 Fig. Changes in BQCV abundance in honey bee colonies throughout the sampling period.
S10 Fig. Correlations between changes in frame count in response to individual pathogen abundances after accounting for sampling date.