All colonies were started from 2 to 2.5 lb. “packages” made on 3 May 2018 from 10 previously established Italian colonies at the USDA Honey Bee Breeding, Genetics, and Physiology Research Laboratory, Baton Rouge, LA (30°22′56″N, 91°10′40″W). Naturally-mated queens (n = 3 per stock) from the five genetic stocks were sourced from the USDA Bee Lab [Pol-Line and Russian, Saelao et al. (41)], a Canadian collaborating breeder (Saskatraz), or purchased from commercial suppliers (Carniolan and Italian). To increase worker acceptance of new queens, all queens were placed into new colonies on 4 May 2018 using plastic cages with a candy plug blocking the entrance that the queen and workers would chew through over a period of days for final release (52). Colonies were checked for the presence of the queen and the number of brood frames; those with fewer than three frames of bees were supplemented with additional brood frames on 17 May 2018 to equalize starting populations. All colonies were maintained following standard management practices in three apiaries within 6 km of each other (with Carniolan, Italian, and Saskatraz sharing an apiary, while Pol-Line and Russian colonies were kept in two separate apiaries). Colonies were not sampled until after 6 weeks post brood supplementation to allow time for population turnover to reflect queen genetics.

To obtain the DWV viral solution, 20 adult bees with DWV infection symptoms were collected and ground in liquid nitrogen to a fine powder, homogenized in 10 ml PBS, and centrifuged at 5,000 rpm for at 4°C for 20 min following established protocols (29, 53, 54). The resulting supernatant containing viruses was filtered through a 0.2-micron filter (milex-GS syringe filter unit #SLGS033SS, Millipore Sigma, Burlington, MA, USA) to remove small tissue debris, fungi, and bacteria. qPCR was conducted to test for the presence of non-target viruses (Acute Bee Paralysis Virus, Black Queen Cell Virus, Chronic Bee Paralysis Virus, Israeli Acute Paralysis Virus, Kashmir Bee Virus, and Lake Sinai Virus) using the methods described below (primers in Supplementary Table 1). Viral quantification for non-specific DWV was performed by absolute quantification using the Standard Curve Method. All methods were previously established based on standard protocols (27, 34). One sample stock solution for DWV (measured using non-specific DWV primers) was selected based on negative results for non-target viruses and used to create the injection stock solution. Stock solutions were diluted to 10⁵, a biologically relevant, sublethal functional titer level for adult bees (1).

From July to October 2018, frames with emerging adult bees were brought into the lab where emerging bees from each colony were immediately uncapped and removed from the frame. Bees were inspected for Varroa mites and those with mites on them or in their pupal cells were discarded to try to mitigate the impacts of prior feeding and DWV transmission. To simulate the vectoring of DWV through the feeding of Varroa mites in a standardized way, bees were injected with 3.0 μL of DWV inoculum (DWV treatment). Three microliters of 1X PBS injection (PBS treatment) and no injection (included as bees may have had naturally occurring DWV infections, hereafter referred to as “control”) were implemented as controls. To reduce movement during injection, bees (including controls) were placed on ice for 2 min then injected (PBS and DWV treatments) using an UltraMicroPump with a SYS-Micro4 Controller (World Precision Instruments, Sarasota, FL, USA) with an infusion flow rate of 0.1 μL/s, following manufacturer's parameters. For the injection, a 30G needle (Hamilton Company, Reno, NV, USA) was inserted into the lateral abdomen between the fourth and fifth tergites, based on established protocols (54, 55). Following injections, bees were housed in cages (maximum of n = 30 bees) of the same treatment and colony according to standard methods (56). All cages were provided 50% sucrose solution and pollen substitute to ensure hypopharyngeal gland development; food was replenished as needed or when desiccated (56, 57). Cages were maintained in an incubator at 34°C and 85% relative humidity. A subset of three bees from each of the three treatments and 15 colonies were sacrificed at 1, 2-, 4-, 7-, and 10-days post-injection for a total of 675 bees and stored in sterile 1.5 mL centrifuge tube at −80°C. Given inherent mortality differences, more bees were treated as needed to obtain the necessarily time point samples.

To determine virus movement within the body of honey bees over time, the three-bee subset was dissected from each colony/treatment/stock/time point combination. Dissections were conducted over dry ice with each bee dissected with a new, sterilized blade. For each bee, the body was separated into legs, head, and abdomen with other body segments removed. The head was then embedded into beeswax (replaced for each dissection) and the hypopharyngeal gland removed according to previously published methods (58). Dissected tissues were stored in separate sterilized tubes on dry ice during dissection and long-term at −80°C until RNA extraction.

RNA was extracted for a single rear leg, the head (sans hypopharyngeal gland), the hypopharyngeal gland, and the abdomen for each of the three bees representing each combination of colony/treatment/stock/time (n = 2,700 RNA extractions in total). To extract virus RNA from hypopharyngeal glands, samples were placed in 30 μL lysis buffer and 30 μL Maxwell Homogenization buffer and vortexed. For legs (cut into pieces), heads, and abdomens, samples were placed in 200 μL lysis buffer and 200 μL Maxwell homogenization buffer, manually ground with a pestle (Sigma-Aldrich, St. Louis, Missouri, USA), and vortexed. All samples were then incubated for 90 min at 4°C. After incubation, 320 μL Maxwell Homogenization buffer was added to hypopharyngeal gland samples. Total RNA of each tissue sample was extracted using the Maxwell RSC 48 cartridges (Promega Corporation, Madison, Wisconsin, USA) according to standard procedures of Maxwell RSC simplyRNA tissue extraction kits and program (Promega Corporation, Madison, Wisconsin, USA). RNA was quantified via NanoDrop One (Thermo-Fisher Scientific Inc., Waltham, Massachusetts, USA) twice using 1 μL of sample each time. The mean ng/μL NanoDrop One readings were calculated per sample then used to determine the quantities of sample and nuclease-free water required to reach a sample concentration of 250 ng of RNA. RNA was stored in 0.6 mL elution tubes wrapped in parafilm (Bemis NA, Neenah, Wisconsin, USA) at −80°C until cDNA synthesis (59, 60).

Frozen RNA samples were thawed on −20°C metal beads, briefly vortexed, then centrifuged. cDNA was then synthesized in two steps using Qiagen QuantiTect Reverse Transcription kits (Thermo-Fisher Scientific Inc., Waltham, Massachusetts, USA). For step one, 2 μL of gDNA wipeout solution was added to the mix of RNA and water for a total reaction volume of 14 μL per sample. Samples were incubated at 42°C for 2 min in a Bio-Rad T100 Thermal Cycler (Bio-Rad, Hercules, California, USA). Samples were briefly vortexed and centrifuged before the addition of 4 μL 5X Buffer, 1 μL of RT Primer mix, and 1 μL of RT enzyme per sample. Samples were again briefly vortexed and centrifuged then placed into the Bio-Rad T100 Thermal Cycler (42°C for 25 min then 95°C for 3 min). Synthesized cDNA was stored in strips tubes wrapped in parafilm at −80°C until RT-qPCR.

To quantify DWV-A and DWV-B levels (1, 22, 61, 62), each sample was replicated two times per primer pair for RT-qPCR analyses. All RT-qPCRs consisted of 5 μL SsoFast Universal SYBR Green supermix (Bio-Rad, Hercules, California, USA), 3 μL nuclease-free water, 0.5 μL forward primer, 0.5 μL reverse primer, and 1 μL cDNA from the sample. All reactions were run in Bio-Rad CFC 96 or Connect Thermal Cyclers (Bio-Rad, Hercules, California, USA) with all reactions of a specific primer occurring in the same machine. All samples were tested with DWV-A and DWV-B primers to determine infection levels (primers in Supplementary Table 1). The PCR cycling protocol for DWV-A was 95°C for 1 min followed by 40 cycles of 95°C for 10 s and 60°C for 15 s then 65°C for 5 s; while the protocol for DWV-B was 95°C for 5 min followed by 40 cycles of 95°C for 5 s and 52.5°C for 10 s then 72°C for 10 s. The thermal protocols included a melt-curve dissociation analysis to confirm product size. DWV-A and DWV-B results were quantified using the Standard Curve Method using linearized plasmid constructs. Quantified virus titer levels were log-transformed for analyses.

We analyzed DWV levels in four tissue types per bee for three individual bees per timepoint (n = 5) and treatment (n = 3). This was replicated simultaneously for three colonies for each of five honey bee stocks for a total of 2,700 RNA extractions and 675 individual bees (Figure 1). The factors influencing levels of DWV-A and DWV-B were conducted using two general linear mixed models (one per virus genotype) or GLMMs (lme4 package) in R v 3.6.1 (63, 64). For all models, colony and bee ID were used as random effects to account for variation among colonies and individual bees (for tissue type) that were unaccounted for by stock. The following variables were considered for the fixed effects: treatment (no injection control, PBS injection, and DWV injection), tissue type (abdomen, head, hypopharyngeal gland, and legs), genetic stock (Carniolan, Italian, Pol-Line, Russian, and Saskatraz), time since treatment (1, 2, 4, 7, and 10 days), Log DWV (genotype B for the genotype A model and vice versa), all double and triple interaction effects of the previous variables, and a quadratic time term. The no injection control treatment (given that there was naturally occurring DWV infection present), the abdominal tissue (site of Varroa feeding), and the Italian bee stock (commercial standard) were specified as the intercept values; model results indicated in tables are relative to these values. Model fit was evaluated using estimated AICc and BIC scores using the ANOVA function in lmer; those with significantly lower scores (ΔAICc > 4) were used. When the scores were not different, the model containing more variables was used for ease of comparison between DWV-A and DWV-B models. P-values were estimated for all models using Satterthwaite's method in the lmerTest package (65). Post-hoc Tukey tests with a Sidak correction were conducted using emmeans package (66). All graphs were made using ggplot2 (67).

In both virus genotype models, the level of the alternate virus genotype was significant. Note, strains were not further characterized and these data are based on the primers (Supplementary Table 1). In the DWV-A model, both the model (Table 1) and fixed effects (Supplementary Table 2) indicated that greater levels of DWV-B were positively associated with an increase in DWV-A levels (Figure 2). The same occurred in the DWV-B model when DWV-A was included. In general (all treatments and tissue types combined), DWV-B levels were higher than DWV-A (ANOVA, DF = 1, 5,230, F = 1,125.482, and p < 0.001). Furthermore, DWV-B levels were higher than DWV-A when split out by day 1 (ANOVA, DF = 1, 1,037, F = 105.883, and p < 0.001) and day 10 (ANOVA, DF = 1, 1,064, F = 246.818, and p < 0.001). The higher DWV-B levels compared to DWV-A for all treatments are potentially due to higher levels of naturally occurring DWV-B infections. Bees (based on control treatments) started with naturally occurring DWV infections with an average of 10^6.08 for DWV-A and 10^6.92 for DWV-B, which continued to develop throughout the observation period.

There was a strong treatment effect for both DWV-A and DWV-B (Table 1). Of the treatments, only the DWV injection treatment was significantly related to DWV-A levels in the fixed effects (Supplementary Table 2). Control and PBS treatments were similar to each other, and both had significantly lower DWV-A levels than the DWV treatment due to naturally occurring DWV infections (Table 2). There was a slightly different pattern for the DWV-B model—both PBS and DWV injection treatments were significant predictors of DWV-B copy number (Supplementary Table 2). Using Tukey comparisons, PBS had significantly greater levels of DWV-B than did the control but still lower levels than the DWV treatment (Table 2).

In the DWV-A model, neither the time nor the quadratic time term was statistically significant without interaction terms (in the corresponding sections below) in explaining levels of DWV-A (Tables 1, 2). However, there was a significant positive relationship of time (and negative with the quadratic term) to DWV-B levels (Table 1; Supplementary Table 2). There was only marginal evidence for a treatment by time interaction in the DWV-A model (Table 1). In the DWV-A model fixed effects, neither PBS nor DWV treatments were significant (Supplementary Table 2). Nor we did not find any evidence, even marginal, for a treatment by time interaction in the DWV-B model (Table 1; Supplementary Table 2). However, virus levels varied with other variables and interaction combinations over time (in the corresponding sections below).

When we looked at overall DWV-A and DWV-B viral levels (controlling for treatment, stock, and time), the tissue type variable was significant (Figure 3; Supplementary Figure 2; Table 1), with all tissue types having significant fixed effects (Supplementary Table 2). However, we saw opposite interactions between the two virus genotypes. In the DWV-A model, legs had the greatest virus levels while abdomens had the lowest levels among all examined tissues (Table 2). Alternatively, for DWV-B, abdomens had the greatest level while legs had the least (Table 2). For each virus genotype, head tissue consistently had relatively high levels of the virus compared to other tissues, making it potentially useful as an indicator of both DWV virus genotypes.

When tissue type by treatment interactions were analyzed, we found significant relationships for both virus models (Figure 3; Table 1), primarily due to the DWV treatment (Supplementary Table 2). The DWV treatment had higher levels of DWV-A than did PBS and control treatments for each tissue type (post-hoc tests not in Table 2). A similar pattern occurred for DWV-B; however, tissue types with PBS treatments had greater virus levels than did their associated control treatments (post-hoc tests not in Table 2).

We did not find any tissue type by bee stock interactions for either DWV-A or DWV-B models (Table 1; Supplementary Table 2) though there was a significant interaction of tissue type with time (Figure 3; Table 1). For the DWV-A model, legs were driving this trend (Supplementary Table 2), with the ending virus levels (day 10) being even lower than the starting levels (day 1). There was also an earlier (though not significantly so) peak at day 4 (Table 2). For the DWV-B model, all tissue types were significant (Supplementary Table 2), with legs having a significantly later peak than hypopharyngeal glands and heads (Table 2). We did not see any significant triple interactions for either virus using tissue type in combination with treatment or time in relation to stock; this makes sense given the lack of an initial interaction of stock and tissue type.

When we compared genetic stocks using only the control treatment (pooled across time), we found no differences in DWV-A levels (ANOVA, DF = 4, 883, F = 1.961, and p = 0.098), though numerically, Russian bees had the highest DWV-A levels (6.589 ± 0.105 log-transformed functional titer level) while Pol-Line (6.205 ± 0.105) had the lowest levels. When we separated the control treatment bees further by time, we found stock differences in DWV-A levels only for day 1 (ANOVA, DF = 4, 178, F = 10.645, and p < 0.001) but not for any other time point (ANOVA, Day 2: DF = 4, 171, F = 1.189, and p = 0.318; Day 4: DF = 4, 175, F = 1.040, and p = 0.388; Day 7: DF = 4, 178, F = 2.024, and p = 0.093; Day 10: DF = 4, 177, F = 1.756, and p = 0.140). On day 1, Russian bees had significantly higher DWV-A levels than all other stocks (Figure 4; Tukey HSD, α = 0.05). When we analyzed DWV-B levels in control bees for all times combined, we found significant differences among stocks (ANOVA, DF = 4, 855, F = 15.236, and p < 0.001). Using Tukey HSD comparisons combining all times (α = 0.05, letters denote significant differences), Russian (7.864 ± 0.156 a) and Saskatraz bees (7.677 ± 0.155 a) had significantly lower DWV-B levels than Carniolan (8.704 ± 00.157 b), Italian (8.717 ± 0.155 b), and Pol-Line bees (9.085 ± 0.156 b). When we separated the control treatment bees by time, we found stock differences in DWV-B levels for all days (Tukey HSD, α = 0.05). On day 1, Russian bees had higher DWV-B levels than Pol-Line and Saskatraz. However, by day 2, Pol-Line bees had the highest levels but were only significantly higher than Saskatraz. This trend further solidified on day 4, where Pol-Line and Italian bees had higher levels of DWV-B than Carniolan, Russian, and Saskatraz bees. Day 7 was very similar, but with Carniolan bees joining the ranks of Pol-Line and Italian. On the final day, day 10, Pol-Line had similar DWV-B levels to Carniolan but higher levels than all other stocks and Carniolan had higher levels than Russian bees.

Without considering interaction terms, when we grouped all treatments, the genetic stock of the bees had no impact on the levels of either DWV-A or DWV-B (Table 1). When looking at the fixed effects of the DWV-A model (Supplementary Table 2), Russian bees were marginally positively correlated with increased DWV-A virus levels. This trend does not hold for DWV-B. In general, the Russian stock tended to have higher levels of DWV-A while Pol-Line had the lowest with Italian, Carniolan, and Saskatraz stocks middling (Table 2), though none were significantly different. The opposite (though again not statistically significant) occurred for DWV-B, with Pol-Line having the higher trending levels and Russian bees having the lowest. This suggests potential viral genotype by host genotype interactions between viral genotypes.

There was a treatment by stock interaction in DWV-A (Table 1), primarily driven in the fixed effects (Supplementary Table 2) by the Pol-Line and DWV treatment when compared to the intercept (Italian Control). However, we did not see any treatment by stock interactions in the DWV-B model (Figure 5). Even in the model fixed effects, no stock/treatment combinations were even marginally significant. For the stock by time interaction variable in the DWV-A model, the Russian stock exhibited significant fixed effects within the model (Supplementary Table 2) though the overall interaction variable was not significant (Table 1; Figure 5; Supplementary Figures 3–6). Alternatively, the DWV-B did have a significant stock by time interaction (Table 1) where DWV-B levels peaked at different times based on stocks (Table 2).

When we considered the stock, treatment, and time interactions, DWV-A had significant effects while DWV-B did not (Table 1). When we looked at the DWV-A model fixed effects, only stock interactions with the DWV treatment were significant over time (Supplementary Table 2), with all stocks but Carniolan being significant. When looking at Figure 5A, we saw few stock differences in DWV-A levels over time within the control treatment. In the PBS treatment, we observed that Pol-Line and Saskatraz appeared to increase over time with a positive peak; but Russian, Italian, and Carniolan stocks had a smaller increase or even a decrease in DWV-A over time. For DWV treatments, we observed that all stocks but Pol-Line peak and start decreasing in DWV-A over time. Pol-Line started on the lower end of DWV-A levels and increased linearly over time, with day 10 levels still being at the lower end of DWV-A levels. When we looked at the same three-way interaction term for DWV-B levels (Figure 5B), we did not see a significant relationship (Table 1), with no significant interactions in the fixed effects (Supplementary Table 2).