A derived honey bee stock confers resistance to Varroa destructor and associated viral transmission

Colonies

Colony setup occurred prior to initiation of the study, between March and May 2017, in Mississippi, USA. Using established methods, queenless colony divisions, obtained from a large commercial beekeeping operation, were equalised to an average calculated population size of ~ 7000 workers¹¹², and housed in 10-frame Langstroth hives (Table S1). After acclimatisation for 24–48 h, they each received an imminently emerging queen cell, containing a queen from one of two stocks, added to the same worker baseline. The stocks used consisted of an Italian ‘Commercial’ stock, propagated from collaborator established breeder queens, and thus representative of the industry standard, and the Varroa-resistant ‘Pol-line’ stock⁵⁴. To ensure consistency, all queens were reared in the same ‘cell builder’ colonies, based at the USDA Honey Bee Breeding, Genetics and Physiology Laboratory, in Baton Rouge, Louisiana, USA. Colonies from each stock were held in independent apiaries, 80 km apart to maintain physical isolation; and to control genetic fidelity, virgin queens were open mated to drones of the same stock via drone saturation. Fourteen days after queen emergence, colonies were inspected, and mated queens were marked with paint on the thorax, to assist with identification, with white corresponding to Commercial, and blue to Pol-line. Colonies were allowed to acclimatise for six weeks before sampling began, and those that failed to achieve mating success, or had unacceptably high [≥ 3.0 ‘mites per hundred bees’ (MPHB)] Varroa levels, were removed, normalising the average between-stock Varroa difference to < 0.1 MPHB from the initiation of the study (Table S1). Each colony was then assigned a random ID number, for a total of 366 colonies; 193 with Commercial queens, and 173 with Pol-line queens. All colonies were provided with supplemental sucrose solution, and soya-based commercial pollen supplements ad libitum, in accordance with standard industry practices.

Experimental design

In order to evaluate the effects of migratory pollination procedures, colonies were randomly assigned to one of three migration route treatments. The migration treatments used were as follows; a ‘California’ group, consisting of 156 colonies; 80 Commercial and 76 Pol-line, moved to South Dakota in May, followed by California in October for overwintering and almond pollination, and back to Mississippi at the end of the study; a ‘Mississippi’ group, consisting of 156 colonies; 80 Commercial and 76 Pol-line, moved to South Dakota in May, and back to Mississippi via California in October for overwintering; and a ‘Stationary’ group, consisting of 54 colonies; 32 Commercial and 22 Pol-line, that remained in Mississippi for the duration of the study (Fig. 1). While in South Dakota, colonies were distributed across four apiaries, and in California and Mississippi, they were held in a single apiary. This practice mirrored standard industry protocol, and provided a more representative sample by accounting for localised differences in climate and forage.

In addition to migration route, colonies were divided by the frequency of acaricide treatment applied, to better elucidate the effects of varying infestation strata. Within each migration treatment, and stock, half of the colonies received a ‘high’ acaricide treatment, and half a ‘low’ acaricide treatment. Colonies in the high treatment group were treated twice, in September and December; while those in the low treatment group were treated only once, in December. All treatments were conducted using the acaricide amitraz.

Sampling

Colonies were sampled five times during the course of the study, in May, June, September, December, and February. For each colony, the queen was examined, and frames of workers, nectar, pollen, and brood were quantified visually by percentage cover¹¹³. Upon observation, if the queen did not have an existing mark, a single red paint mark was applied to her thorax, to indicate supersedure status. Then, ensuring that the queen was sequestered to prevent unintentional sampling, two samples of ~ 300 workers were taken in resealable bags; one being flash frozen with dry ice, and stored at − 80 °C for pathogen analyses; and the other being frozen with conventional ice, for Varroa analysis. To achieve a representative age sample, workers were removed from two brood frames, containing sealed and unsealed brood. Workers were evenly mixed using a 20 L bucket, before collection via a cup-measure, and subsequent assignation to one of the two sample types. The queen was then returned to the hive, along with any unsampled workers.

Varroa analysis

Varroa analyses were conducted for all colonies, at all five time points. To quantify the level of Varroa infestation in colonies, a detergent-wash method was employed¹¹⁴. The sample of bees was first placed in a mesh-partitioned cup with a solid base, to allow mite retrieval. Water was then added until the bees were submerged, and a surfactant-based detergent (Procter & Gamble) mixed in, ensuring an even coating. The lid of the cup was secured, and it was placed into a Model E5850 reciprocal shaker (Eberbach), to be agitated for a period of 60 min at 120 rpm. After agitation, the lower faction of the cup was removed, and the contents poured into a shallow tray, to count the number of mites falling through the mesh. This process was repeated until two consecutive zero counts were recorded, indicating that all mites had been removed. The bees were then separated and counted by tally counter, to produce an exact sample size. Finally, the total Varroa count was transformed to produce a standardised ‘mites per hundred bees’ (MPHB) value, as follows:

Pathogen analyses

Four specific viral targets were chosen: DWV-A, DWV-B, CBPV, and BQCV. Analyses were carried out on a randomised subset of 92 colonies from the California and Mississippi migration groups, divided evenly across stock, and survival outcome. Four major time points were selected: June, September, December, and February. All analyses were conducted using RT-qPCR, following previously established methods^94,115.

RNA extraction

Per colony and time point, pools of 65 bees were randomly selected from the master samples of ~ 300, placed into 30 ml 19-6358Z bead tubes (Omni), and stored at − 80 °C, for homogenisation using a Bead Ruptor Elite (Omni). Immediately after homogenisation, 5 ml of homogenisation solution (Promega), and 4 ml sterile 1xPBS buffer, each at 4 °C, were added to sample tubes, before vortexing for 20 s using a Vortex-Genie 2 (Scientific Industries). 1.8 ml of each sample was then transferred to 2 ml tubes (Eppendorf), and centrifuged for 60 s at 5000 rpm, at 4 °C, using a 5430-R centrifuge (Eppendorf). Per sample, 400 μl of the resultant supernatant was then transferred to Maxwell RSC-48 cartridges (Promega), prepared in accordance with the manufacturer’s instructions, and RNA extractions were carried out using the Maxwell RSC simplyRNA tissue protocol (Promega) (Maxwell RSC simplyRNA Tissue Kit, AS1340, TM416, 2019).

Sample purity and yield were then calculated in duplicate, using a NanoDrop One microvolume UV–Vis spectrophotometer (Thermo Fisher Scientific), and appropriate RNA dilutions were completed to ensure sample standardisation to 100 ng/μl. Following standardisation, cDNA synthesis was achieved using QuantiTect Reverse Transcription Kits (Qiagen), utilising a T100 thermal cycler (Bio-Rad), running custom QuantiTect cDNA synthesis protocols (Qiagen) (QuantiTect Reverse Transcription Handbook, 2009). Resultant samples were stored at − 20 °C, prior to initiation of RT-qPCR analyses.

RT-qPCR

RT-qPCR was performed on 1 μl aliquots of each sample, in triplicate, in a total reaction volume of 10 μl, utilising SsoAdvanced Universal SYBR Green Supermix (Bio-Rad), on Multiplate 96-well optical PCR plates (Bio-Rad). The primers used for quantification of DWV-A¹¹⁶, DWV-B¹¹⁷, CBPV⁹⁶, and BQCV¹¹⁵ have been reported previously, and validated for target specificity. For full sequence details, see (Table S2). Negative controls were included for each target, consisting of 1 μl of nuclease-free H₂O (Promega), again run in triplicate. Additionally, to enable viral quantity mean determination, standard curves were produced for all targets, via tenfold serial dilutions of known viral quantities, covering 8 orders of magnitude. Linearity (r²), and reaction efficiency (E), were maintained at ≥ 0.990, and ≥ 92.5% respectively, for all assays (Table S2). All analyses were run on CFX Connect Real-Time PCR Detection Systems (Bio-Rad), using previously optimised thermal protocols (Table S3).

Data transformation

Samples were analysed in triplicate, to form a mean Ct value. Any triplicates with a standard deviation ≥ 1.0 were examined, and the divergent replicate removed, resulting in 26 technical replicate removals. If this failed to bring the standard deviation to < 1.0, the sample itself was replaced, and removed from further analyses. Ct means were then quantified against the standard curves for each target, to calculate absolute viral titres, using the following equation:

Due to the wide distribution of quantities inherent to viral replication dynamics¹¹⁸, these values were then log-transformed, to produce suitable data for further analyses. Where undetected values were present, a constant Ct mean of 45 was assumed, to allow for subsequent transformation of the data¹⁰⁵.

Range of assessment factors

Analyses were broadly divided into measurements of colony survival, the extraneous influences of Varroa levels (MPHB), and viral titres (log₁₀ viral quantity mean per μl cDNA) upon them, and how stock, mite treatment, and migration route modulated these interactions. Crucially, we examined both the general effect of Varroa levels and viral titres upon colony survival, and the predictive power of such measures at seasonally-relevant time points. We additionally determined the explanatory power of viral titres when colonies were matched for Varroa level, thus elucidating the relative additive influences of Varroa-transmissible viruses.

Colony and queen survival over time, population size, and honey production

First, the survival of colonies and queens over time, along with the population sizes of surviving colonies, and honey production, were compared between the two stocks using time series data collected throughout the course of the study. Colonies were classified as ‘dead’ when the worker frame-count dropped to < 1.0, and they were confirmed as vacant during any sequential sampling. This definition of colony death remained consistent across all further analyses. Queens were classified as ‘dead’ when the existing queen could not be found, and a new unmarked queen, or queen cells, were present, indicating supersedure¹¹⁹. Population sizes were compared using worker frame-counts from all colonies classified as ‘alive’, that is, those with ≥ 1.0 frames at the end of the study, in February. Honey production was quantified via net weight (kg) of honey extracted per colony during September, excluding colonies in the stationary migration group, for which extraction was not recorded.

Factors influencing colony survival

We then defined the influences of Varroa levels and viral titres upon colony survival, using data from all five and four time points, respectively. The effect of each factor, along with stock, mite treatment, migration route, and their concordant interactions, was quantified. By tracking individual colonies over time, in concert with their Varroa levels, and viral titres, we were able to determine factor effect sizes, and the interactions leading to death or survival in February. February was chosen as the time point of interest, as it fell immediately prior to almond pollination, and thus gave a representative, and commercially relevant, measure of colony strength going into this system. As viral titre measures were taken for a subset of colonies, separate analyses were conducted for these and Varroa measures, in order to preserve maximum sample sizes in each case. These data were then used to inform subsequent predictive analyses, based on Varroa levels and viral titres.

Predictive differences in Varroa levels and viral titres

To assess initial predictive differences in Varroa levels and viral titres, across survival outcomes and stocks, we linked survival status in February to two predictive time points: June, and September. Notably, the former of these encompassed the principle period of colony population growth, and the latter, of honey harvesting, prior to colony overwintering. These points constitute key junctures at which management decisions must be made, and thus are practically relevant as predictive benchmarks. In analyses, colonies were classified by stock, or their status as being either alive or dead in February. Varroa levels and viral titres were then compared between groups, thus representing a basic test of the usefulness of single time point data in defining outcomes.

Prognostic power of Varroa levels and viral titres

To quantify the epidemiological significance of Varroa levels and viral titres in determining colony survival, and therefore their predictive power, we conducted relative risk (RR) analyses for these factors, utilising colony death by February as the outcome of choice. RR is an established epidemiological measure, used to determine the probability of a given health outcome when exposed to a specific risk factor, and thus compare risk magnitudes¹²⁰. RR analyses were conducted for both Varroa levels, and viral titres, again for the predictive time points of June and September. Significant interactions were further validated with attributable risk (AR) analyses, to provide a measure of effect size; and complementary Bayesian RR analyses, to visualise the risk probability distribution. Analyses were pooled across stocks, in order to better encompass the full range of infestation and infection strata. Additionally, they were stratified by Varroa level, to determine how changes in infestation severity influenced the risk of colony death.

We then isolated the additive effect of viral titres as a predictive factor, independent of Varroa. We considered this an important test of the stand-alone influence of Varroa-associated pathogens in Varroa-mediated colony loss, as ordinarily, these factors are difficult to decouple. To achieve this, a subset of 60 colonies were matched by September Varroa level (mean MPHB, < ± 0.1 between survival outcomes), and left unstratified in terms of viral titres. September was chosen as the key matching point for this analysis, as colonies showed the highest Varroa levels and viral titres at this stage, and thus these values were considered most relevant to overwinter survival outcomes. Colonies were subsequently divided into alive and dead cohorts, based on their status in February, to produce groups that differed in their survival outcomes, and viral titres, but that had tightly coupled Varroa level distributions. We then conducted RR analyses for this subset, examining viral titres in both June and September as predictive factors for colony survival, hence constituting a measure of viral effects, decoupled from Varroa.

Varroa model

To demonstrate the relationship between Varroa level and colony mortality, and generate a predictive model for treatment prognoses in both stocks, we pooled all colony data into Varroa infestation strata, derived from Varroa level in September, and assigned resultant percentage mortalities based on survival rates in each stratum. We then fit curve models to the data for the two stocks, thus enabling a predictive mortality outcome for the Varroa level and stock in question, based on an extensive empirical dataset.

Statistical analyses

For all pairwise comparison data in which analyses would assume a normal distribution, we performed Shapiro–Wilk tests to check for normality, and hence inform the application of appropriate statistical tests. In all cases, the data were not normally distributed, and thus independent sample Mann–Whitney U-tests were employed. For these analyses, we opted to report mean ranks, rather than medians, as based on descriptive statistics, we did not assume identical frequency distributions between sample groups. Due to the large sample sizes used, we also considered it important to include a measure of effect size, and this was achieved via Eta-squared (η²) post hoc tests¹²¹. Additionally, for analyses that were close to the threshold of significance (0.05–0.06), we cross-verified results by applying parametric methods to transformed data. This occurred in one case, and after transformation via reflection of the square root, the significance of the test was unaltered.

Cumulative colony and queen survival differences between stocks were assessed using mixed-effects Cox proportional hazards models, via the R packages ‘coxme’¹²² and ‘multcomp’¹²³. The proportional hazards assumption was verified in each case by confirming independence between the Schoenfeld residuals and time, for both the fixed factors and model, using the R package ‘survival’¹²⁴.

To assess differences between stocks in colony population size and honey production, we utilised generalised linear models (GLMs). The population size model used February as the response time point, while the honey production model used September. Model selection was based on AIC, and validated using omnibus tests, along with assessment of the deviance to degrees of freedom ratios.

To determine how Varroa levels and viral titres influenced colony survival outcome, and how stock, mite treatment, and migration route modulated these interactions, we employed generalised linear mixed models (GLMMs) with repeated measures. Separate models were generated for Varroa level and viral titre analyses, in order to maximise sample sizes, and minimise model complexity. Varroa level measurements were repeated over five time points, and viral titre measurements over four. Model selection was based on AIC, beginning with the full model and interactions. In all cases, model fit was validated via evaluation of the binned standardised residuals.

When testing the predictive power of Varroa levels and viral titres at single time points, we used relative risk (RR) analyses. All RR analyses made use of an ‘outgroup’ for the epidemiological factor in question, and then compared the mortality rate in this group to that of an ‘exposure’ group, to ascertain the relative change in risk. The exposure group in all viral titre analyses was defined as colonies with titres ≥ 1 × 10⁷, a cut-off value indicative of an epidemiologically severe infection^1,96,98, while the outgroup was composed of colonies with titres < 1 × 10⁷. The exposure groups in the Varroa level analyses were stratified to cover the range of MPHB levels present in colonies, ≥ 1–2.5, > 2.5–5, > 5–7.5, and > 7.5, and the outgroup was maintained at < 1, thus providing a graded measure of the change in risk with increasingly severe infestations. RR values greater than one indicated an associated increase in mortality risk, while those less than one indicated an associated reduction. In accordance with standard interpretation, RR values were deemed not significant if the resultant 95% confidence intervals overlapped the value of one¹²⁰. For significant RR values, we also calculated the AR, to demonstrate the magnitude of absolute effect size, as is recommended for binary epidemiological outcomes¹²⁵. Significant risk probability distributions were then visualised via parallel Bayesian RR analyses, using the R package ‘brr’¹²⁶. Across tests, we confirmed requisite sample sizes to provide a minimum power (1 − β) of 0.80, at an alpha (α) of 0.05, using standard deviation (σ) and mean difference (δ) values from the data. All statistical analyses were performed in PASS (release v. 21.0.2), SPSS (release v. 25.0.0.0), and R (release v. 4.0.2)¹²⁷.

Colony and queen survival over time

We used mixed-effects Cox proportional hazards models to assess differences in colony and queen survival between the two stocks over the course of the study. These utilised stock, mite treatment, and migration route as fixed factor predictors, and colony ID as a random factor. As sampling occurred every two months, cumulative survival was updated on a bi-monthly basis, over a 10 month period.

Population size and honey production

The GLM assessing the effect of stock on colony population size used frame-count of surviving colonies in February as a gamma response variable with a log link, and stock, mite treatment, and migration route as fixed factor predictors. The GLM assessing the effect of stock on honey production used net honey weight in September as a gamma response variable with a log link, stock and mite treatment as fixed factor predictors, and omitted migration route as honey was extracted prior to transport.

Factors influencing colony survival

The GLMM assessing the effect of Varroa levels on colony survival used colony survival status in February as a binomial response variable with a probit link, Varroa MPHB, stock, mite treatment, and migration route as fixed factor predictors, their two-way interactions, and colony ID as a random factor. The GLMM assessing the effect of viral titres on colony survival used colony survival status in February as a binomial response variable with a logit link, log₁₀ DWV-A titre, log₁₀ DWV-B titre, log₁₀ CBPV titre, log₁₀ BQCV titre, stock, mite treatment, and migration route as fixed factor predictors, their two-way interactions, and colony ID as a random factor.

Predictive differences in Varroa levels and viral titres

Independent sample Mann–Whitney U-tests were utilised to compare differences in mean Varroa levels and viral titres, based on colony survival outcome and stock, for June and September time points.

Prognostic power of Varroa levels and viral titres

To quantify the relative predictive powers of Varroa levels and viral titres in June and September when determining colony mortality, and that of viral titres when matched for Varroa level, we used RR, and AR, analyses. To visualise risk probability distributions, we used Bayesian RR analyses, utilising a Poisson response, and the Jeffreys prior as an objective reference prior.

Varroa model

A curve estimation analysis was employed to fit curve equations to the Varroa model data, and curve selection was evaluated via the resultant model outputs, based on F-values.

Source: Ecology - nature.com