Additional details regarding methods are available in the Supplemental Information (Online Resource 1).

Longitudinal monitoring and sampling of commercially managed honey bee colonies

Three Montana-based (Broadwater, Yellowstone, and Treasure counties) commercial beekeeping operations that transport their honey bee colonies ~1,200 miles to California (Merced and Stanislaus counties) each winter for the almond bloom provided honey bee samples before (October–December 2013), during (February 2014), and after (March/April and June 2014) almond pollination (Figure 1 and Supplemental Table S2). Colony health, using colony population size as a proxy, was assessed by the number of frames covered with honey bees (frame counts) at each sampling event (Delaplane and van der Steen 2013; OSU 2011). Colony strength was defined as follows: weak colonies (<5 frames covered with bees), average colonies (6–8 frames covered with bees), and strong colonies (>9 frames covered with bees). Live honey bee samples (~100 per sample) were obtained from the top of the frames in the middle of the colony. Samples were composed of female bees of mixed age, including nurse, worker, and forager bees. The samples were collected on ice or dry ice, stored at −20 °C, shipped on dry ice, and transferred to −80 °C prior to analysis. At the onset of the study in November 2013, each beekeeper identified 15–20 colonies of differential health. Specifically, operation 3 initiated the study with five weak, five average, and five strong colonies and provided samples at three time points; operation 2 initiated the study with five weak, 13 average, and two strong colonies and provided samples at four time points; and operation 1 initiated the study with five weak, four average, and ten strong colonies and provided samples at four time points (Supplemental Table S2). A total of 176 honey bee samples with corresponding colony strength observations were obtained and analyzed, four observations of colony strength lacked corresponding samples, and eight of the original colonies died during the course of this study (Supplemental Table S2).

Figure 1. Longitudinal monitoring of commercially managed honey bee colonies, before, during, and after the 2014 almond pollination season. Honey bee colonies (n = 54) from three Montana-based commercial beekeeping operations were monitored before (October–November 2013), during (February 2014), and after (March–April); after 2 (June) almond pollination in California. Colony strength was measured at each sampling event. PCR was utilized to detect pathogens associated with each sample, and qPCR was utilized to determine the abundance of pathogens associated with a subset of samples (Supplemental Table S2). Full size image

Honey bee samples

Five female bees from each sample were used for RNA extraction, cDNA synthesis, pathogen-specific PCR, and qPCR. The objective for pathogen screening was to identify the most prevalent pathogens associated with honey bees sampled from individual colonies at each sampling event. Based on empirical data, literature values, and practical sample handling considerations, we assayed five bees per colony per sampling event. 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 (Pirk et al. 2013); this sample size has been proven sufficient for the pathogen-specific PCR detection of highly prevalent pathogens (Daughenbaugh et al. 2015; Runckel et al. 2011).

RNA isolation

Bee samples were homogenized in water using beads (3 mm) and a TissueLyzer (Qiagen) at 30 Hz for 2 min. Samples were centrifuged for 12 min at 12,000×g at 4 °C to pellet debris, and RNA from supernatants was extracted using TRIzol reagent (Life Technologies) according to the manufacturer’s instructions (Runckel et al. 2011).

Reverse transcription/cDNA synthesis

cDNA synthesis reactions were performed by incubating 1,000–2,000 ng total RNA, Moloney murine leukemia virus (M-MLV) reverse transcriptase (Promega), and 500 ng random hexamer primers (IDT) for 1 h at 37 °C, according to the manufacturer’s instructions (Runckel et al. 2011).

Polymerase chain reaction (PCR)

PCR was performed according to standard methods using the primers listed in Supplemental Table S1 (Runckel et al. 2011). In brief, 1 μL cDNA template was combined with 10 pmol of each forward and reverse primer and amplified with ChoiceTaq polymerase (Denville) according to the manufacturer’s instructions using the following cycling conditions: 95 °C for 5 min; 35 cycles of 95 °C for 30 s, 57 °C for 30 s, and 72 °C for 30 s, followed by final elongation at 72 °C for 4 min. The PCR products were visualized by gel electrophoresis/fluorescence imaging.

Quantitative PCR (qPCR)

Quantitative PCR was used to analyze the relative abundance of the most prevalent pathogens, which were all RNA viruses, in select samples to investigate the relationship between virus abundance and honey bee colony health. Five hundred nanograms of RNA from each of these samples was reverse transcribed with M-MLV as described above. All qPCR reactions were performed in triplicate with a CFX Connect Real Time instrument (BioRad); reaction conditions and equations for determining the relative abundance based on standard curves are provided in supplemental methods (Online Resource 1).

Statistical analysis of PCR

For this study, we use “pathogen prevalence” to refer to the total number of pathogens detected by PCR out of a target list of 16. Though our interest was in the relationship between strength rating and pathogen prevalence, graphical analyses indicated that there were likely relationships between pathogen prevalence and sampling time as well as between strength and sampling time. Thus, we used a Poisson log-linear regression model and accounted for an interaction between sample date (time period), beekeeping operation, colony strength, and pathogen prevalence. Observations with average strength rating were not included in some analyses to simplify the inferences between strong (S) and weak (W). The natural logarithm (ln) of the pathogen prevalence data was used in comparisons between each beekeeping operation and time period combination (Pirk et al. 2013). For the model, we used beekeeping operation 1, before almond pollination (time period 1), and weak colonies as the base level.

In all, our model can be expressed

$$ \begin{array}{c}\hfill {y}_i\sim \mathrm{Poisson}\left({\mu}_i\right)\hfill \\ {}\hfill \log \left({\mu}_i\right)={\beta}_0+{\beta}_1\times \mathrm{operation}\ {2}_i+\hfill \\ {}\hfill {\beta}_2\times \mathrm{operation}\ {3}_i+{\beta}_3\times {\left(S:\mathrm{period}\ 1\right)}_i+\hfill \\ {}\hfill {\beta}_4\times {\left(\mathrm{W}:\mathrm{period}\ 2\right)}_i+{\beta}_5\times {\left(\mathrm{S}:\mathrm{period}\ 2\right)}_i+\hfill \\ {}\hfill {\beta}_6\times {\left(\mathrm{W}:\mathrm{period}\ 3\right)}_i+{\beta}_7\times {\left(\mathrm{S}:\mathrm{period}\ 3\right)}_i+\hfill \\ {}\hfill {\beta}_8\times {\left(\mathrm{W}:\mathrm{period}\ 4\right)}_i+{\beta}_9\times {\left(\mathrm{S}:\mathrm{period}\ 4\right)}_i\hfill \end{array} $$

where y i = the total abundance/prevalence for the i th observation i = 1, 2, …, 180.

Operation2 i = 1 if observation i came from beekeeping operation 2 and 0 otherwise.

Operation3 i = 1 if observation i came from beekeeping operation 3 and 0 otherwise.

Period 2 i = 1 if observation i was taken during and 0 otherwise.

Period 3 i = 1 if observation i was taken after pollination and 0 otherwise.

Period 4 i = 1 if observation i was taken in the second after pollination sampling time and 0 otherwise.

A i = 1 if observation i was average (colony strength) and 0 otherwise.

S i = 1 if observation i was strong (colony strength) and 0 otherwise.

Statistical analysis of qPCR

The relationship between colony strength rating and pathogen abundance was evaluated using a log-normal regression to model the total abundance with the predictor of interest, strength rating, while also accounting for the different beekeeping operations and the sampling time period. To evaluate if the relative abundance of the most prevalent pathogens, which were all (+) sense RNA viruses (i.e., BQCV, SBV, LSV1, and LSV2), we utilized virus copy number as an indication of infection severity, though the relationship between virus copy number and virus associated disease or affects on the honey bee host are largely unknown (de Miranda et al. 2013). Pathogen abundance was defined as the summed abundance of the RNA virus genome copy numbers, which were measured by qPCR. Samples that did not test positive for the virus by PCR were not assessed by qPCR and given a value of zero. In total, there were 53 observations of total abundance after inputting zeros for negative PCR tests. A log-linear regression was used to model the total abundance with the predictor of interest, strength rating, while also accounting for the different beekeeping operations and the sampling time period; 1 was added to each observation since some observations had 0 total abundance. Some bee colonies were measured multiple times in the 53 observations. We accounted for these repeated measures on colonies with a random effect for colony, but found the variance between colonies to be minimal compared to the overall variance.

$$ \begin{array}{c}\hfill \log \left({y}_i\right)={\beta}_0+{\beta}_1\times \mathrm{operation}\ {2}_i+\hfill \\ {}\hfill {\beta}_2\times \mathrm{operation}\ {3}_i+{\beta}_3\times \mathrm{period}\ {2}_i+\hfill \\ {}\hfill {\beta}_4\times \mathrm{period}\ {3}_i+{\beta}_5\times \mathrm{period}\ {4}_i+\hfill \\ {}\hfill {\beta}_6\times {S}_i+{\beta}_7\times {A}_i+{\gamma}_{j(i)}+{\epsilon}_i\hfill \end{array} $$