Honey Bee Genetic Stock Determines Deformed Wing Virus Symptom Severity but not Viral Load or Dissemination Following Pupal Exposure

Honey bees exposed to Varroa mites incur substantial physical damage in addition to potential exposure to vectored viruses such as Deformed wing virus (DWV) that exists as three master variants (DWV-A, DWV-B, and DWV-C) and recombinants. Although mite-resistant bees have been primarily bred to mitigate the impacts of Varroa mites, mite resistance may be associated with increased tolerance or resistance to the vectored viruses. The goal of our study is to determine if five honey bee stocks (Carniolan, Italian, Pol-Line, Russian, and Saskatraz) differ in their resistance or tolerance to DWV based on prior breeding for mite resistance. We injected white-eyed pupae with a sublethal dose (10⁵) of DWV or exposed them to mites and then evaluated DWV levels and dissemination and morphological symptoms upon adult emergence. While we found no evidence of DWV resistance across stocks (i.e., similar rates of viral replication and dissemination), we observed that some stocks exhibited reduced symptom severity suggestive of differential tolerance. However, DWV tolerance was not consistent across mite-resistant stocks as Russian bees were most tolerant, while Pol-Line exhibited the most severe symptoms. DWV variants A and B exhibited differential dissemination patterns that interacted significantly with the treatment group but not bee stock. Furthermore, elevated DWV-B levels reduced adult emergence time, while both DWV variants were associated with symptom likelihood and severity. These data indicate that the genetic differences underlying bee resistance to Varroa mites are not necessarily correlated with DWV tolerance and may interact differentially with DWV variants, highlighting the need for further work on mechanisms of tolerance and bee stock–specific physiological interactions with pathogen variants.

1 Introduction

Honey bee colony health is substantially impacted by the presence of the parasitic mite Varroa destructor (Anderson and Trueman), known to vector honey bee viruses such as Deformed wing virus (DWV) (Conte et al., 2010; Rosenkranz et al., 2010; Nazzi and Conte, 2016). Overt morphological symptoms of DWV, of which there are several prevailing variants (A, B, and C) and recombinants, include malformed wings and bloated abdomens (Moore et al., 2011; Zioni et al., 2011; Martin et al., 2012; Mordecai et al., 2016a; Mordecai et al., 2016b). Symptomatic adults die quickly upon emergence, while asymptomatic (i.e., no visible morphological changes) adults exhibit learning difficulties and diminished flight ability (Javaid and Uli, 2007; de Miranda and Genersch, 2010; Wells et al., 2016; Benaets et al., 2017). Such effects on individual bees may accumulate to impact whole colonies, with DWV infection decreasing the probability of colony survival (Berthoud et al., 2010; Dainat et al., 2012; Martin and Brettell, 2019).

One strategy to mitigate Varroa mite impacts on honey bee colony health and survival is breeding for mite resistance (Buchler, 1994; Harbo and Harris, 1999a; Rinderer et al., 2010; Mondet et al., 2020). Previous honey bee breeding programs have selected several parameters such as hygienic and grooming behaviors meant to decrease the mite load in the colony (Harbo and Harris, 1999b; Guzman-Novoa et al., 2012; Bąk and Wilde, 2015; Guichard et al., 2020; Mondet et al., 2020). Hygienic behavior, known to reduce the levels of diseases such as American and European foulbrood by the removal of damaged pupae, has also been shown to reduce Varroa mite loads in European honey bee colonies (Rothenbuhler, 1964; Spivak and Reuter, 2001a; Spivak and Reuter, 2001b). Varroa-sensitive hygienic bees (VSH) including the Pol-Line stock are a sub-selected group exhibiting hygienic behavior that targets mite-infested pupal cells rather than damaged pupae, removing reproductive Varroa mites, minimizing mite population growth (Harris, 2007; Rinderer et al., 2010; Danka et al., 2013; Danka et al., 2016; Mondet et al., 2016; Saelao et al., 2020; Spivak and Danka, 2021). Another mite-resistant stock, the Russian honey bee, was bred starting from a population historically associated with Varroa mites and was selected more generally for low mite population growth, which involves several mechanisms of defense such as brood breaking and heightened grooming (de Guzman et al., 2007; Rinderer et al., 2010). Similarly, the Saskatraz stock has also undergone selection primarily for survivorship and low mite population growth without targeting a single, specific mechanism (Robertson et al., 2014).

Furthermore, the genotypic differences underlying mite-resistant behaviors may be associated with reduced levels of mite-vectored viruses such as DWV (Locke et al., 2014; Khongphinitbunjong et al., 2016; de Guzman et al., 2017; Mendoza et al., 2020; Weaver et al., 2021; O’Shea-Wheller et al., 2022). Such bee genotype × DWV interactions may be the result of enhanced virus resistance or tolerance in certain bee stocks, with prior data indicating a greater potential for tolerance (Strauss et al., 2013; Locke et al., 2014; Khongphinitbunjong et al., 2015; Khongphinitbunjong et al., 2016; Thaduri et al., 2018; Thaduri et al., 2019; Locke et al., 2021). Resistance is the ability of the bee to prevent an infection from establishing or increasing after exposure, while tolerance is the ability of the individual to maintain health and functionality (i.e., reduced symptoms), even while having an infection (Locke et al., 2014; Mordecai et al., 2016a; Burgan et al., 2018). Host genotype × virus interactions may be direct through individual host immune responses to the virus or indirect via colony-level mite-removal traits (Chen and Siede, 2007; Khongphinitbunjong et al., 2015; de Guzman et al., 2019; Bouuaert et al., 2021; Weaver et al., 2021). For instance, genotype may mediate the bee’s ability to detect signals such as mite-related cuticular hydrocarbons; then, these detection differences can translate into differential hygienic responses resulting in cascading differences in mite and vectored virus levels (Mondet et al., 2021; Wagoner et al., 2019). In addition, virus levels may potentially feed back into this system by further disrupting the bee’s mite detection ability (Wagoner et al., 2019; Mondet et al., 2015). If direct mechanisms drive genotype-associated differences in virus levels, we expect to observe genotype differences in laboratory studies of virus loads as the need for mite-removal behavior is not necessary for virus reduction (Locke et al., 2014). Complicating this matter are potential differences among DWV variants, with conflicting reports indicating context-dependent variant virulence and likelihood of inducing symptoms or mortality (Mordecai et al., 2016a; McMahon et al., 2016; Brettell et al., 2017; Natsopoulou et al., 2017; Gisder et al., 2018; Barroso-Arévalo et al., 2019; Kevill et al., 2019; Tehel et al., 2019; Dubois et al., 2020; Norton et al., 2020).

In this study, we evaluated potential resistance or tolerance to DWV in five honey bee genetic stocks (Carniolan, Italian, Pol-Line, Russian, and Saskatraz), with varying levels of resistance to the Varroa mite using a laboratory environment to control for indirect behavioral mechanisms. The Italian and Carniolan stocks are commonly used throughout the commercial beekeeping industry and have been bred for honey production and colony size; whereas Russian, Pol-Line, and Saskatraz bees have been specifically bred for Varroa mite resistance (de Guzman et al., 2007; Robertson et al., 2014; Danka et al., 2016; Caron, 2021; O’Shea-Wheller et al., 2022). These genetic stock differences represent different genotypes in a broad sense as they have clear differences in phenotype and are derived from different breeding populations. The specific genotypic information was not further assessed as part of this study though it can be found, in part, in prior studies (Jiang et al., 2016; Wang, 2016; Saelao et al., 2020). To determine potential interactions of bee stock and DWV, we exposed honey bee pupae, the life stage that is most commonly infested with mites (Ifantidis, 1983, 1984; Donzé and Guerin, 1994), to Varroa mites or injected them with a sublethal dose (10⁵) of DWV sourced from symptomatic (e.g., deformed wings) adult bees. Upon adult emergence, we determined DWV levels and dissemination throughout different tissue types [abdomen, head, hypopharyngeal gland, and a rear leg as in the study by Penn et al. (2021)], the number of days until adult emergence, and DWV symptom presence and severity. The tissue types selected were chosen as legs have been used to indicate viral dissemination in other arthropods (Boncristiani et al., 2009; Diagne et al., 2015); the head has been an indicator of overt bee infections (Yue and Genersch, 2005; Möckel et al., 2011); hypopharyngeal glands may provide possible transmission by food trophallaxis (Chen Y. et al., 2006, Chen et al., 2006 Y. P.; Möckel et al., 2011); and the abdomen as the site of mite feeding and our injection treatment groups. Since there are two common variants (A and B) of DWV in the United States and the potential virulence and mortality effects may differ based on interactions between host and viral genotypes (Dainat et al., 2012; Ryabov et al., 2017; Kevill et al., 2019; Loope et al., 2019; Grindrod et al., 2021), we tested for RNA copy levels of both DWV-A and DWV-B variants.

2 Materials and Methods

2.1 Source Colonies

All honey bee colonies were started from 0.90 to 1.10 kg “packages” made on 3 May 2018 from 10 established colonies of the USDA Honey Bee Breeding, Genetics and Physiology Research Unit in Baton Rouge, LA, United States (30°22′56″ N, 91°10′40″ W). Naturally mated queens from five genetic stocks were sourced from the USDA laboratory (Pol-Line and Russian), a Canadian collaborating breeder (Saskatraz), or purchased from commercial California suppliers (Carniolan and Italian). Queens were queen candy released into colonies on 4 May 2018. Colonies (N = 3 colonies/stock) were not sampled until 6 weeks post queen introduction and supplementation to allow time for population turnover to reflect queen genetics. All colonies were similarly maintained in three yards near the USDA laboratory (with Carniolan, Italian, and Saskatraz sharing one yard, while Pol-Line and Russian colonies were maintained in two separate yards). To allow for direct comparison, the same colonies were used in a complementary study following DWV levels over 10 days in newly emerged adult bees (Penn et al., 2021).

2.2 Viral Isolation and Mite Sourcing

To obtain the DWV viral solution for injection, 20 adult honey bees with overt DWV symptoms were frozen at −80°C, ground to a fine powder, homogenized in 10 ml 1X PBS, and centrifuged at 5,000 rpm for at 4°C for 20 min. The resulting supernatant containing viruses was filtered through a 0.2-micron filter (milex-GS syringe filter unit #SLGS033SS, Millipore Sigma, Burlington, MA, United States) to remove small tissue debris, fungi, and bacteria. qPCR was conducted to test for 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 mentioned later. DWV quantification using general DWV (not variant-specific) primers (Supplementary Table S1) was performed by absolute quantification using the standard curve method. All methods were previously established based on standard protocols (Simone-Finstrom et al., 2018). One sample stock solution was selected based on negative results for non-target viruses and used to create the injection stock solution. The stock solution was diluted to 10⁵ viral copies of DWV, a biologically relevant but sublethal to adult bees (Gisder et al., 2009). For the mite inoculation treatment group, Varroa mites were obtained from non-study hives using powdered sugar rolls to dislodge live mites from nurse bees (Macedo et al., 2002). Dislodged mites were removed from powdered sugar using a paint brush and stored in a petri dish containing a maximum of 100 mites, allowed to feed on 10 previously un-infested bee pupae (exchanged daily), and maintained in the incubator at 34°C and 85% relative humidity until use within 24–48 h (Dietemann et al., 2013).

2.3 Pupal Assay

From mid-July to early October 2018, frames with white-eyed pupae were removed from each colony and brought back to the laboratory where the Varroa mite-free (no mites in pupal cells or on body) pupae were removed, placed on a folded filter paper in petri dishes, and stored in an incubator at 34°C and 85% relative humidity. Pupae exhibiting any discoloration due to damage during handling after 2 h were removed from the experiment, while healthy pupae (17 treatment/colony) were assigned to one of four treatment groups—1) no manipulation (included since bees had naturally occurring DWV infections, referred to as “control”), 2) 3.0 µl 1X PBS injection to control for injection damage (PBS), 3) 3.0 µl 10⁵ DWV injection to simulate the vectoring of DWV without mite presence (DWV), or 4) Varroa mite inoculation (mite). Control bees were placed into an individual, size “1” gel capsules (Capsule Connection, Prescott, AZ, United States), with a small hole in the top created using an insect pin (Nazzi and Milani, 1994; Piou et al., 2016; Posada-Florez et al., 2019). All pupae in the PBS and DWV treatment groups were injected using an UltraMicroPump with an SYS-Micro4 controller (World Precision Instruments, Sarasota, FL, United States) with an infusion flow rate of 1.0 μl/s, following the manufacturer’s parameters. For each injection, a 30G needle (Hamilton Company, Reno, NV, United States) was inserted into the lateral abdomen between the fourth and fifth pleurites, based on established protocols (Simone-Finstrom et al., 2018). After injection, pupae were then transferred to individual gel capsules as mentioned earlier. For the mite treatment group, pupae were transferred to an individual gel capsule followed by an individual mite (Dietemann et al., 2013).

All pupae were incubated at 34°C and 85% relative humidity until adult emergence, with all bees checked daily and necrotic individuals removed. Mite treatments where the mite was dead or did not defecate (an indication that no feeding occurred) by adult bee emergence were not used in subsequent analysis but were replaced with replicates meeting the requirements. Upon emergence, all adult bees were evaluated for DWV symptoms and rated on a scale of severity from 0 to 3 (Supplementary Figure S1), with 0 indicating normal wings, 1 indicating slight malformation, 2 indicating major malformations but with wings present, and 3 indicating completely malformed. Mites were removed from respective emerged bees; then all mites and bees were placed into individual sterile 1.5-ml centrifuge tubes and stored at −80°C. The first emerging bees per colony per stock per treatment group combination were evaluated in this manner (n = 12–17 bees/colony/treatment group/stock, N = 981 total emerged bees) and then a random subset of three bees was selected for tissue dissection and DWV quantification (Figure 1).

FIGURE 1

FIGURE 1. Experimental design: experimental design to determine DWV symptoms and levels based on bee stock and treatment of white-eyed pupae (no injection control/naturally occurring infection, PBS sham injection, Varroa mite exposure, or sublethal DWV injection). DWV symptoms were analyzed in emerged adults (N = 12–17 bees per stock/treatment combination); then a subset of three bees was dissected into four tissue types (heads, hypopharyngeal glands, abdomen, and rear legs) subsequently used for viral analyses.

2.4 Adult Assay

The adult assay was part of a complementary study on DWV dissemination in injected adult bees over time (Penn et al., 2021). The same colonies were used for the adult injection as the pupal injection experiment and throughout the same timeframe (July through October 2018). Varroa mite–free newly emerged bees were treated similarly, except for the no mite treatment group, and all bees were placed on ice for 2 m prior to injection to reduce movement. Bees were injected using the same methods as aforementioned with a 3.0 µl of DWV inoculum from an aliquot of the same inoculum as used for the pupae experiment (DWV) or 3 µl 1X PBS (PBS) or had no injection (control) were implemented as controls. The bees were housed in cages of maximum 30 bees and provided with 50% sucrose solution and pollen substitute. The bees were maintained in an incubator at 34°C and 85% relative humidity and then sacrificed at 1, 2, 4, 7, and 10 days post injection. Only day seven data (N = 135 bees and 540 associated RNA extractions) were included within this study as 7 days allows for a similar period of viral replication compared to pupae data.

2.5 Dissection and RNA Isolation

To determine virus dissemination through the bee body over time, the randomly selected three-bee subset from each colony/treatment group/stock combination was dissected in the same manner as the complementary study on adult bees (Penn et al., 2021). Dissections were conducted over dry ice with individuals dissected with a new, sterilized blade. The body was separated into legs, wings, head, thorax, and abdomen; the head was embedded into beeswax (new per individual bee) and the hypopharyngeal gland was removed according to the previously published methods (Corby-Harris and Snyder, 2018; Penn et al., 2021). Dissected tissues were stored in separate sterilized tubes on dry ice during dissection and long term at −80°C.

Total RNA was extracted for a single rear leg, the head (sans hypopharyngeal gland), hypopharyngeal gland, and abdomen for each of the three bees representing each combination of colony/treatment group/stock (N = 180 bees total and N = 720 RNA extractions in total) and all mites (N = 45). Individual hypopharyngeal glands and mites were placed in 30 µl lysis buffer and 30 μl Maxwell homogenization buffer and vortexed. The leg (cut into pieces), head, and abdomen were placed in 200 µl lysis buffer and 200 μl Maxwell homogenization buffer (Promega Corporation, Madison, Wisconsin, United States), manually ground with a pestle (Sigma-Aldrich, St. Louis, Missouri, United States), and vortexed. All samples were then incubated for 90 min at 4°C. After incubation, 320 μl Maxwell homogenization buffer was added to the hypopharyngeal gland and mite samples. The samples were then extracted using the Maxwell RSC 48 cartridges (Promega Corporation, Madison, Wisconsin, United States). The total RNA was extracted according to standard procedures using RSC simplyRNA tissue extraction kits and program (Promega Corporation, Madison, Wisconsin, United States). RNA was stored in 0.6-ml elution tubes wrapped in parafilm (Bemis NA, Neenah, Wisconsin, United States) at −80°C.

2.6 cDNA Synthesis and RT-PCR

Frozen RNA samples were thawed on −20°C metal beads, briefly vortexed, and then centrifuged. Each RNA sample was nano-dropped (NanoDrop One, Thermo-Fisher Scientific Inc., Waltham, Massachusetts, United States) twice using 1 μl of sample. The mean ng/µl NanoDrop One readings (for 260/280 and 260/230) were calculated per sample and then used to determine adequate sample purity and the quantities of sample and nuclease-free water required to dilute each sample to a concentration of 250 ng of RNA. cDNA was then synthesized in two steps using Qiagen QuantiTect Reverse Transcription kits (Thermo-Fisher Scientific Inc., Waltham, Massachusetts, United States). For step 1, 2 µl of gDNA wipeout buffer was added to the mix of RNA and water for a total reaction volume of 14 µl per sample. The samples were incubated at 42°C for 2 min in a Bio-Rad T100 Thermal Cycler (Bio-Rad, Hercules, California, United States). The samples were then 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. The samples were again vortexed and centrifuged and then placed into the Bio-Rad T100 Thermal Cycler (42°C for 25 min then 95°C for 3 min). cDNA was stored in strip tubes wrapped in parafilm at −80°C.

All samples were tested with DWV-A and DWV-B primers using qRT-PCR to determine infection levels (primers in Supplementary Table S1), and each sample was replicated two times per primer pair. All qRT-PCRs consisted of 5 µl SsoFast Universal SYBR Green Supermix (Bio-Rad, Hercules, California, United States), 3 µl nuclease-free water, 0.5 µl forward primer, 0.5 µl reverse primer, and 1 µl cDNA from each sample. All reactions were run in Bio-Rad CFC 96 or Connect thermal cyclers (Bio-Rad, Hercules, California, United States), with all reactions of a specific primer pair occurring in the same machine. 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 the product size. DWV-A and DWV-B results were quantified using the standard curve method using linearized plasmid constructs (ranging from 10⁵–10¹²). Quantified virus RNA copy levels (DWV RNA equivalents per 10 ng of RNA) were log-transformed for analyses.

2.7 Statistical Analyses

All statistical analyses were conducted in R 4.0.2 (R Core Team, 2020), and all graphs were plotted using ggplot2 (Wickam, 2016). To determine factors influencing DWV-A and DWV-B levels and dissemination to different tissue types, the three bee sub-samples with virus data were used with all tissue types included in the analyses. Separate general linear mixed models (GLMMs) were generated for each DWV variant (DWV-A and B) using the lme4 (glmer function with Gaussian distribution) package (Bates et al., 2015). Variables for model selection included bee stock, treatment group, tissue type, all interaction terms of those three variables, and the levels of the alternate DWV variant. Colony and individual bees were used as random effects to account for multiple tissue types coming from the same individual and individuals from the same colony. The Italian bee stock, control treatment group, and head tissue were specified as model intercept values. Backward model selection was conducted using a combination of minimum AICc and BIC values; models with significantly lower scores (Δ > 4) were used. It is to be noted that the stock × treatment group interaction term was tested in model selection and was neither significant nor did it contribute to model fit so was not included in the final model. Final model information can be found in Supplementary Table S2. Significance values were estimated using Satterthwaite approximation with the lmerTest package (Kuznetsova et al., 2017). Post hoc comparisons were conducted with Kenward–Roger estimation in the emmeans package (Lenth, 2020).

This experiment also allowed us the unique opportunity to evaluate DWV levels and dissemination differences between this experiment and prior work with injected newly emerged adults from the same colonies at the same time of the same year of collection (Penn et al., 2021). To determine if differences occurred between these two experiments, the mite treatment group was dropped from this pupae experiment as it was not included as a treatment group in the adult experiment. Similarly, we only evaluated the day 7 time-point from the adult experiment as that matched the pupae experiment in terms of number of days post injection. Day 7 comparison models were conducted as aforementioned except that experiment and associated interaction terms were included in the model and the adult experiment was used in the model intercept.

The time (days) until adult emergence was modeled using a GLMM using the lme4 (glmer function with Gaussian distribution) package (Bates et al., 2015). To avoid redundant sampling per bee, only head tissue (sans hypopharyngeal glands) data were used as this particular tissue type was the best representation of infections for both DWV variants (Schurr et al., 2019; Dubois et al., 2020; Penn et al., 2021). Variables considered for model selection included bee stock, treatment group, stock × treatment group, and DWV-A and -B levels; colony was used as a random effect. The Italian bee stock and the control treatment group were specified as model intercept values. Model selection and post hoc comparisons were conducted as mentioned earlier. The presence of DWV symptoms was analyzed similarly using a GLMM using the lme4 (glmer function with binomial distribution) package but with a logit link function within a binomial logistic regression (presence/absence of symptoms).

DWV symptom severity (scale of 0–3) was modeled using an ordered logistic regression using the polr function (logistic) in MASS (Ripley et al., 2018). Marginal effects and predicted probabilities for the treatment group and stock comparisons were calculated using the effects package (Fox, 2003; Fox and Hong, 2009; Fox and Weisberg, 2019). Variables and methods for model selection were as mentioned previously (time to emergence and symptom presence). We further synthesized the results of the symptom presence and severity models using classification trees. The two-response categorical variable for DWV symptom presence (0, 1) and the four-response ordinal variable for DWV symptom severity (wind deformity on a scale of 0–3) were used to construct two classification trees with the same predictor variables as the associated logistic models. The classification trees were created and plotted using recursive partitioning provided in the rpart and rpart plot packages (Milborrow, 2017; Therneau et al., 2018). All predictor variables were maintained in the final classification trees with each having a complexity parameter >0.01 after pruning. The relationships among bee stock, treatment group, DWV levels, and the DWV symptom severity scale were further modeled using multiple correspondence analysis (MCA) and factor analysis of mixed data (FAMD) using the FactoMineR package (Husson et al., 2008). MCA allows for the analysis of multiple qualitative variables (similar to PCA with quantitative data), whereas FAMD analyzes the relationships between the combinations of qualitative and quantitative variables. Variables presented closer together on the plot are more similar to each other than those further away (like on the opposite side of the origin).

For ease of graphical representation, PBS, mite, and DWV pupae treatment groups were also compared against the control treatment group with regards to log-transformed DWV-A and DWV-B levels, days until emergence, and symptom severity. This was performed by calculating the mean for the control data for each of the three replicates per tissue type, treatment group, and colony. The mean control value was assumed to be the baseline (y = 0 value on Figure 2) value for that colony × tissue type × treatment group combination and was then subtracted from each of the three replicates per treatment group for that colony and tissue type. Statistical significance represented in Figure 2 was determined by using post hoc comparisons for the treatment group, described earlier.

FIGURE 2

FIGURE 2. Treatment values relative to naturally occurring levels. Comparison of (A) log-transformed DWV-A and (B) DWV-B levels in head tissues indicates that the pupae DWV injection treatment sourced from symptomatic adults successfully increased DWV-A but not DWV-B levels in all bee stocks. (C) Treatment × bee stock interactions decreased the number of days from pupal treatment until adult emergence compared to the control. (D) DWV symptom severity generally increased in the DWV treatment but interacted with bee stock. Values were made relative to the untreated control bees from the same colony with naturally occurring infections (y = 0). Boxplots are in the style of Tukey where the box limits represent the lower 25% quantile and upper 75% quantile with the line representing the median. Individual points indicate individual bees (N = 9 bees/treatment/stock) with shapes associated with colony replicate (ColonyRep) and color indicating stock. * Denotes significant differences between the treatment and control values (Dunnett’s test; p < 0.05).

Graphical representation of DWV dissemination (Figures 3, 4) was conducted in Microsoft 365 Excel and PowerPoint. The mean log-transformed DWV-A and DWV-B values (across all bee stocks) per treatment group were calculated and copied and pasted into a column in an Excel spreadsheet with an associated treatment group and variant labeling in the prior column. Whole numbers (as represented by the color scale in Figures 3, 4) were placed through this column and encompassed the minimum and maximum values of both the DWV-A and DWV-B values. The DWV values were then selected and colored using the conditional formatting tool where the lowest whole number was assigned blue (RGB: 0, 112, 192; Hex #0070C0) and the largest whole number was assigned red (RGB: 255, 0, 0; Hex #FF0000). This color scale and precise values based on the tissue type × treatment group means were pasted into PowerPoint as an image to maintain the appropriate color scale. Graphical outlines of bees were generated by the authors in PowerPoint, and each tissue type per treatment group was colored using the corresponding DWV values with the eyedropper tool.

FIGURE 3

FIGURE 3. DWV dissemination in newly emerged adult bees. Mean log-transformed DWV-A and DWV-B values (color values scaled for overall DWV levels) in each sampled tissue type (abdomen, head, hypopharyngeal gland, and rear leg) for each treatment (uninjected control, PBS injection, Varroa mite exposure, and DWV injection). Dissemination patterns differed between the two DWV variants but was not related to bee stock (displayed DWV values were averaged across bee stocks).

FIGURE 4

FIGURE 4. Comparison of DWV dissemination in treated pupae and adults. Mean log-transformed DWV-A and DWV-B values (color values scaled for overall DWV levels) in each sampled tissue type (abdomen, head, hypopharyngeal gland, and rear leg) for each treatment (uninjected control, PBS injection, and DWV injection) overlapping between treated pupae and treated newly emerged adult experiments. The pupae experiment represents data described in this manuscript (bees injected as white-eyed pupae but sampled upon adult emergence about 7–8 days post-treatment), while the adult data represent samples of injected newly emerged adults collected 7 days post treatment Supplementary Figure S2 for the comparison of DWV for the combination of stock, treatment, and experiment.

3 Results

3.1 Deformed Wing Virus Levels and Dissemination

3.1.1 Treatment Impacts on Viral Load

When we analyzed DWV levels in the newly emerged adult bees, we found that the treatment group was significant but differed with the DWV variant (Table 1; Supplementary Table S3). Overall, DWV-A levels were higher in the DWV-injected bees than in the untreated control bees with naturally occurring infection (Kenward–Roger estimation, t = 9.279, p < 0.001), PBS (t = 6.061, p < 0.001), or mite treatments (t = 7.595, p < 0.001) (Figure 2). Unlike DWV-A, PBS injection and mite infestation resulted in higher DWV-B levels than those of the untreated bees (Kenward–Roger estimation, compared to PBS: t = 5.208, p < 0.001; compared to mite: t = 4.463, p < 0.001) and the DWV injection treatment group (compared to PBS: t = 3.218, p = 0.008; compared to mite: t = 2.489, p = 0.065). Moreover, the DWV treatment group did not differ in DWV-B levels from the untreated control group (Kenward-Roger estimation, t = 1.664, p = 0.346). Mites used for the infestation treatment group had significantly lower levels of DWV-A (10^{2.42 ± 0.34}) than the DWV-B level (10^{4.82 ± 0.20}) (ANOVA, F_{(1, 87)} = 36.735, p < 0.001), but mites expressed similar levels of each variant regardless of bee stock that they were infesting (DWV-A: ANOVA, F_{(4, 43)} = 0.642, p = 0.636; DWV-B: ANOVA, F_{(4, 43)} = 0.556, p = 0.696).

TABLE 1

TABLE 1. DWV dissemination ANOVAs. ANOVA summary of the generalized linear mixed-model fixed effects for the log-transformed DWV-A and DWV-B RNA copy levels disseminated throughout the newly emerged adult tissues (N = 180 bees total and N = 720 RNA extractions in total). Supplementary Table S1 for model parameter estimates.

3.1.2 Deformed Wing Virus Levels in Different Tissues

Bee stock did not significantly impact overall DWV-A levels (Table 1; Supplementary Table S3; Figure 2), indicating that genetic stocks are not resistant to DWV, given the exposure circumstances of this experiment. However, both DWV variants varied significantly but differentially among the tissue types (Table 1; Supplementary Table S3; Figure 3). Pooling all treatment groups, for DWV-A, heads and legs had higher levels than abdomens (Kenward–Roger estimation, compared to heads: t = −7.224, p < 0.001; compared to legs: t = −9.279, p < 0.001) and hypopharyngeal glands (compared to heads: t = −7.250, p < 0.001; compared to legs: t = −7.406, p < 0.001). With DWV-B levels, heads had higher levels than all other sampled tissues (compared to abdomen: t = 7.750, p < 0.001; compared to hypopharyngeal gland: t = 4.756, p < 0.001; compared to leg: t = 4.832, p < 0.001), and legs had higher levels than abdomens (t = 3.325, p = 0.005) but were not different from hypopharyngeal glands (t = 0.784, p = 0.862). Overall, heads and legs had the highest viral loads of both DWV genotypes.

Tissue type significantly interacted with the treatment group for both DWV-A and B levels (Table 1; Supplementary Table S3; Figure 3). Within the control treatment group, heads, legs, and hypopharyngeal glands had similar levels of DWV-A; but heads and legs exhibited higher DWV-A levels than those of the abdomens (Table 2). For DWV-B in the control treatment group, heads had the greatest virus levels and legs had higher levels than those of abdomens but were not different from hypopharyngeal glands (Table 2). In the PBS treatment group, both heads and legs had higher DWV-A levels than abdomens and hypopharyngeal glands (Table 2). Abdomens had lower DWV-A levels than all other tissue types in the mite treatment group (Table 2). For DWV-B in the PBS treatment group, the only significant difference was between heads and hypopharyngeal glands; and there were no observed differences in DWV-B levels among tissue types in the mite treatment group (Table 2). In the DWV injection treatment group, hypopharyngeal glands had lower DWV-A levels than all other tissue types (Table 2). For DWV-B levels in the DWV treatment group, abdomens had lower levels than all other dissection components levels, which were not different from each other (Table 2).

TABLE 2

TABLE 2. Tissue-type contrasts. Post hoc contrasts comparing DWV-A (A) and DWV-B (B) levels among tissue types [head, hypopharyngeal gland (hyp. gland), rear leg (leg), and abdomen] for each treatment [naturally occurring infection (control), PBS injection, DWV injection and Varroa mite exposure (mite)]. Statistical significance is also indicated by “*”.

3.1.3 Deformed Wing Virus Levels in Injected Pupae and Adults After 7 Days

In order to compare DWV levels and dissemination differences based on the life stage (pupae or newly emerged adult) at the time of treatment, we modeled DWV-A and DWV-B levels for this experiment in conjunction with data from a complementary study on adult bees (Penn et al., 2021). While different cohorts of bees were used within each colony for the pupae and adult studies, the sampling of these cohorts occurred within the same months. In the adult bee study, newly emerged bees from the same colonies and time of year were injected with different aliquots of DWV from the same inoculum at the same concentration. Comparisons of data from the newly emerged bees treated as pupae in the current study were made to the data from 7 days post injection in the adult study as this timeline allows for a similar period of viral replication between the two datasets. We found that for both DWV-A and DWV-B, the experiment (adult vs. pupae injection) was significant, interacted with the treatment group and tissue type individually (Figure 4; Supplementary Figure S2), and interacted with the treatment group and tissue type in a three-way interaction (Table 3). Within each treatment group, pupae had significantly greater levels of DWV-A than adults (Kenward–Roger estimation, control: t = −2.921, p = 0.004; PBS: t = −3.151, p = 0.002; DWV: t = −4.760, p < 0.001). However, this differed for DWV-B levels, where adults had higher amounts than pupae within all treatment groups but was only significant in the DWV treatment group (control: t = 1.795, p = 0.074; PBS: t = 1.349, p = 0.179; DWV: t = 8.150, p < 0.001). The experiment also exhibited a three-way interaction with bee stock and treatment for DWV-A but did not interact with bee stock for DWV-B (Table 3).

TABLE 3

TABLE 3. Day 7 comparison ANOVAs. ANOVA summaries for log DWV-A and DWV-B linear mixed models analyzing day 7 values for pupal and adult-injected bees (experiment). For pupal injections, the mite treatment was excluded; all other treatments remained in the analysis (control, PBS injection, and sublethal DWV injection). Each experiment was conducted on the same three colonies of each of the five stocks with the same four tissue types analyzed.

3.2 Time Until Adult Emergence

After treatment, pupae were checked every 24 h for adult emergence. The time (days) until adult emergence was impacted by the treatment group, stock × treatment group interactions, and DWV-B levels (Tables 4, 5) and was monitored. In the model, the mite treatment group was significantly negatively associated with the number of days until emergence (Table 5; Supplementary Table S4). Pol-Line bees assigned to the mite treatment group took less than half a day longer to emerge than Italian (Kenward–Roger estimation, t = 3.332, p = 0.022), Carniolan bees (t = 3.227, p = 0.027), and Saskatraz (t = 2.822, p = 0.065) bees exposed to mites (Table 5; Supplementary Table S4). In contrast, higher DWV-B levels decreased the number of days until adult emergence (Table 5; Supplementary Table S4), so bees with higher DWV-B developed more quickly.

TABLE 4

TABLE 4. Mean days to emergence. The mean number of days until adult emergence (± standard error of the mean) for all observed bees and the nine-bee subset (per stock and treatment combination) used for viral analyses.

TABLE 5

TABLE 5. Days to emergence ANOVA. ANOVA summary of the generalized linear mixed model fixed effects for the number of days until adult emergence. All viral data were collected from head tissues (N = 180 bees and N = 180 RNA extractions in total). Supplementary Table S1 for model parameter estimates.

3.3 Differences in Morphologically Deformed Wing Virus Symptoms

3.3.1 Symptom Presence

The treatment group and the DWV variant impacted the probability of emerging bees exhibiting any DWV symptoms (i.e., wing deformities) (Table 6). The treatment group had a significant impact on the presence of symptoms, with DWV injections (associated with increased DWV-A levels) inducing symptoms in a higher number of bees relative to all other treatment groups (compared to control: z = 3.962, p < 0.001; compared to PBS: z = 3.037, p = 0.013; compared to mite: z = 2.814, p = 0.025) (Table 6). Symptomatic bees had higher DWV-A (10^{9.71 ± 0.20}) and DWV-B (10^{10.14 ± 0.12}) levels than asymptomatic bees (DWV-A: 10^{8.45 ± 0.14}; DWV-B: 10^{9.57 ± 0.11}) (Figure 5). When DWV levels were compared between symptomatic and asymptomatic bees split out by bee stock, we found that Italian and Pol-Line bees exhibited higher DWV-A levels in symptomatic bees, while only Pol-Line bees exhibited higher DWV-B levels in symptomatic bees (Supplementary Figure S3).

TABLE 6

TABLE 6. DWV symptom presence logit. Binomial logistic regression fixed effects for the presence of any DWV symptoms (wing deformities). All viral data were collected from head tissues (N = 180 bees and N = 180 RNA extractions in total).

FIGURE 5

FIGURE 5. DWV levels given symptom presence. The log-transformed DWV-A and DWV-B levels in head tissues for all sampled bees (N = 180) indicate that bees exhibiting any morphological symptoms had higher levels of both DWV-A and DWV-B. The lighter shade represents asymptomatic bees and the darker shade represents symptomatic bees. Boxplots are in the style of Tukey where the box limits represent the lower 25% quantile and upper 75% quantile, with the line representing the median. * denotes significant differences between asymptomatic and symptomatic values for each DWV variant (Kruskal–Wallis test; p < 0.05). For stock-specific comparisons, Supplementary Figure S3.

3.3.2 Symptom Severity

When the severity of DWV symptoms (Supplementary Figure S1) was analyzed (Supplementary Table S5), bee stock, treatment group, and stock × treatment group were significant. All bee stocks were significant within the model, indicating potential stock-based differences in DWV tolerance compared to Italian bees–Carniolan (ordered logit parameter estimation, t = 15.934, p < 0.001), Pol-Line (t = 12.980, p < 0.001), Russian (t = 19.370, p < 0.001), and Saskatraz (t = 20.278, p < 0.001). Stocks varied both in the proportion of observations per symptom severity category (Figure 6) and in their associated DWV-A and -B levels (Supplementary Table S6). Within the ordered logit model (Supplementary Table S5), all treatment groups were significant relative to the control treatment group–mite (t = 20.849, p < 0.001), PBS (t = 20.240, p < 0.001), and DWV (t = 18.982, p < 0.001). DWV injections had the greatest level of symptom severity followed by the PBS and mite treatment groups. No manipulation controls had low-level symptom severity. These same trends mostly held for stock × treatment group interactions but with severity levels differing among stocks within each treatment group (Figure 6). Overall trends from the ordered logit model did a fairly good job predicting symptom severity per treatment and stock; but interpreting precise model results needs to be carried out with care and in conjunction with the following multivariate analyses for a more complete picture as no symptoms were observed for Italian bees in control treatment and Carniolan bees in the PBS treatment group, skewing the ordered logit model values.

FIGURE 6

FIGURE 6. Observed rates of DWV symptom severity. The observed rates of symptom severity for each stock and treatment combination for the subset of bees used for viral analyses (N = 9 bees per stock × treatment combination, N = 180 bees total) indicate significant bee stock × treatment interactions. Contrary to expectations, these interactions are not consistent for all Varroa mite–resistant stocks (Pol-Line, Russian, and Saskatraz). Zero represents no symptoms and 3 represents debilitating morphological symptoms based on Supplementary Figure S1. Supplementary Table S4 for mean DWV levels for each bee stock and symptom severity categorization.

3.3.3 Multivariate Analyses of Symptom Outcomes

To better understand the relative importance of bee stock, treatment group, and DWV levels for DWV symptoms, we used classification trees to analyze the probability of symptom presence and severity. We found that DWV levels of both variants in combination with bee stock and the treatment group played a critical role in predicting the probability of symptom presence (Figure 7A). In the first decision node of the symptom presence (1)/absence (0) tree, bees with DWV injection had the highest probability of symptoms (0.89). The second decision node split probabilities using DWV-B levels (cut off of 10^10.91), where bees whose heads contained levels above this cutoff were more likely to exhibit symptoms. The next two nodes indicated equal importance of DWV-A levels (cut off of 10^10.92) and bee stock. This node grouped Italian, Pol-Line, and Saskatraz bees together as being more likely to exhibit symptoms (>78% of bees likely to show symptoms) and Carniolan and Russian bees together as less likely to exhibit symptoms (>53% of bees likely to not exhibit morphological symptoms).

FIGURE 7

FIGURE 7. Classification trees of DWV symptoms. Classification trees to analyze (A) DWV symptom presence (presence = 1, absence = 0) and (B) DWV symptom severity (ordinal scale 0–3, 0 represents no symptoms, and 3 represents debilitating symptoms) were created using bee stock, treatment, log-transformed DWV-A (DWV.A) and DWV-B (DWV.B) levels as predictor variables through recursive partitioning and with a complexity parameter of >0.01 after pruning. The predicted probabilities per decision node are shown beneath each variable categorization (i.e., 0 = asymptomatic, 1 = symptomatic in A; or the scale of 0–3 where 0 represents no symptoms and 3 represents debilitating morphological symptoms based on Supplementary Figure S1). Virus data were based on head tissues (N = 180 bees).

The classification tree for symptom severity (Figure 7B) was more simplified than that for symptom presence. Like the symptom presence tree, the first decision node in the severity tree was the DWV injection treatment group and the second was the level of DWV-B (with the same cut off of 10^10.91). The third node was solely based on honey bee stock, again grouping Italian, Pol-Line, and Saskatraz bees together, indicating that these stocks have an increased probability of having more severe symptoms than Carniolan and Russian stocks. DWV-A levels did not explicitly play a part in the symptom severity tree, but DWV-A levels may have been partially accounted for through the DWV injection treatment as DWV-A levels increased in this treatment group whereas DWV-B levels did not (Table 1; Figure 2).

Using multiple correspondence analysis (MCA) to further validate the relationships among stock, treatment group, and symptom severity (Figure 8A), we found that treatment groups fell along a gradient of symptom severity on dimension 1 (Dim.1: 0.685). DWV injections were represented at the far positive end of Dim.1, and no manipulation controls at the far negative end with mite and PBS injection treatment groups being intermediate along this dimension. Symptom severity (Dim.1: 0.790 and Dim.2: 0.602) loaded onto both dimensions/axes while bee stock loaded primarily on Dim.2 (Dim.1: 0.105, Dim.2: 0.548). Notably, Carniolan, Italian, and Russian stocks loaded on the negative end of Dim.1, while Pol-Line and Saskatraz bees load onto the positive end of Dim.1. Furthermore, the Russian bees were the only stock to load positively onto Dim.2. Overall, Dim.1 explained 15.8% of the variation in the data, while Dim.2 accounted for 12.1%. Dimensions 3–5 also saw significant loadings of the bee stock and treatment group variables.

FIGURE 8

FIGURE 8. Multivariate visualization of bee stock, treatments, and symptom severity. (A) Multiple correspondence analysis to evaluate the interactions among the following categorical variables: bee stocks, treatments, and DWV symptom severity (scale 0–3, 0 represents no symptoms, and 3 represents debilitating morphological symptoms based on Supplementary Figure S1). (B) Factor analysis of mixed data to evaluate the interactions among five bee stocks, four treatment groups, and DWV symptom severity with the addition of log-transformed levels (virus) of both DWV variants (A and B). Both figures represent the correlations between variables, where variables with similar profiles group together. Sample sizes were based on bees with DWV quantification from head tissues (N = 180 bees).

When we used factor analysis of mixed data (FAMD) (Figure 8B) to add DWV-A and DWV-B levels to the multivariate analysis, Dim.1 accounted for 18.5% of the variation in the data, while Dim.2 accounted for 12.1%, with the variables of DWV-A (Dim.1: 30.88 and Dim.2: 6.79) and DWV-B (Dim.1: 5.43 and Dim.2: 35.54) being the primary contributors. Dimensions 3 and 4 had significant loadings of the bee stock (Dim.3: 48.66 and Dim.4: 61.53) and symptom severity (Dim.3: 48.82 and Dim.4: 32.73). In FAMD, the grouping by stock shifted such that Italian bees joined the Pol-Line and Saskatraz grouping with more severe symptoms (positive along Dim.1), while Carniolan and Russian bees still grouped together with less severe symptoms (negative along Dim.1). This grouping of stocks was consistent with the classification tree results for both symptom presence and severity. However, the stocks and treatment groups can also be grouped by DWV levels, where Pol-Line, Russian, and Saskatraz bees (the three mite-resistant stocks within the study) and mite and PBS treatment groups were grouped with DWV-B on the positive end of Dim.2, while Carniolan and Italian bees loaded on the negative end of Dim.2 with DWV-A.

4 Discussion

The goal of this study was to determine if Varroa mite–resistant honey bee stocks (Pol-Line, Russian, and Saskatraz) exhibited differential resistance or tolerance to DWV compared to mite-susceptible stocks (Carniolan and Italian). We exposed white-eyed pupae to mites, injected them with a sublethal dose (10⁵) of DWV obtained from symptomatic adults, injected them with PBS as sham injection, or did not manipulate them at all as a control. After adult emergence, we evaluated dissemination of DWV variants A and B throughout four tissue types and adult emergence time and symptom severity. The DWV variants differed in their dissemination patterns and associations with treatment groups and, potentially, in their contributions to symptom severity. While bee stocks did not differ in overall virus RNA copy levels or dissemination patterns associated with either DWV variant (not indicative of resistance), bee stock interacted with the treatment group and the DWV level to impact symptom severity (indicative of tolerance). However, these stock interactions were not consistently correlated with known mite resistance.

The levels of both DWV variants differed with the treatment group but not in relation to bee stock, indicating that, in this study, stocks were not resistant to DWV unlike the complementary study on injected adults where some stocks appeared to exhibit DWV resistance (Penn et al., 2021). Sublethal pupal injections of DWV (inoculum obtained from symptomatic adults screened using general DWV primers) induced elevated levels of DWV-A but not DWV-B, which is also dissimilar to responses in DWV-injected adults from different cohorts collected within the same timeframe from the same colonies using an aliquot from the same inoculum (both variants were elevated with DWV injections) (Penn et al., 2021). Conversely, PBS injection and mite exposure pupal treatment groups induced higher levels of DWV-B but not DWV-A (Figure 2). The elevated levels of DWV-B in the PBS and mite treatment groups are also similar to prior responses seen in adult bees but not necessarily in pupae (Wu et al., 2017; Dubois et al., 2020; Norton et al., 2020; Penn et al., 2021; Ray et al., 2021). The differential responses to viral infection between pupae and adults have been previously documented both in infection levels and in the ability of dsRNA treatments to mitigate infections (Möckel et al., 2011; Desai et al., 2012; Thaduri et al., 2019). Such ontogenetic responses to viral infections are likely related to differences in immune gene regulation among life stages, with the introduction of infection or stress at certain points such as the pupal stage being more critical due to lower levels of immune response (Wilson-Rich et al., 2008; Laughton et al., 2011; Koleoglu et al., 2017). These results highlight the importance of including life stage as a factor when investigating or selecting for viral resistance.

Despite having no symptoms and no prior exposure to Varroa mites, the untreated control group bees exhibited some DWV infection, possibly from vertical transmission or horizontal transmission in larval food as mites were not present within the pupal cells (Yue et al., 2007; Möckel et al., 2011). Therefore, the increased levels of DWV in the PBS group compared to those of the control group are most likely due to injection trauma–inducing viral replication (Gusachenko et al., 2020; Penn et al., 2021). The Varroa mites used in this study had greater DWV-B levels (10^{4.82 ± 0.20}) than DWV-A (10^{2.42 ± 0.34}), with similar DWV levels observed in previous studies (Gisder et al., 2009; Wu et al., 2017; Tentcheva et al., 2018). However, the observed mite DWV levels were not correlated with the associated bee levels for either DWV variant (unlike Wu et al., 2017). Mites may lose their DWV infection quite rapidly, particularly if mites feed on bees with low virus levels (Posada-Florez et al., 2019); and different mite loads do not necessarily increase bee DWV levels but may be associated with certain DWV variants (Ryabov et al., 2014). Given these data and that the mites were fed on pupae exchanged daily for 24–48 h prior to experimental use, the similar responses to the mite and PBS injection treatment groups in this study are not unprecedented.

Patterns of DWV dissemination to different tissue types following pupal exposure varied with the virus variant and the treatment group but not bee stock, reflecting complementary data on injected adult bees from the same colonies (Penn et al., 2021). For dissemination, DWV infection of both variants was readily evident in bee heads, consistent with work showing high DWV infections within the head particularly in overt infections (Shah et al., 2009; Zioni et al., 2011; Penn et al., 2021). The levels of DWV in other tissue types were variant-dependent. Unless directly injected with DWV, abdomens had lower levels of DWV-A than other tissue types (Figure 3) but had levels of DWV-B matching that of other tissues in the mite and PBS treatment groups (Penn et al., 2021). Treatment group differences in dissemination may indicate that replication of latent viral infections may not only be induced by non-viral stressors but also depend on virus variant and bee immune responses (Anderson and Gibbs, 1988; Yang and Cox-Foster, 2005; Boncristiani et al., 2013; Kuster et al., 2014; Khongphinitbunjong et al., 2015; Annoscia et al., 2019; Remnant et al., 2019; Tehel et al., 2019; Mookhploy et al., 2021). For instance, phenol oxidase and antimicrobial peptides have been shown to be upregulated when bees were injected with DWV but not PBS compared to controls, whereas hemocyte counts increased following both DWV and PBS injections (Ryabov et al., 2016; Millanta et al., 2019).

Both DWV variants contributed to symptom presence – DWV-B levels and DWV injection (associated with elevated DWV-A levels) were significant in the model (Table 5) and the symptom presence/absence classification tree (Figure 7A), similar to prior work (Tehel et al., 2019; Dubois et al., 2020). Symptom severity was linked to all non-control treatment groups in the model (DWV-A associated with DWV injection and DWV-B associated with PBS injection and mite exposure) (Supplementary Table S5). Though previous data on how DWV levels impact symptom severity have had mixed results (Möckel et al., 2011; Desai et al., 2012; Tehel et al., 2019; Dubois et al., 2020; Yañez et al., 2020; Mookhploy et al., 2021), we observed that the DWV treatment group (elevated DWV-A) and higher levels of DWV-B were associated with the most severe symptoms (Figure 7). The two variants may also differ in which symptom severity they cluster with (Figure 8B) where DWV-A (and the DWV treatment group) clusters with the severity of 3 (most severe), whereas DWV-B clusters closer severities 1 and 2. Variant-related differences in mortality (not measured here)—where DWV-A-infected bees had higher mortality than DWV-B infected bees—have also been documented in Australian, DWV-naïve bee populations (Norton et al., 2020). Using the classification tree data (Figure 7), we found that DWV levels of 10^10.91 appeared to be the initial cutoff point for symptom presence and severity. However, these values varied more with the variant in the symptom logit model; symptomatic bees had an average of 10^{9.71 ± 0.20} for DWV-A and 10^{10.14 ± 0.12} for DWV-B compared to asymptotic bees with an average of 10^{8.45 ± 0.14} DWV-A and 10^{9.57 ± 0.11} DWV-B. The two DWV variants were also significantly positively correlated with each other as seen in similar studies of both injected pupae and adults (Dubois et al., 2020; Penn et al., 2021), potentially reflecting the presence of variant recombinants that will need to be addressed in future work (Gisder et al., 2018; Jamnikar-Ciglenecki et al., 2019; Martin and Brettell, 2019).

Bee stock differences could exacerbate symptoms and infection levels after viral exposure as immune responses and rates of both pupal and adult exposure have been shown to shift with host genetics (Santillán-Galicia et al., 2010; Laughton et al., 2011; Möckel et al., 2011; Chang et al., 2021; Hinshaw et al., 2021; Posada-Florez et al., 2021; Weaver et al., 2021). For instance, hygienic bees may have higher exposure rates to DWV than non-hygienic bees since hygienic bees exhibit increased rates of cannibalism on mite-infested pupae (Posada-Florez et al., 2021). When we considered adult emergence time and DWV symptom presence and severity, bee stocks appeared to differ in their tolerance to DWV; though these results were not necessarily consistent with stocks bred for mite resistance (Khongphinitbunjong et al., 2016; Locke et al., 2021). In the multivariate and classification tree analyses (Figures 7, 8), the Russian stock most consistently grouped with Carniolan as being DWV-tolerant (high level of the virus with fewer symptoms and/or lower severity), while Pol-Line grouped with Italian and Saskatraz as more DWV susceptible, particularly when considering symptom severity (Dim.1 axis on Figures 8A,B). However, when focusing on DWV-B levels (Figure 8B), the stock groupings were consistent with mite resistance (positive on Dim.2 axis). Russian, Pol-Line, and Saskatraz grouped near the higher DWV-B loading and the mite and PBS treatment groups on the positive side of Dim.2 in the FAMD (Figure 8B). As prior work has shown that mites are more likely to be infected with DWV-B compared to DWV-A (Gisder and Genersch, 2021), we might expect that mite-resistant stocks would group together when based on variant B.

Differentiation of mite-resistant stocks from each other is not surprising based on some limited prior research (Khongphinitbunjong et al., 2016; Penn et al., 2021) and given their different selection regimens (Rinderer et al., 2010; Robertson et al., 2014; Danka et al., 2016; Saelao et al., 2020). The separation of Pol-Line from Russian and Carniolan bees is similar to the grouping by genetic sequencing conducted by Saelao et al. (2020), where Pol-Line splits into a separate group from Carniolan, Italian, and Russian stocks (Saelao et al., 2020). Although Saelao et al. (2020) did not include Saskatraz bees in their analysis and our Italian bees did not segregate similarly to theirs in our analyses, it appears that the Pol-Line separation from Carniolan and Russian bees may be due, in part, to differential susceptibility to DWV in addition to genetic clustering for other traits (e.g., the strong selection for Varroa sensitive hygienic behavior). When we compared results from this study and the complementary study on injected adult bees from the same colonies (Penn et al., 2021), injected Pol-Line adults exhibited the greatest resistance to DWV-A compared to the other bee stocks but had the lowest tolerance following DWV injections of pupae in this study (Figures 6–8). Russian bees appeared tolerant to DWV overall as injected adults had the highest levels of DWV-A in the other study and the greatest tolerance in this study. These data indicate that there might be tradeoffs in breeding for mite resistance, DWV resistance, and symptom expression or tolerance. Regardless of the particular mechanisms underpinning these differences, for example, physiological responses to DWV or prevention of Varroa parasitism and DWV exposure (Navajas et al., 2008; Khongphinitbunjong et al., 2015; Posada-Florez et al., 2021), stock-related virus resistance or tolerance may help account for differences in colony survival in the field (Locke et al., 2014; Thaduri et al., 2019).

5 Conclusion

These data suggest that the genetic stock of the honey bee host is important for viral tolerance when pupae are exposed to DWV and other stressors but is not always perfectly aligned with breeding for mite resistance. Furthermore, we found that DWV variants differ in their dissemination in newly emerged adults and contribute significantly to DWV symptom presence though differentially to symptom severity. More severe symptoms were associated with DWV injections that elevated DWV-A levels, while less severe symptoms were associated with DWV-B, injection injury, and mite feeding. Given prior work on the same colonies showing that certain stocks of injected adults bees are more resistant to DWV-A infection over time, further study into the timing of infection and stock-mediated immune responses is vital to unraveling how genetic stocks interact with DWV. Comparison of these studies also indicates that more study of isolated viral variants rather than the naturally occurring combination of variants from overtly symptomatic adults is needed to parse out the full physiological impacts of each variant. More broadly, this research continues to highlight the importance of considering the combination of host genotype and pathogen genotype in epidemiological studies.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

MS-F, KH, and YC designed the research. HP performed the experiments, conducted the statistics, analyzed the data, and wrote the manuscript. All authors edited the manuscript, contributed to this article, and approved the submitted version.

Funding

This work is possible from the USDA NIFA Grant 2017-69004-26515 and the USDA ARS.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank Phil Tokaraz for his assistance in creating injection inoculum and for troubleshooting extraction protocols and Lilia de Guzman for her expertise on maintaining the mite treatment group. We would like to thank Evan Bramlet for help in pulling pupae, Hunter Martin for pupal injections, and Natalie Martin and RaeDiance Fuller for tissue dissections, RNA extraction, and RT-PCR. We also thank Albert Robertson for supplying Saskatraz queens. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture. The USDA is an equal opportunity provider and employer.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.909392/full#supplementary-material

References

Anderson, D. L., and Gibbs, A. J. (1988). Inapparent Virus Infections and Their Interactions in Pupae of the Honey Bee (Apis mellifera Linnaeus) in Australia. J. General Virology 69, 1617–1625. doi:10.1099/0022-1317-69-7-1617

CrossRef Full Text | Google Scholar

Annoscia, D., Brown, S. P., Di Prisco, G., De Paoli, E., Del Fabbro, S., Frizzera, D., et al. (2019). Haemolymph Removal by Varroa Mite Destabilizes the Dynamical Interaction between Immune Effectors and Virus in Bees, as Predicted by Volterra's Model. Proc. R. Soc. B 286, 20190331. doi:10.1098/rspb.2019.0331

PubMed Abstract | CrossRef Full Text | Google Scholar

Bąk, B., and Wilde, J. (2015). Grooming Behavior by Worker Bees of Various Subspecies of Honey Bees to removeVarroa Destructormites. J. Apic. Res. 54, 207–215. doi:10.1080/00218839.2016.1147791

CrossRef Full Text | Google Scholar

Barroso-Arévalo, S., Fernández-Carrión, E., Goyache, J., Molero, F., Puerta, F., and Sánchez-Vizcaíno, J. M. (2019). High Load of Deformed Wing Virus and Varroa Destructor Infestation Are Related to Weakness of Honey Bee Colonies in Southern Spain. Front. Microbiol. 10, 1331. doi:10.3389/fmicb.2019.01331

PubMed Abstract | CrossRef Full Text | Google Scholar

Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015). Fitting Linear Mixed-Effects Models Usinglme4. J. Stat. Soft. 67, 1–48. doi:10.18637/jss.v067.i01

CrossRef Full Text | Google Scholar

Benaets, K., Van Geystelen, A., Cardoen, D., De Smet, L., de Graaf, D. C., Schoofs, L., et al. (2017). Covert Deformed Wing Virus Infections Have Long-Term Deleterious Effects on Honeybee Foraging and Survival. Proc. R. Soc. B 284, 20162149. doi:10.1098/rspb.2016.2149

PubMed Abstract | CrossRef Full Text | Google Scholar

Berthoud, H., Imdorf, A., Haueter, M., Radloff, S., and Neumann, P. (2010). Virus Infections and Winter Losses of Honey Bee Colonies (Apis mellifera). J. Apic. Res. 49, 60–65. doi:10.3896/IBRA.1.49.1.08

CrossRef Full Text | Google Scholar

Boncristiani, H. F., Di Prisco, G., Pettis, J. S., Hamilton, M., and Chen, Y. P. (2009). Molecular Approaches to the Analysis of Deformed Wing Virus Replication and Pathogenesis in the Honey Bee, Apis mellifera. Virol. J. 6, 221. doi:10.1186/1743-422X-6-221

PubMed Abstract | CrossRef Full Text | Google Scholar

Boncristiani, H. F., Evans, J. D., Chen, Y., Pettis, J., Murphy, C., Lopez, D. L., et al. (2013). In Vitro infection of Pupae with Israeli Acute Paralysis Virus Suggests Disturbance of Transcriptional Homeostasis in Honey Bees (Apis mellifera). PLOS ONE 8, e73429. doi:10.1371/journal.pone.0073429

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouuaert, D. C., De Smet, L., and de Graaf, D. C. (2021). Breeding for Virus Resistance and its Effects on Deformed Wing Virus Infection Patterns in Honey Bee Queens. Viruses 13, 1074. doi:10.3390/v13061074

PubMed Abstract | CrossRef Full Text | Google Scholar

Brettell, L. E., Mordecai, G. J., Schroeder, D. C., Jones, I. M., da Silva, J. R., Vicente-Rubiano, M., et al. (2017). A Comparison of Deformed Wing Virus in Deformed and Asymptomatic Honey Bees. Insects 8, 28. doi:10.3390/insects8010028

PubMed Abstract | CrossRef Full Text | Google Scholar

BüChler, R. (1994). Varroa Tolerance in Honey Bees-Occurrence, Characters and Breeding. Bee World 75, 54–70. doi:10.1080/0005772X.1994.11099201

CrossRef Full Text | Google Scholar

Burgan, S. C., Gervasi, S. S., Johnson, L. R., and Martin, L. B. (2019). How Individual Variation in Host Tolerance Affects Competence to Transmit Parasites. Physiological Biochem. Zoology 92, 49–57. doi:10.1086/701169

PubMed Abstract | CrossRef Full Text | Google Scholar

Caron, D. M. (2021). “Honey Bee Strains,” in Honey Bee Medicine for the Veterinary Practitioner (Hoboken, New Jersey: John Wiley & Sons), 73–80. doi:10.1002/9781119583417.ch6

CrossRef Full Text | Google Scholar

Chang, Z.-T., Huang, Y.-F., Chen, Y.-W., Yen, M.-R., Hsu, P.-Y., Chen, T.-H., et al. (2021). Transcriptome-level Assessment of the Impact of Deformed Wing Virus on Honey Bee Larvae. Sci. Rep. 11, 15028. doi:10.1038/s41598-021-94641-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Evans, J., and Feldlaufer, M. (2006a). Horizontal and Vertical Transmission of Viruses in the Honey Bee, Apis mellifera. J. Invertebr. Pathology 92, 152–159. doi:10.1016/j.jip.2006.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y. P., Pettis, J. S., Collins, A., and Feldlaufer, M. F. (2006b). Prevalence and Transmission of Honeybee Viruses. Appl. Environ. Microbiol. 72, 606–611. doi:10.1128/AEM.72.1.606-611.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y. P., and Siede, R. (2007). Honey Bee Viruses. Adv. Virus Res. 70, 33–80. doi:10.1016/S0065-3527(07)70002-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Corby-Harris, V., and Snyder, L. A. (2018). Measuring Hypopharyngeal Gland Acinus Size in Honey Bee (Apis mellifera) Workers. JoVE 139, e58261. doi:10.3791/58261

PubMed Abstract | CrossRef Full Text | Google Scholar

Dainat, B., Evans, J. D., Chen, Y. P., Gauthier, L., and Neumann, P. (2012). Dead or Alive: Deformed Wing Virus and Varroa Destructor Reduce the Life Span of Winter Honeybees. Appl. Environ. Microbiol. 78, 981–987. doi:10.1128/AEM.06537-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Danka, R. G., Harris, J. W., and Dodds, G. E. (2016). Selection of VSH-Derived "Pol-Line" Honey Bees and Evaluation of Their Varroa-Resistance Characteristics. Apidologie 47, 483–490. doi:10.1007/s13592-015-0413-7

CrossRef Full Text | Google Scholar

Danka, R. G., Harris, J. W., Villa, J. D., and Dodds, G. E. (2013). Varying Congruence of Hygienic Responses to Varroa Destructor and Freeze-Killed Brood Among Different Types of Honeybees. Apidologie 44, 447–457. doi:10.1007/s13592-013-0195-8

CrossRef Full Text | Google Scholar

de Guzman, L., Frake, A., and Simone-Finstrom, M. (2017). Comparative Flight Activities and Pathogen Load of Two Stocks of Honey Bees Reared in Gamma-Irradiated Combs. Insects 8, 127. doi:10.3390/insects8040127

CrossRef Full Text | Google Scholar

de Guzman, L. I., Rinderer, T. E., and Frake, A. M. (2007). Growth of Varroa Destructor (Acari: Varroidae) Populations in Russian Honey Bee (Hymenoptera: Apidae) Colonies. Ann. Entomol. Soc. Am. 100, 187–195. doi:10.1603/0013-8746(2007)100[187:GOVDAV]2.0.CO;2

CrossRef Full Text | Google Scholar

de Guzman, L. I., Simone-Finstrom, M., Frake, A. M., and Tokarz, P. (2019). Comb Irradiation Has Limited, Interactive Effects on Colony Performance or Pathogens in Bees, Varroa Destructor and Wax Based on Two Honey Bee Stocks. Insects 10, 15. doi:10.3390/insects10010015

CrossRef Full Text | Google Scholar

de Miranda, J. R., and Genersch, E. (2010). Deformed Wing Virus. J. Invertebr. Pathology 103, S48–S61. doi:10.1016/j.jip.2009.06.012

CrossRef Full Text | Google Scholar

Desai, S. D., Eu, Y.-J., Whyard, S., and Currie, R. W. (2012). Reduction in Deformed Wing Virus Infection in Larval and Adult Honey Bees (Apis mellifera L.) by Double-Stranded RNA Ingestion. Insect Mol. Biol. 21, 446–455. doi:10.1111/j.1365-2583.2012.01150.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Diagne, C. T., Diallo, D., Faye, O., Ba, Y., Faye, O., Gaye, A., et al. (2015). Potential of Selected Senegalese Aedes Spp. Mosquitoes (Diptera: Culicidae) to Transmit Zika Virus. BMC Infect. Dis. 15, 492. doi:10.1186/s12879-015-1231-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Dietemann, V., Nazzi, F., Martin, S. J., Anderson, D. L., Locke, B., Delaplane, K. S., et al. (2013). Standard Methods for Varroa Research. J. Apic. Res. 52, 1–54. doi:10.3896/IBRA.1.52.1.09

CrossRef Full Text | Google Scholar

Donzé, G., and Guerin, P. M. (1994). Behavioral Attributes and Parental Care of Varroa Mites Parasitizing Honeybee Brood. Behav. Ecol. Sociobiol. 34, 305–319. doi:10.1007/BF00197001

CrossRef Full Text | Google Scholar

Dubois, E., Dardouri, M., Schurr, F., Cougoule, N., Sircoulomb, F., and Thiéry, R. (2020). Outcomes of Honeybee Pupae Inoculated with Deformed Wing Virus Genotypes A and B. Apidologie 51, 18–34. doi:10.1007/s13592-019-00701-z

CrossRef Full Text | Google Scholar

Fox, J. (2003). Effect Displays inRfor Generalised Linear Models. J. Stat. Soft. 8, 1–27. doi:10.18637/jss.v008.i15

CrossRef Full Text | Google Scholar

Fox, J., and Hong, J. (2009). Effect Displays inRfor Multinomial and Proportional-Odds Logit Models: Extensions to theeffectsPackage. J. Stat. Soft. 32, 1–24. doi:10.18637/jss.v032.i01

CrossRef Full Text | Google Scholar

Fox, J., and Weisberg, S. (2019). An R Companion to Applied Regression. 3rd edition. Thousand Oaks, CA: SAGE. Available at: https://socialsciences.mcmaster.ca/jfox/Books/Companion/index.html.

Google Scholar

Gisder, S., Aumeier, P., and Genersch, E. (2009). Deformed Wing Virus: Replication and Viral Load in Mites (Varroa Destructor). J. Gen. Virol. 90, 463–467. doi:10.1099/vir.0.005579-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Gisder, S., and Genersch, E. (2021). Direct Evidence for Infection of Varroa Destructor Mites with the Bee-Pathogenic Deformed Wing Virus Variant B, but Not Variant A, via Fluorescence In Situ Hybridization Analysis. J. Virol. 95, e01786–20. doi:10.1128/JVI.01786-20

CrossRef Full Text | Google Scholar

Gisder, S., Möckel, N., Eisenhardt, D., and Genersch, E. (2018). In Vivoevolution of Viral Virulence: Switching of Deformed Wing Virus between Hosts Results in Virulence Changes and Sequence Shifts. Environ. Microbiol. 20, 4612–4628. doi:10.1111/1462-2920.14481

PubMed Abstract | CrossRef Full Text | Google Scholar

Grindrod, I., Kevill, J. L., Villalobos, E. M., Schroeder, D. C., and Martin, S. J. (2021). Ten Years of Deformed Wing Virus (DWV) in Hawaiian Honey Bees (Apis mellifera), the Dominant DWV-A Variant Is Potentially Being Replaced by Variants with a DWV-B Coding Sequence. Viruses 13, 969. doi:10.3390/v13060969

PubMed Abstract | CrossRef Full Text | Google Scholar

Guichard, M., Dietemann, V., Neuditschko, M., and Dainat, B. (2020). Advances and Perspectives in Selecting Resistance Traits against the Parasitic Mite Varroa Destructor in Honey Bees. Genet. Sel. Evol. 52, 71. doi:10.1186/s12711-020-00591-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Gusachenko, O. N., Woodford, L., Balbirnie-Cumming, K., Campbell, E. M., Christie, C. R., Bowman, A. S., et al. (2020). Green Bees: Reverse Genetic Analysis of Deformed Wing Virus Transmission, Replication, and Tropism. Viruses 12, 532. doi:10.3390/v12050532

PubMed Abstract | CrossRef Full Text | Google Scholar

Guzman-Novoa, E., Emsen, B., Unger, P., Espinosa-Montaño, L. G., and Petukhova, T. (2012). Genotypic Variability and Relationships between Mite Infestation Levels, Mite Damage, Grooming Intensity, and Removal of Varroa Destructor Mites in Selected Strains of Worker Honey Bees (Apis mellifera L.). J. Invertebr. Pathology 110, 314–320. doi:10.1016/j.jip.2012.03.020

CrossRef Full Text | Google Scholar

Harbo, J. R., and Harris, J. W. (1999a). Heritability in Honey Bees (Hymenoptera: Apidae) of Characteristics Associated with Resistance to Varroa jacobsoni(Mesostigmata: Varroidae). J. Econ. Entomol. 92, 261–265. doi:10.1093/jee/92.2.261

CrossRef Full Text | Google Scholar

Harbo, J. R., and Harris, J. W. (1999b). Selecting Honey Bees for Resistance to Varroa Jacobsoni. Apidologie 30, 183–196. doi:10.1051/apido:19990208

CrossRef Full Text | Google Scholar

Harris, J. W. (2007). Bees with Varroa Sensitive Hygiene Preferentially Remove Mite Infested Pupae Aged ≤ Five Days Post Capping. J. Api. Res. 46, 134–139. doi:10.1080/00218839.2007.1110138310.3896/ibra.1.46.3.02

CrossRef Full Text | Google Scholar

Hinshaw, C., Evans, K. C., Rosa, C., and López-Uribe, M. M. (2021). The Role of Pathogen Dynamics and Immune Gene Expression in the Survival of Feral Honey Bees. Front. Ecol. Evol. 8, 505. doi:10.3389/fevo.2020.594263

CrossRef Full Text | Google Scholar

Ifantidis, M. D. (1983). Ontogenesis of the MiteVarroa Jacobsoniin Worker and Drone Honeybee Brood Cells. J. Apic. Res. 22, 200–206. doi:10.1080/00218839.1983.11100588

CrossRef Full Text | Google Scholar

Ifantidis, M. D. (1984). Parameters of the Population Dynamics of theVarroaMite on Honeybees. J. Apic. Res. 23, 227–233. doi:10.1080/00218839.1984.11100637

CrossRef Full Text | Google Scholar

Iqbal, J., and Mueller, U. (2007). Virus Infection Causes Specific Learning Deficits in Honeybee Foragers. Proc. R. Soc. B 274, 1517–1521. doi:10.1098/rspb.2007.0022

PubMed Abstract | CrossRef Full Text | Google Scholar

Jamnikar-Ciglenecki, U., Pislak Ocepek, M., and Toplak, I. (2019). Genetic Diversity of Deformed Wing Virus from Apis mellifera Carnica (Hymenoptera: Apidae) and Varroa Mite (Mesostigmata: Varroidae). J. Econ. Entomol. 112, 11–19. doi:10.1093/jee/toy312

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, S., Robertson, T., Mostajeran, M., Robertson, A. J., and Qiu, X. (2016). Differential Gene Expression of Two Extreme Honey Bee (Apis mellifera) Colonies Showing Varroa Tolerance and Susceptibility. Insect Mol. Biol. 25, 272–282. doi:10.1111/imb.12217

PubMed Abstract | CrossRef Full Text | Google Scholar

Kevill, J. L., de Souza, F. S., Sharples, C., Oliver, R., Schroeder, D. C., and Martin, S. J. (2019). DWV-A Lethal to Honey Bees (Apis mellifera): A Colony Level Survey of DWV Variants (A, B, and C) in England, Wales, and 32 States across the US. Viruses 11, 426. doi:10.3390/v11050426

PubMed Abstract | CrossRef Full Text | Google Scholar

Khongphinitbunjong, K., de Guzman, L. I., Rinderer, T. E., Tarver, M. R., Frake, A. M., Chen, Y., et al. (2016). Responses of Varroa -resistant Honey Bees ( Apis mellifera L.) to Deformed Wing Virus. J. Asia-Pacific Entomology 19, 921–927. doi:10.1016/j.aspen.2016.08.008

CrossRef Full Text | Google Scholar

Khongphinitbunjong, K., de Guzman, L. I., Tarver, M. R., Rinderer, T. E., Chen, Y., and Chantawannakul, P. (2015). Differential Viral Levels and Immune Gene Expression in Three Stocks of Apis mellifera Induced by Different Numbers of Varroa Destructor. J. Insect Physiology 72, 28–34. doi:10.1016/j.jinsphys.2014.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Koleoglu, G., Goodwin, P. H., Reyes-Quintana, M., Hamiduzzaman, M. M., and Guzman-Novoa, E. (2017). Effect of Varroa Destructor, Wounding and Varroa Homogenate on Gene Expression in Brood and Adult Honey Bees. PLOS ONE 12, e0169669. doi:10.1371/journal.pone.0169669

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuster, R. D., Boncristiani, H. F., and Rueppell, O. (2014). Immunogene and Viral Transcript Dynamics during Parasitic Varroa Destructor Mite Infection of Developing Honey Bee (Apis mellifera) Pupae. J. Exp. Biol. 217, 1710–1718. doi:10.1242/jeb.097766

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuznetsova, A., Brockhoff, P. B., and Christensen, R. H. B. (2017). Lmertest Package: Tests in Linear Mixed Effects Models. J. Stat. Soft. 82, 1–26. doi:10.18637/jss.v082.i13

CrossRef Full Text | Google Scholar

Laughton, A. M., Boots, M., and Siva-Jothy, M. T. (2011). The Ontogeny of Immunity in the Honey Bee, Apis mellifera L. Following an Immune Challenge. J. Insect Physiology 57, 1023–1032. doi:10.1016/j.jinsphys.2011.04.020

CrossRef Full Text | Google Scholar

Le Conte, Y., Ellis, M., and Ritter, W. (2010). Varroamites and Honey Bee Health: canVarroaexplain Part of the Colony Losses? Apidologie 41, 353–363. doi:10.1051/apido/2010017

CrossRef Full Text | Google Scholar

Lê, S., Josse, J., and Husson, F. (2008). FactoMineR: AnRPackage for Multivariate Analysis. J. Stat. Soft. 25, 1–18. doi:10.18637/jss.v025.i01

CrossRef Full Text | Google Scholar

Lenth, R. V. (2020). Emmeans: Estimated Marginal Means, Aka Least-Squares Means. Available at: https://CRAN.R-project.org/package=emmeans (Accessed April 16, 2021).

Google Scholar

Locke, B., Forsgren, E., and de Miranda, J. R. (2014). Increased Tolerance and Resistance to Virus Infections: A Possible Factor in the Survival of Varroa Destructor-resistant Honey Bees (Apis mellifera). PLOS ONE 9, e99998. doi:10.1371/journal.pone.0099998

PubMed Abstract | CrossRef Full Text | Google Scholar

Locke, B., Thaduri, S., Stephan, J. G., Low, M., Blacquière, T., Dahle, B., et al. (2021). Adapted Tolerance to Virus Infections in Four Geographically Distinct Varroa Destructor-resistant Honeybee Populations. Sci. Rep. 11, 12359. doi:10.1038/s41598-021-91686-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Loope, K. J., Baty, J. W., Lester, P. J., and Wilson Rankin, E. E. (2019). Pathogen Shifts in a Honeybee Predator Following the Arrival of the Varroa Mite. Proc. R. Soc. B 286, 20182499. doi:10.1098/rspb.2018.2499

PubMed Abstract | CrossRef Full Text | Google Scholar

Macedo, P. A., Wu, J., and Ellis, M. D. (2002). Using Inert Dusts to Detect and Assess Varroa Infestations in Honey Bee Colonies. J. Apic. Res. 41, 3–7. doi:10.1080/00218839.2002.11101062

CrossRef Full Text | Google Scholar

Martin, S. J., and Brettell, L. E. (2019). Deformed Wing Virus in Honeybees and Other Insects. Annu. Rev. Virol. 6, 49–69. doi:10.1146/annurev-virology-092818-015700

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, S. J., Highfield, A. C., Brettell, L., Villalobos, E. M., Budge, G. E., Powell, M., et al. (2012). Global Honey Bee Viral Landscape Altered by a Parasitic Mite. Science 336, 1304–1306. doi:10.1126/science.1220941

PubMed Abstract | CrossRef Full Text | Google Scholar

McMahon, D. P., Natsopoulou, M. E., Doublet, V., Fürst, M., Weging, S., Brown, M. J. F., et al. (2016). Elevated Virulence of an Emerging Viral Genotype as a Driver of Honeybee Loss. Proc. R. Soc. B 283, 20160811. doi:10.1098/rspb.2016.0811

PubMed Abstract | CrossRef Full Text | Google Scholar

Mendoza, Y., Tomasco, I. H., Antúnez, K., Castelli, L., Branchiccela, B., Santos, E., et al. (2020). Unraveling Honey Bee-Varroa Destructor Interaction: Multiple Factors Involved in Differential Resistance between Two Uruguayan Populations. Veterinary Sci. 7, 116. doi:10.3390/vetsci7030116

PubMed Abstract | CrossRef Full Text | Google Scholar

Milborrow, S. (2017). Rpart.plot: Plot Rpart Models. An Enhanced Version of plot.Rpart. Available at: http://www.milbo.org/rpart-plot (Accessed November 23, 2021).

Google Scholar

Millanta, F., Sagona, S., Mazzei, M., Forzan, M., Poli, A., Felicioli, A., et al. (2019). Phenoloxidase Activity and Haemolymph Cytology in Honeybees Challenged with a Virus Suspension (Deformed Wings Virus DWV) or Phosphate Buffered Suspension (PBS). Cienc. Rural. 49, e20180726. doi:10.1590/0103-8478cr20180726

CrossRef Full Text | Google Scholar

Mockel, N., Gisder, S., and Genersch, E. (2011). Horizontal Transmission of Deformed Wing Virus: Pathological Consequences in Adult Bees (Apis mellifera) Depend on the Transmission Route. J. General Virology 92, 370–377. doi:10.1099/vir.0.025940-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mondet, F., Alaux, C., Severac, D., Rohmer, M., Mercer, A. R., and Le Conte, Y. (2015). Antennae Hold a Key to Varroa-Sensitive Hygiene Behaviour in Honey Bees. Sci. Rep. 5, 10454. doi:10.1038/srep10454

PubMed Abstract | CrossRef Full Text | Google Scholar

Mondet, F., Beaurepaire, A., McAfee, A., Locke, B., Alaux, C., Blanchard, S., et al. (2020). Honey Bee Survival Mechanisms against the Parasite Varroa Destructor: A Systematic Review of Phenotypic and Genomic Research Efforts. Int. J. Parasitol. 50, 433–447. doi:10.1016/j.ijpara.2020.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Mondet, F., Blanchard, S., Barthes, N., Beslay, D., Bordier, C., Costagliola, G., et al. (2021). Chemical Detection Triggers Honey Bee Defense against a Destructive Parasitic Threat. Nat. Chem. Biol. 17, 524–530. doi:10.1038/s41589-020-00720-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Mondet, F., Kim, S. H., de Miranda, J. R., Beslay, D., Le Conte, Y., and Mercer, A. R. (2016). Specific Cues Associated with Honey Bee Social Defence against Varroa Destructor Infested Brood. Sci. Rep. 6, 25444. doi:10.1038/srep25444

PubMed Abstract | CrossRef Full Text | Google Scholar

Mookhploy, W., Krongdang, S., and Chantawannakul, P. (2021). Effects of Deformed Wing Virus Infection on Expressions of Immune- and Apoptosis-Related Genes in Western Honeybees (Apis mellifera). Insects 12, 82. doi:10.3390/insects12010082

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, J., Jironkin, A., Chandler, D., Burroughs, N., Evans, D. J., and Ryabov, E. V. (2011). Recombinants between Deformed Wing Virus and Varroa Destructor Virus-1 May Prevail in Varroa Destructor-infested Honeybee Colonies. J. General Virology 92, 156–161. doi:10.1099/vir.0.025965-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mordecai, G. J., Brettell, L. E., Martin, S. J., Dixon, D., Jones, I. M., and Schroeder, D. C. (2016a). Superinfection Exclusion and the Long-Term Survival of Honey Bees in Varroa-Infested Colonies. ISME J. 10, 1182–1191. doi:10.1038/ismej.2015.186

PubMed Abstract | CrossRef Full Text | Google Scholar

Mordecai, G. J., Wilfert, L., Martin, S. J., Jones, I. M., and Schroeder, D. C. (2016b). Diversity in a Honey Bee Pathogen: First Report of a Third Master Variant of the Deformed Wing Virus Quasispecies. ISME J. 10, 1264–1273. doi:10.1038/ismej.2015.178

PubMed Abstract | CrossRef Full Text | Google Scholar

Natsopoulou, M. E., McMahon, D. P., Doublet, V., Frey, E., Rosenkranz, P., and Paxton, R. J. (2017). The Virulent, Emerging Genotype B of Deformed Wing Virus Is Closely Linked to Overwinter Honeybee Worker Loss. Sci. Rep. 7, 5242. doi:10.1038/s41598-017-05596-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Navajas, M., Migeon, A., Alaux, C., Martin-Magniette, M., Robinson, G., Evans, J., et al. (2008). Differential Gene Expression of the Honey Bee Apis mellifera Associated with Varroa Destructor Infection. BMC Genomics 9, 301. doi:10.1186/1471-2164-9-301

PubMed Abstract | CrossRef Full Text | Google Scholar

Nazzi, F., and Le Conte, Y. (2016). Ecology of Varroa Destructor, the Major Ectoparasite of the Western Honey Bee, Apis mellifera. Annu. Rev. Entomol. 61, 417–432. doi:10.1146/annurev-ento-010715-023731

PubMed Abstract | CrossRef Full Text | Google Scholar

Nazzi, F., and Milani, N. (1994). A Technique for Reproduction of Varroa Jacobsoni Oud under Laboratory Conditions. Apidologie 25, 579–584. doi:10.1051/apido:19940608

CrossRef Full Text | Google Scholar

Norton, A. M., Remnant, E. J., Buchmann, G., and Beekman, M. (2020). Accumulation and Competition Amongst Deformed Wing Virus Genotypes in Naïve Australian Honeybees Provides Insight into the Increasing Global Prevalence of Genotype B. Front. Microbiol. 11, 620. doi:10.3389/fmicb.2020.00620

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Shea-Wheller, T. A., Rinkevich, F. D., Danka, R. G., Simone-Finstrom, M., Tokarz, P. G., and Healy, K. B. (2022). A Derived Honey Bee Stock Confers Resistance to Varroa Destructor and Associated Viral Transmission. Sci. Rep. 12, 4852. doi:10.1038/s41598-022-08643-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Penn, H. J., Simone-Finstrom, M., Lang, S., Chen, J., and Healy, K. (2021). Host Genotype and Tissue Type Determine DWV Infection Intensity. Front. Insect Sci. 1. doi:10.3389/finsc.2021.756690

CrossRef Full Text | Google Scholar

Piou, V., Tabart, J., Urrutia, V., Hemptinne, J.-L., and Vétillard, A. (2016). Impact of the Phoretic Phase on Reproduction and Damage Caused by Varroa Destructor (Anderson and Trueman) to its Host, the European Honey Bee (Apis mellifera L.). PloS one 11, e0153482. doi:10.1371/journal.pone.0153482

PubMed Abstract | CrossRef Full Text | Google Scholar

Posada-Florez, F., Childers, A. K., Heerman, M. C., Egekwu, N. I., Cook, S. C., Chen, Y., et al. (2019). Deformed Wing Virus Type A, a Major Honey Bee Pathogen, Is Vectored by the Mite Varroa Destructor in a Non-propagative Manner. Sci. Rep. 9, 12445. doi:10.1038/s41598-019-47447-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Posada-Florez, F., Lamas, Z. S., Hawthorne, D. J., Chen, Y., Evans, J. D., and Ryabov, E. V. (2021). Pupal Cannibalism by Worker Honey Bees Contributes to the Spread of Deformed Wing Virus. Sci. Rep. 11, 8989. doi:10.1038/s41598-021-88649-y

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team (2020). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Available at: http://www.R-project.org/.

Google Scholar

Ray, A. M., Davis, S. L., Rasgon, J. L., and Grozinger, C. M. (2021). Simulated Vector Transmission Differentially Influences Dynamics of Two Viral Variants of Deformed Wing Virus in Honey Bees (Apis mellifera). J. Gen. Virol. 102, 001687. doi:10.1099/jgv.0.001687

CrossRef Full Text | Google Scholar

Remnant, E. J., Mather, N., Gillard, T. L., Yagound, B., and Beekman, M. (2019). Direct Transmission by Injection Affects Competition Among RNA Viruses in Honeybees. Proc. R. Soc. B 286, 20182452. doi:10.1098/rspb.2018.2452

PubMed Abstract | CrossRef Full Text | Google Scholar

Rinderer, T. E., Harris, J. W., Hunt, G. J., and de Guzman, L. I. (2010). Breeding for Resistance toVarroa Destructorin North America. Apidologie 41, 409–424. doi:10.1051/apido/2010015

CrossRef Full Text | Google Scholar

Ripley, B., Venables, B., Bates, D. M., Hornik, K., Gebhardt, A., and Firth, D. (2018). MASS: Support Functions and Datasets for Venables and Ripley’s MASS. Available at: https://CRAN.R-project.org/package=MASS (Accessed February 18, 2019).

Google Scholar

Robertson, A. J., Trost, B., Scruten, E., Robertson, T., Mostajeran, M., Connor, W., et al. (2014). Identification of Developmentally-specific Kinotypes and Mechanisms of Varroa Mite Resistance through Whole-Organism, Kinome Analysis of Honeybee. Front. Genet. 5, 139. doi:10.3389/fgene.2014.00139

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenkranz, P., Aumeier, P., and Ziegelmann, B. (2010). Biology and Control of Varroa Destructor. J. Invertebr. Pathology 103, S96–S119. doi:10.1016/j.jip.2009.07.016

CrossRef Full Text | Google Scholar

Rothenbuhler, W. C. (1964). Behavior Genetics of Nest Gleaning in Honey Bees. IV. REsponses of F1 and Backcross Generations to Disease-Killed Brood. Am. Zool. 4, 111–123. doi:10.1093/icb/4.2.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryabov, E. V., Childers, A. K., Chen, Y., Madella, S., Nessa, A., vanEngelsdorp, D., et al. (2017). Recent Spread of Varroa Destructor Virus-1, a Honey Bee Pathogen, in the United States. Sci. Rep. 7, 17447. doi:10.1038/s41598-017-17802-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryabov, E. V., Fannon, J. M., Moore, J. D., Wood, G. R., Evans, D. J., Moore, J. D., et al. (2016). The Iflaviruses Sacbrood Virus and Deformed Wing Virus Evoke Different Transcriptional Responses in the Honeybee Which May Facilitate Their Horizontal or Vertical Transmission. PeerJ 4, e1591. doi:10.7717/peerj.1591

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryabov, E. V., Wood, G. R., Fannon, J. M., Moore, J. D., Bull, J. C., Chandler, D., et al. (2014). A Virulent Strain of Deformed Wing Virus (DWV) of Honeybees (Apis mellifera) Prevails after Varroa Destructor-mediated, or In Vitro, Transmission. PLOS Pathog. 10, e1004230. doi:10.1371/journal.ppat.1004230

PubMed Abstract | CrossRef Full Text | Google Scholar

Saelao, P., Simone-Finstrom, M., Avalos, A., Bilodeau, L., Danka, R., de Guzman, L., et al. (2020). Genome-wide Patterns of Differentiation within and Among U.S. Commercial Honey Bee Stocks. BMC Genomics 21, 704. doi:10.1186/s12864-020-07111-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Santillán-Galicia, M. T., Ball, B. V., Clark, S. J., and Alderson, P. G. (2010). Transmission of Deformed Wing Virus and Slow Paralysis Virus to Adult Bees (Apis melliferaL.) byVarroa Destructor. J. Apic. Res. 49, 141–148. doi:10.3896/IBRA.1.49.2.01

CrossRef Full Text | Google Scholar

Schurr, F., Tison, A., Militano, L., Cheviron, N., Sircoulomb, F., Rivière, M.-P., et al. (2019). Validation of Quantitative Real-Time RT-PCR Assays for the Detection of Six Honeybee Viruses. J. Virological Methods 270, 70–78. doi:10.1016/j.jviromet.2019.04.020

CrossRef Full Text | Google Scholar

Shah, K. S., Evans, E. C., and Pizzorno, M. C. (2009). Localization of Deformed Wing Virus (DWV) in the Brains of the Honeybee, Apis mellifera Linnaeus. Virol. J. 6, 182. doi:10.1186/1743-422X-6-182

PubMed Abstract | CrossRef Full Text | Google Scholar

Simone-Finstrom, M., Aronstein, K., Goblirsch, M., Rinkevich, F., and de Guzman, L. (2018). Gamma Irradiation Inactivates Honey Bee Fungal, Microsporidian, and Viral Pathogens and Parasites. J. Invertebr. Pathology 153, 57–64. doi:10.1016/j.jip.2018.02.011

CrossRef Full Text | Google Scholar

Spivak, M., and Danka, R. G. (2021). Perspectives on Hygienic Behavior in Apis mellifera and Other Social Insects. Apidologie 52, 1–16. doi:10.1007/s13592-020-00784-z

CrossRef Full Text | Google Scholar

Spivak, M., and Reuter, G. S. (2001a). Resistance to American Foulbrood Disease by Honey Bee Colonies Apis mellifera Bred for Hygienic Behavior. Apidologie 32, 555–565. doi:10.1051/apido:2001103

CrossRef Full Text | Google Scholar

Spivak, M., and Reuter, G. S. (2001b). Varroa Destructor Infestation in Untreated Honey Bee (Hymenoptera: Apidae) Colonies Selected for Hygienic Behavior. ec 94, 326–331. doi:10.1603/0022-0493-94.2.326

PubMed Abstract | CrossRef Full Text | Google Scholar

Strauss, U., Human, H., Gauthier, L., Crewe, R. M., Dietemann, V., and Pirk, C. W. W. (2013). Seasonal Prevalence of Pathogens and Parasites in the Savannah Honeybee (Apis mellifera Scutellata). J. Invertebr. Pathology 114, 45–52. doi:10.1016/j.jip.2013.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Tehel, A., Vu, Q., Bigot, D., Gogol-Döring, A., Koch, P., Jenkins, C., et al. (2019). The Two Prevalent Genotypes of an Emerging Infectious Disease, Deformed Wing Virus, Cause Equally Low Pupal Mortality and Equally High Wing Deformities in Host Honey Bees. Viruses 11, 114. doi:10.3390/v11020114

PubMed Abstract | CrossRef Full Text | Google Scholar

Tentcheva, D., Gauthier, L., Bagny, L., Fievet, J., Dainat, B., Cousserans, F., et al. (2006). Comparative Analysis of Deformed Wing Virus (DWV) RNA inApis melliferaandVarroa Destructor. Apidologie 37, 41–50. doi:10.1051/apido:2005057

CrossRef Full Text | Google Scholar

Thaduri, S., Locke, B., Granberg, F., and de Miranda, J. R. (2018). Temporal Changes in the Viromes of Swedish Varroa-Resistant and Varroa-Susceptible Honeybee Populations. PLOS ONE 13, e0206938. doi:10.1371/journal.pone.0206938

PubMed Abstract | CrossRef Full Text | Google Scholar

Thaduri, S., Stephan, J. G., de Miranda, J. R., and Locke, B. (2019). Disentangling Host-Parasite-Pathogen Interactions in a Varroa-Resistant Honeybee Population Reveals Virus Tolerance as an Independent, Naturally Adapted Survival Mechanism. Sci. Rep. 9, 6221. doi:10.1038/s41598-019-42741-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Therneau, T., Atkinson, B., and Ripley, B. (2018). Package “rpart.”. Available at: https://cran.r-project.org/web/packages/rpart/rpart.pdf (Accessed May 3, 2018).

Google Scholar

Wagoner, K., Spivak, M., Hefetz, A., Reams, T., and Rueppell, O. (2019). Stock-specific Chemical Brood Signals Are Induced by Varroa and Deformed Wing Virus, and Elicit Hygienic Response in the Honey Bee. Sci. Rep. 9, 1–14. doi:10.1038/s41598-019-45008-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J. (2016). Identification of Potential Biomarker Genes for Selecting Varroa Tolerant Honey Bees (Apis mellifera) and Biochemical Characterization of a Differentially Expressed Carboxylesterase Gene in Response to Mite Infestation (Saskatoon (Saskatchewan, Canada): University of Saskatchewan). master’s thesis.

Weaver, D. B., Cantarel, B. L., Elsik, C. G., Boncristiani, D. L., and Evans, J. D. (2021). Multi-tiered Analyses of Honey Bees that Resist or Succumb to Parasitic Mites and Viruses. BMC Genomics 22, 720. doi:10.1186/s12864-021-08032-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wells, T., Wolf, S., Nicholls, E., Groll, H., Lim, K. S., Clark, S. J., et al. (2016). Flight Performance of Actively Foraging Honey Bees Is Reduced by a Common Pathogen. Environ. Microbiol. Rep. 8, 728–737. doi:10.1111/1758-2229.12434

PubMed Abstract | CrossRef Full Text | Google Scholar

Wickam, H. (2016). ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag. Available at: https://ggplot2.tidyverse.org.

Google Scholar

Wilson-Rich, N., Dres, S. T., and Starks, P. T. (2008). The Ontogeny of Immunity: Development of Innate Immune Strength in the Honey Bee (Apis mellifera). J. Insect Physiology 54, 1392–1399. doi:10.1016/j.jinsphys.2008.07.016

CrossRef Full Text | Google Scholar

Wu, Y., Dong, X., and Kadowaki, T. (2017). Characterization of the Copy Number and Variants of Deformed Wing Virus (DWV) in the Pairs of Honey Bee Pupa and Infesting Varroa Destructor or Tropilaelaps Mercedesae. Front. Microbiol. 8, 1558. doi:10.3389/fmicb.2017.01558

PubMed Abstract | CrossRef Full Text | Google Scholar

Yañez, O., Chávez-Galarza, J., Tellgren-Roth, C., Pinto, M. A., Neumann, P., and de Miranda, J. R. (2020). The Honeybee (Apis mellifera) Developmental State Shapes the Genetic Composition of the Deformed Wing Virus-A Quasispecies during Serial Transmission. Sci. Rep. 10, 5956. doi:10.1038/s41598-020-62673-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., and Cox-Foster, D. L. (2005). Impact of an Ectoparasite on the Immunity and Pathology of an Invertebrate: Evidence for Host Immunosuppression and Viral Amplification. Proc. Natl. Acad. Sci. U.S.A. 102, 7470–7475. doi:10.1073/pnas.0501860102

PubMed Abstract | CrossRef Full Text | Google Scholar

Yue, C., and Genersch, E. (2005). RT-PCR Analysis of Deformed Wing Virus in Honeybees (Apis mellifera) and Mites (Varroa Destructor). J. Gen. Virol. 86, 3419–3424. doi:10.1099/vir.0.81401-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Yue, C., Schröder, M., Gisder, S., and Genersch, E. (2007). Vertical-transmission Routes for Deformed Wing Virus of Honeybees (Apis mellifera). J. Gen. Virol. 88, 2329–2336. doi:10.1099/vir.0.83101-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zioni, N., Soroker, V., and Chejanovsky, N. (2011). Replication of Varroa Destructor Virus 1 (VDV-1) and a Varroa Destructor Virus 1-deformed Wing Virus Recombinant (VDV-1-DWV) in the Head of the Honey Bee. Virology 417, 106–112. doi:10.1016/j.virol.2011.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar