Immunogene and viral transcript dynamics during parasitic Varroa destructor mite infection of developing honey bee (Apis mellifera) pupae

Ectoparasites can harm their hosts directly or indirectly by vectoring diseases. The vector-borne diseases and their vectors often share the common interest of attenuating the host's immune response to facilitate feeding and reproduction. Selection for an effective suppression of host immunity presumably has led to the functional integration of polydnaviruses into the genome of parasitoid wasps (Herniou et al., 2013) and immunosuppressive proteins in the saliva of various ectoparasitic arthropods (Wikel, 1999; Zhao et al., 2009). Immunosuppressive effects have also been reported for Varroa destructor Anderson & Trueman 2000, an ectoparasitic mite that feeds on honey bee hemolymph and significantly impacts honey bee health in conjunction with associated viruses (Rosenkranz et al., 2010; Dainat et al., 2012a; Nazzi et al., 2012; Francis et al., 2013).

Novel pathogens and pathogen combinations have been proposed as a potentially important contributor to the recent, dramatic honey bee declines (Evans and Schwarz, 2011; Cornman et al., 2012). Honey bees (Apis mellifera L.) are the most important commercial pollinators worldwide (Potts et al., 2010; Calderone, 2012), and health declines may be due to multiple factors that differ regionally (Ellis et al., 2010; Neumann and Carreck, 2010; Williams et al., 2010; Pettis et al., 2012). However, in Europe and North America, the association between honey bee viruses and V. destructor features prominently among other candidate factors, such as general stress, nutrition and pesticides (Le Conte et al., 2010; Williams et al., 2010).

After its initial transmission to A. mellifera in southeast Asia, V. destructor has followed the distribution of its novel host across the world, resulting in declining natural populations and health problems in managed populations (Sammataro et al., 2000; Rosenkranz et al., 2010). The life cycle of V. destructor involves a phoretic dispersal phase before the female enters a brood cell of a developing honey bee to feed, lay its eggs and raise its offspring. Mite reproduction occurs in close synchrony with the development of its host (Martin, 1994). Shortly after the cell is capped, the female mite wounds the larva and establishes a feeding site that is shared with all subsequent offspring. Within the first 72 h post-capping (HPC), when the host undergoes its final pre-pupal molt, the V. destructor female deposits the first, male egg near the top of the cell (Ifantidis, 1983; Martin, 1994). Subsequently, the female mite deposits diploid eggs that develop into daughters approximately every 30 h. Mating occurs among siblings before emergence of all female mites from the brood cell with their adult host at the end of the bee's developmental period (Martin, 1994).

During their life cycle, V. destructor can effectively vector viruses and other pathogens that reside in the bees' hemolymph within and between honey bee colonies. This vertical transmission drastically changes host–pathogen dynamics and evolution (Martin et al., 2012) and may lead to fatal viral outbreaks (Nazzi et al., 2012; Neumann et al., 2012). At least 18 viruses have been described for A. mellifera, and V. destructor vectors some of the most widespread, such as deformed wing virus (DWV), sacbrood virus and acute bee paralysis virus (Chen and Siede, 2007). In association with V. destructor, the pathogenicity of diverse viruses may increase. For example, the symptoms of DWV infections strongly depend on the co-infection with V. destructor (Yang and Cox-Foster, 2005; Nazzi et al., 2012).

Most putative functions of the honey bee immune system are inferred from studies of homologous genes in Drosophila (Evans et al., 2006). However, several Drosophila immune genes are missing from A. mellifera (Evans et al., 2006), suggesting that some aspects of immunity may not be comparable between these taxa. The functioning of the honey bee immune system in response to viral infection is particularly understudied and no response (Azzami et al., 2012), immunosuppression (Nazzi et al., 2012) or upregulation of immunogenes (Boncristiani et al., 2013) has been reported. Honey bees respond to pattern recognition of bacterial and fungal molecules through the rapid expression of antimicrobial peptides. These inducible immune-effectors include abaecin, hymenoptaecin, apidaecin and defensin, and another inducible effector is prophenoloxidase, which plays a role in melanization (Evans et al., 2006). These are presumably regulated by the main Toll and Imd immune signaling pathways, which rely on several intermediaries, such as relish, for signal transmission and are in turn activated by diverse pattern recognition receptors, such as PGRP-LC, PGRP-S1 and PGRP-S2 (Evans et al., 2006). However, most of these and other immune functions and genes lack direct experimental verification.

Even without a complete functional understanding of the honey bee immune system, expression assays of known immune genes have been used to study the interaction between V. destructor, DWV and honey bees. Varroa parasitism has been suggested to increase the pathogenicity of viruses by suppressing honey bee immunity (Gregory et al., 2005; Yang and Cox-Foster, 2005; Navajas et al., 2008; Nazzi et al., 2012), but the findings are difficult to interpret. Mite presence in naturally infected brood cells correlated with lower transcript abundance of defensin and abaecin, although this effect was only found when few mites were found on a developing bee and not with numerous mites (Gregory et al., 2005). Pupal mite infestation also correlated with diminished responses of multiple immune genes when adult bees were challenged by bacterial injection after emergence (Yang and Cox-Foster, 2005). However, it is unclear whether the results of this study are due to a specific immunosuppression or a general weakening of the bees. Additionally, a delayed suppression of honey bee immunity by V. destructor would not benefit the mites' feeding, which has been suggested as the ultimate reason for the immunosuppression (Gregory et al., 2005). Adult bees from hives with higher mite infestation levels showed decreased levels of the immune system regulator dorsal but the downregulation was associated with infection of DWV instead of V. destructor (Nazzi et al., 2012). The common association between DWV and V. destructor (Yang and Cox-Foster, 2005; Nazzi et al., 2012) may confound effects of V. destructor on honey bee immunity. In addition, an in vitro study of the effects of V. destructor infestation on developing worker bees showed increased expression of several immune genes (Gregorc et al., 2012). Thus, the immune response of the honey bee host to Varroa mite infection warrants further research, particularly in the context of DWV, development and quantitative mite effects.

Here, we describe such a systematic study of quantitative mite effects, simultaneously assessing DWV titers. Specifically, we attempt to capture the temporal dynamics of the three-way interactions between virus, mite and bee. We compare the gene expression of several known immune genes in experimentally mite-infected and control honey bees across the mites' reproductive development in their natural colony context. Furthermore, we describe DWV titers and mite reproductive patterns that are associated with the gene expression patterns. Our results demonstrate a nuanced response to V. destructor infection and demonstrate overall that mite parasitism does not result in a simple downregulation of honey bee immunity.

The manual introduction of adult phoretic mites into newly capped cells of fifth instar honey bee larvae was generally successful, with over 95% of the 740 introduced adult mites surviving until cells were opened again for sampling. Our experimental procedure also led to successful mite reproduction in most cases: for 72 HPC and all later sampling points, 77.5% of all mite-infested cells containing live adults contained at least one offspring. Overall, the number of mite offspring produced in cells increased with time, and cells that were inoculated with one foundress (single mite cells) contained fewer mite offspring than cells that were inoculated with two (double mite cells) or more (3+ mite cells) foundresses (Fig. 1A). In contrast, the number of offspring produced per foundress declined from single mite cells and double mite cells to 3+ mite cells (Fig. 1B).

Overall, the expression for all transcripts was affected by age and treatment, except for abaecin (Table 1). However, treatment and time significantly interacted, except for abaecin and DWV (Table 1), indicating different treatment effects at different time points (Fig. 2).

At 24 HPC, the wounded bees exhibited a higher expression of abaecin than the single-mite infected bees and of PGRP-S1 than the control and 3+ mite infected bees. The expression of apidaecin and relish was significantly higher in the wounded bees than in any other group. The expression of defensin, hymenoptaecin, PPOAct, PRGP-LC and PRGP-S2 was not significantly different among any treatment groups. At 72 HPC, the control group expressed abaecin less than the 3+ mite group and PRGP-S2 less than the wounded and the 3+ mite groups. For apidaecin, defensin, hymenoptaecin and relish, the values of the control group were significantly lower than any other group. PGRP-LC, PGRP-S1 and PPOAct expression levels did not show any significant treatment effects (Fig. 2).

At 120 HPC, the wounded bees expressed significantly more PPOAct than both mite-infected groups and more hymenoptaecin than the control and the single-mite groups. Apidaecin, defensin, PGRP-LC, PRGP-S1, PGRP-S2 and relish were more expressed in the wounded group than in any other group. Abaecin expression was not significantly affected by treatment at this time point. At 192 HPC, the control group showed less expression of apidaecin than the wounded and the single-mite groups, less expression of defensin and relish than the single-mite infected group, and less expression of hymenoptaecin than both mite-infected groups. No significant differences among groups were found at 192 HPC for abaecin, PPOAct, PGRP-LC, PGRP-S1 and PGRP-S2 (Fig. 2).

At 240 HPC, defensin expression was significantly higher in the 3+ mite-infected bees than in any other group. Apidaecin expression was highest in the wounded individuals, followed by the 3+ mite, single-mite and control bees. PGRP-S2 was more expressed in the wounded bees than in either kind of mite-infected bees, and PGRP-LC showed higher expression in the wounded bees than in any other group. Relish was more expressed in the 3+ mite-infected bees than in the control or wounded bees. PGRP-S1 was significantly more expressed in the 3+ mite-infected group than in the control. Abaecin, hymenoptaecin and PPOAct expression were not significantly affected (Fig. 2).

DWV titers were significantly affected by treatment at all sampling time points. At 24 HPC, the DWV titers were lower in the control group than in any other group (Fig. 3). At 72 HPC, DWV titers in the control group were lowest, followed by the wounded, single-mite and 3+ mite-infected groups. At 120 HPC, the control group was significantly lower than the 3+ mite group and at 192 HPC the control group was significantly lower than both mite groups. At 240 HPC, DWV transcripts were least abundant in the control group followed by the wounded, single-mite and 3+ mite-infected groups (Fig. 3).

The comparisons of the raw cycle threshold (C_t) values instead of the ΔC_t values resulted in several different outcomes (supplementary material Fig. S1), particularly with regard to the wounded bees at 24 HPC and 120 HPC. Similar differences among the experimental groups were found for transcript abundance for all genes except relish at 192 HPC and for all but abaecin, PGRP-S2 and relish at 72 HPC. At 24 HPC, outcomes differed between relative and absolute quantification for all genes, at 120 HPC only apidaecin and PPOAct were consistent, and at 240 HPC, abaecin, defensin, PGRP-LC and PGRP-S2 exhibited similar differences among the experimental groups. The significant increases of raw C_t values from control to wounded to the single-mite and 3+ mite-infected bees at all time points were also consistent with the results of the ΔC_t analysis.

Across all experimental groups, the relative gene expression (ΔC_t values) of most immune genes was correlated: 30 of 36 possible pairwise correlations were significant, with correlation coefficients ranging from −0.02 to 0.72. A hierarchical cluster analysis grouped the three pattern recognition genes together, and suggested pairwise connections between defensin and relish and between apidaecin and hymenoptaecin, while abaecin and PPOAct were not included in any particular cluster (Fig. 4). DWV abundance was most strongly correlated to defensin (r_p=0.44, n=209, P<0.001), followed by relish (r_p=0.38, n=209, P<0.001) and hymenoptaecin (r_p=0.32, n=209, P<0.001).

Analyses of the absolute values (raw C_t values) showed that the expression patterns of the two reference genes and all immune genes were correlated with each other (r_p ranging from 0.15 to 0.76, n=209, P<0.05), except for the pair defensin–apidaecin (r_p=−0.04, n=209, P=0.606). The clustering of the genes based on their absolute transcript abundance was identical to the clustering based on ΔC_t values (Fig. 4). DWV titers were correlated to raw C_t values of defensin (r_p=0.27, n=209, P<0.001) and hymenoptaecin (r_p=0.19, n=209, P=0.007), but not to the raw C_t values of the remaining immune genes or reference genes.

Overall, the relative expression levels (ΔC_t) of several immune genes showed significant, negative associations with the number of offspring per mite in the ANCOVA models. The strongest effects were exhibited by relish (B=−0.21, F_1,81=16.2, P<0.001), PGRP-S1 (B=−0.13, F_1,81=9.3, P=0.003) and hymenoptaecin (B=−0.12, F_1,81=8.9, P=0.004), followed by apidaecin (B=−0.09, F_1,81=8.4, P=0.005), defensin (B=−0.09, F_1,81=5.8, P=0.019) and PPOAct(B=−0.06, F_1,81=5.4, P=0.023). Analyzing time points separately, PGRP-S1 (B=−0.07, F_1,18=5.4, P=0.033) and PPOAct (B=−0.04, F_1,18=7.2, P=0.015) showed significant effects at 72 HPC, hymenoptaecin (B=−0.23, F_1,16=10.1, P=0.006) and relish (B=−0.27, F_1,16=6.3, P=0.023) were significant at 192 HPC, and apidaecin (B=−0.24, F_1,22=6.5, P=0.019) showed a significant association with mite reproduction at 240 HPC.

Overall, the reproductive success of introduced mites was significantly associated with the raw C_t values of α-tubulin (B=0.10, F_1,81=7.5, P=0.008) and RPS5 (B=0.09, F_1,81=8.0, P=0.006). When time points were analyzed separately, only RPS5 raw C_t values showed a significant association at 240 HPC (B=0.16, F_1,22=5.4, P=0.030) and PPOAct at 72 HPC (B=−0.04, F_1,18=6.9, P=0.017).

Understanding the interactions between V. destructor and its honey bee hosts is crucial because of the central role of these mites in honey bee health declines (Dainat et al., 2012b; Nazzi et al., 2012; van Dooremalen et al., 2012). Our experimental design allowed us to assess the immediate impact of V. destructor on honey bee immune gene expression over developmental time and in a quantitative fashion. Our findings are presumably applicable to the normal context because the experimentally introduced mites reproduced normally and their honey bee hosts experienced a typical development. Mite parasitism was quantitatively related to titers of DWV in the developing bees (Nazzi et al., 2012). However, our data showed little evidence for immunosuppression of honey bees by V. destructor parasitism, contrary to earlier studies (Gregory et al., 2005; Yang and Cox-Foster, 2005) but supported by more recent analyses (Aronstein et al., 2012; Gregorc et al., 2012). Instead, a complex picture of experimental effects emerged with differences among specific immune genes and time points. The gene expression patterns indicated some functional relationships among the studied immune genes and could be related to the reproductive success of Varroa.

Similar to other studies, our study of gene expression across developmental stages and strong experimental effects experienced the problem of identifying a good normalization method to account for methodological differences among samples. One solution is to search for suitable reference genes to normalize quantitative PCR results until specific reference genes for an experimental paradigm have been identified (Cameron et al., 2013; Reim et al., 2013). However, large sets of candidate genes may have to be studied without a priori information, particularly in two-factorial experimental designs that cross developmental stages (Cameron et al., 2013), such as ours. To reduce normalization bias, multiple reference genes can be combined (Lourenco et al., 2008). We evaluated multiple reference genes that have often been used in the context of honey bee immunity (Boncristiani et al., 2012; Boncristiani et al., 2013) and combined two reference genes based on our preliminary data. However, even this combined reference for calculating ΔC_t values was significantly different among some experimental groups. As an alternative, we also analyzed our data and present the results based on raw C_t values (Boncristiani et al., 2013). However, raw C_t values of all honey bee transcripts were more correlated with each other than ΔC_t values. Therefore, we preferentially present and discuss the relative differences of immune genes but recommend a wider set of reference genes (but see Boncristiani et al., 2013) or specific sets of reference genes (Cameron et al., 2013; Reim et al., 2013) in the future.

Regardless of the method of quantification, the experimental treatments unambiguously and consistently affected the titer of DWV in the developing bees. As reported before (Nazzi et al., 2012), a positive association between V. destructor mites and DWV was found and our results show that the effect is quantitative: more V. destructor mites parasitizing a developing bee translate into higher DWV titers. This pattern seems to be established at 72 HPC and relative and absolute DWV quantities did not increase beyond this time point. The wounding treatment with a sterile needle also increased DWV titers relative to the negative control, but not as much as the V. destructor introductions. It is possible that the open wounds allowed DWV entry into the bee host from the cell environment (Kanbar and Engels, 2003) or that DWV replication was increased by the wounding trauma itself (Boncristiani et al., 2013). This suggests the intriguing possibility that honey bee viruses may be activated by several different stressors (Podgwaite and Mazzone, 1986), including mite feeding. The bacterial cofactor effect reported earlier (Yang and Cox-Foster, 2005) may represent another general stressor or be due to the injection trauma itself. Thus, V. destructor could increase DWV by vectoring (Martin et al., 2012) and activating amplification without necessarily compromising honey bee immunity.

The increased DWV titers with increased foundress number can be explained by either increased vectoring or multiple wounding sites (Donzé and Guerin, 1994). The persistent differences in DWV titers among the experimental groups without a consistent increase of these titers over time suggest that the initial amount of vectored DWV is important in addition to the wounding effect. No evidence for a negative association between immune gene expression and DWV titers was found, in contrast to an earlier study of adult honey bee workers (Yang and Cox-Foster, 2005), suggesting that an active suppression of DWV by the honey bee immune system is unlikely in this case. Instead, a few immune genes were positively associated with the virus titers (Boncristiani et al., 2013) and most did not show a significant association (Azzami et al., 2012). More research is needed to understand virus interactions with insect immunity (Costa et al., 2009; Flenniken and Andino, 2013).

Overall, we found little evidence for a reduced expression of the studied nine immune genes during mite infection and the effects of mite parasitism on the immune system of their honey bee hosts varied among different genes and time points. Mite-infected bees showed a consistently reduced expression only for PGRP-S2 at 240 HPC and PPOAct at 120 HPC. In contrast, higher expression in the mite-infected groups than the control groups was apparent in at least two time points for defensin, hymenoptaecin and relish, and at 72 HPC for abaecin, although not all group differences were significant. The 72 HPC time point is particularly relevant because it coincides with the initiation of mite reproduction (Ifantidis, 1983). In 17 of the 45 studied cases, no significant effect of the different treatments was detected. These results contradict the initial finding of abaecin and defensin suppression by Varroa during bee development (Gregory et al., 2005) and are more consistent with studies reporting an upregulation or no expression change of specific honey bee immune genes in developing brood during Varroa parasitism (Aronstein et al., 2012; Gregorc et al., 2012). Our findings do not exclude the possibility that Varroa parasitism has a delayed effect on adult honey bees (Yang and Cox-Foster, 2005), but such an effect would be more likely mediated by general weakening of the individual (Aronstein et al., 2012) than by specific immunosuppression. Expression differences between the single-mite and 3+ mite groups were significant only for defensin at 240 HPC, indicating that foundress number did not affect host immune gene expression in general.

Even though we have no direct genetic evidence that all introduced mites reproduced, our careful manual introduction of phoretic mites ensured that reproductively active V. destructor females entered the cells of developing honey bee hosts at the correct time. Timing is crucial for the reproductive activity of Varroa (Steiner et al., 1994), and a stringent time window was employed. The reproductive patterns found were comparable to natural situations (Martin, 1994; Rosenkranz et al., 2010), although our numbers may be slight overestimates because juvenile mites were counted the same as mature offspring. Similar to natural situations, we also encountered only a relatively low number of dead or non-reproductive foundresses. The data provide further evidence for competition among multiple foundresses because the individual reproductive success declined with increasing foundress number (Fuchs and Langenbach, 1989; Martin, 1995). However, overall offspring number per cell increased with higher foundress numbers, suggesting that multiple-mite infestations present a stronger drain of energy and nutrients from the developing honey bee.

Mite reproductive success was negatively correlated with the relative expression of most immune genes across all time points. These negative associations were detected despite the positive overall relationships between mites, DWV and transcript levels of immune genes by statistically controlling for treatment and time point effects. Thus, the level of immunity may influence the likelihood of mite reproduction, as has been suggested previously (Gregory et al., 2005). Alternatively, these overall effects could represent artifacts of the normalization because expression of both reference genes based on raw C_t values was positively associated with mite reproduction. However, a plausible explanation for this association does not exist. Moreover, the reference genes cannot statistically account for the association of immune gene expression and mite reproduction at specific time points. These time-specific patterns are intriguing because they suggest an early role for PPOAct, which may be explained by the wound-healing effect of prophenoloxidases that may delay sufficient feeding during the initiation of the feeding site by the foundress (Ifantidis, 1983; Steiner et al., 1994). The expression of the peptidoglycan recognition gene PGRP-S1 was also negatively correlated with mite reproduction at 72 HPC and may be explained by a possible activation of PPOAct by PGRP-S1. The time point of the negative associations of PPOAct and PGRP-S1 with mite reproductive success coincides with the first reproductive event in the mite reproductive phase (Martin, 1994). The effects do not persist, but later, mite reproduction is negatively affected by apidaecin, hymenoptaecin and relish expression, which are downstream members of the Imd immune pathway (Evans et al., 2006). The antimicrobial effector genes may have negative consequences for bacteria that grow in the Varroa feeding sites (Kanbar and Engels, 2003) and thus may help the wound persistence over time. In general, the phenoloxidase function is believed to act immediately while the inducible components of insect immunity may produce longer-lasting effects (Laughton et al., 2011).

We employed a positive and a negative control treatment to compare the effects of V. destructor mite parasitism. Although puncturing the developing honey bee with a small glass capillary has been reported to inflict wounds similar to mite feeding sites (Herrmann et al., 2005), our data show profound, long-lasting perturbations of the expression of the investigated genes. The patterns of ΔC_t values of immune genes relative to each other (supplementary material Fig. S2) suggested that the wounding treatment was qualitatively dissimilar from the negative control and was quite unique at 120 HPC. However, the absolute C_t values of our reference genes were significantly affected by wounding, particularly at 24 HPC and 120 HPC, resulting in increased ΔC_t values in a number of immune genes. We did not visually assess the damage of these wounds but we may have injured the bees more severely than anticipated. The profound effect of wounding on gene expression renders them less suitable as a control. Nevertheless, it shows the importance of physical wounding and may partly explain the differences between single- and multiple-mite-infested honey bees because each foundress mite typically produces a single feeding site (Donzé and Guerin, 1994).

The expression patterns of the investigated immune genes were positively correlated with each other, a finding that was also found in adult bees, although to a more variable extent (Yang and Cox-Foster, 2005). In the absence of much experimental work to support the Drosophila-inferred functional relationships of honey bee immune genes (Evans et al., 2006), co-expression patterns may be useful to infer functional connections among genes (De Gregorio et al., 2001). Consistent with this conjecture, a tight association between the three peptidoglycan recognition genes resulted from our overall clustering of the gene expression patterns. In addition, close pairwise relationships between defensin and relish and between apidaecin and hymenoptaecin were revealed. The former suggests that the honey bee defensin-1 gene depends mostly on the signaling of the Imd pathway (Evans et al., 2006), but a feedback loop from defensin-1 to relish (Erler et al., 2011) might also explain the observed relationship. The other pairing of the AMPs apidaecin and hymenoptaecin suggests that these two genes experience a close coregulation and are probably influenced by additional factors that do not affect defensin-1. Abaecin exhibited the most different expression pattern of the AMPs, suggesting that the regulation of this unique gene (Casteels et al., 1990) is largely independent. Among the studied genes, PPOAct differed most in the expression pattern, which supports the hypothesis that the phenoloxidase function is regulated independently from the inducible parts of the honey bee immune system (Evans et al., 2006; Laughton et al., 2011).

The developmental dynamics of the expression of immune genes relative to each other (supplementary material Fig. S2) illustrates the strong effect of developmental time, indicated by the results of the full-factorial model (supplementary material Table S1). However, our data do not confirm a general increase of humoral immune response genes (Evans, 2004; Laughton et al., 2011). Across the investigated time points, 120 HPC seems to have the lowest overall expression of immunogenes in the negative control and the mite groups. Wounded bees showed a more homogeneous pattern of immune gene expression across developmental time. The 120 HPC time point is close to the peak of metamorphosis, when internal reorganization of tissues, including the fat body, may require available resources, compromising immune gene expression.

All experiments were carried out in three hives of mixed European descent at the University of North Carolina at Greensboro beeyard in Greensboro, NC, USA, between June and August 2011 to capture the peak production of honey bee worker brood. A clear transparency sheet was pinned to experimental frames and used to mark brood cells containing fourth instar honey bee larvae that were old enough for cell capping (De Ruijter, 1987). After 6 h, the experimental frames were re-evaluated to identify larvae that had been sealed off by adult workers in the past 6 h. The cell caps of these larvae were carefully opened to randomly administer one of the following five experimental treatments before closing the cells again. Negative control larvae were left undisturbed. Larvae in a second control group received a small physical wound between the second thoracic segment and the first abdominal segment with a single-use capillary needle pulled to 50 μm diameter (Herrmann et al., 2005). Individuals in the three mite treatment groups received one, two, or three to five (3+) phoretic V. destructor mites. These mites were collected from nurse bees of two independent donor colonies using powdered sugar shakes (Boecking and Ritter, 1993) to ensure that mites were reproductively active foundresses (Kraus, 1993; Steiner et al., 1994). Mites were cleaned of powdered (confectioner's) sugar using a small paintbrush and water and were then transferred into brood cells using ethanol-washed insect pinning needles (Aumeier and Rosenkranz, 2001).

Pseudo-randomized subsets of experimental cells from the five treatment groups were collected at 24 h intervals over the course of 10 days following cell capping. Thus, 10 sampling groups at 24, 48, 72, 96, 120, 144, 168, 192, 216 and 240 HPC with ~10 individuals in each group (N_total=522) were collected (Table 2). The larval or pupal hosts were separated from any associated mites and were placed in 1.5 ml TRIzol™ reagent (Life Technologies, Carlsbad, CA, USA), briefly homogenized with a sterilized pipette tip, and stored on dry ice until they were transferred into an ultralow-temperature freezer (−80°C) to prevent RNA degradation prior to extraction. The number of adult mites and all mite offspring, ranging from egg to pre-adult nymphal stages, was determined visually to determine mite reproductive success (Ifantidis, 1983; Rosenkranz et al., 2010).

Quantitative molecular data were collected from bees of both control and the single mite and 3+ mite treatment groups at the time points of 24, 72, 120, 192 and 240 HPC to capture the beginning, middle and end of host development as well as the initiation and end of reproduction by the foundress mite (Ifantidis, 1983).

RNA extractions were performed following the TRIzol™ protocol, following the manufacturer's recommendations (Life Technologies). Frozen samples were thawed to room temperature, vortexed for 20 to 30 s, and further homogenized using a pipette tip. After adding 300 μl of ultrapure chloroform, each sample was mixed by manual shaking for 15 s, followed by 2–3 min incubation at room temperature and centrifugation at 12,000 g for 15 min at 4°C for phase separation. The uppermost, aqueous phase containing the isolated RNA was carefully transferred into a clean microcentrifuge tube and 750 μl of 100% isopropanol was added. After 10 min of incubation at room temperature, samples were again centrifuged at 12,000 g for 10 min at 4°C. The supernatant was carefully decanted and the RNA pellet was washed with 1.5 ml of 75% ethanol before a final centrifugation at 7500 g for 5 min. The supernatant was removed and once the pellet had sufficiently air-dried, the RNA was resuspended in 50 μl RNase-free water containing 1% RNaseOUT solution (Life Technologies), heating the samples to 60°C for 15 min before storage at −80°C.

After quantification of the RNA using a Nanodrop™ spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), samples were diluted to 500 ng μl⁻¹ using RNase-free water and an 8 μl aliquot was treated with 1 μl DNase I (2 U μl⁻¹; Life Technologies), 1 μl DNase buffer and 0.2 μl RNaseOUT for 1 h at 37°C to remove DNA contaminants and prevent sample degradation. After 10 min at 75°C to inactivate the enzymes, reverse transcription was initiated by mixing the RNA solution with 0.02 μl poly-dT (n=12–18, 0.5 μg μl⁻¹), 0.5 μl random hexamer mix (50 μmol l⁻¹), 0.2 μl 2 mmol l⁻¹ dNTP and 0.298 μl H₂O. The samples were heated to 65°C for 5 min and then placed on ice for 10 min to anneal the random primers. Superscript II (Life Technologies) was used for cDNA synthesis, and 4 μl buffer, 2 μl dithiothreitol, 0.5 μl enzyme (200 U μl⁻¹) and 3.5 μl H₂O were added to each sample and incubated at 42°C for 50 min, then at 70°C for 15 min to inactivate the enzyme. Finally, samples were diluted 1:5 with molecular-grade water (G-Biosciences, St Louis, MO, USA) for a total of 100 μl cDNA template for quantitative PCR (qPCR).

Target sequences for amplification of individual samples included constitutively expressed ‘housekeeping’ genes, multiple genes in the Toll and Imd pathways, antimicrobial peptides, a precursor to ProPO, and DWV (Table 3). These transcripts were selected from a larger list of putative immunological and pathogen targets (Evans, 2006). For these selected targets, a preliminary screening of the published primers followed by melt-curve analysis was performed to ensure specific primer binding and selective amplification of a unique target sequence.

qPCR was performed using Brilliant SYBR Green™ Master Mix (Life Technologies, Carlsbad, CA, USA). Individual reactions contained 10 μl master mix, 1 μl cDNA template, 0.5 μl forward and reverse target primers (8 μmol l⁻¹), and 8 μl molecular-grade water. Reactions were performed in a StepOnePlus™ (Applied Biosystems) thermocycler following established protocols (Evans, 2006), with a slight adjustment to the initial holding stage and annealing temperature. The following cycling conditions were used: 3 min at 95°C, followed by 40 cycles of 95°C for 20 s, 60°C for 30 s, 72°C for 1 min, and 72°C for another 20 s, during which fluorescence measurements were taken. A final melt curve analysis was included at 95°C for 15 s, 60°C for 1 min, and a final ramp at 0.3°C to 95°C for 15 s.

C_t values were collected based on the default StepOnePlus algorithmic threshold search criteria, using a weighted average threshold across the three plates per gene that were necessary to amplify all 209 samples. Reactions that failed to amplify the target transcript were repeated once. When two qPCR runs indicated no detection of the target, a C_t value of 45 was used instead of a missing value to include non-detection in the subsequent analyses. This surrogate value exceeded the maximum number of cycles in the reaction protocol by five to account for the fact that non-detection could be a true zero but also might represent a value just below the detection limit. Initially, raw C_t values determined from qPCR were compared across the three hives, four treatments and five time points with a mixed ANOVA. Treatment and time point effects were abundant but the hive environment did not significantly affect any transcript levels (supplementary material Table S1). Therefore, all subsequent analyses were performed on the three hives combined.

A preliminary evaluation of several potential reference genes, including β-actin, showed inconsistent melt curves or suggested pronounced expression differences between treatment groups or time points. RPS5 and α-tubulin were found to be most suitable, but their full-scale analysis revealed significant differences among treatment groups and time points (supplementary material Table S1). The combined average of the two reference genes also differed among experimental groups (at 24 and 120 HPC the wounded treatment was greater than all others and at 72 and 240 HPC the 3+ treatment was greater than the control). However, stability analysis with Normfinder v0.953 (Andersen et al., 2004) suggested the use of this average as reference instead of either individual reference gene. Therefore, we used the average C_t value of both selected reference genes to calculate ΔC_t values for each of the 10 other transcripts (Antúnez et al., 2009). Without prior knowledge of the amplification efficiency, we report relative results (ΔC_t). In addition, we evaluated the raw C_t values (supplementary material Fig. S1) to provide absolute expression data of the transcripts based on the tentative assumption that cDNA normalization after Nanodrop™ quantification equalized all qPCR templates (Boncristiani et al., 2013). The overall expression levels of the target sequences were tested in response to time and treatment using a full two-factorial ANOVA. To test for treatment effects at specific time points, simple ANOVAs of ΔC_t values were performed, followed by Tukey's HSD post hoc analysis. Pearson's product-moment correlation was performed on the raw C_t data for all 12 transcripts across the full data set to identify potential co-expression patterns.

Additionally, mite reproductive success was compared among the treatment groups and correlated to immune gene expression. Reproductive success was analyzed separately as either total offspring in a cell or a value of individual reproductive success per mite (offspring per foundress mite). All analyses concerning mite reproduction were performed only on mite-infested cells (treatment groups III–V), omitting the first two sampling time points because the initial mite offspring is produced ~60 HPC (Ifantidis, 1983; Martin, 1994). Expression values of target transcripts (ΔC_t) were correlated to mite reproductive success for each time point separately by ANCOVA, taking treatment into account.

We thank Mike Simone-Finstroem for providing advice on some of the experimental procedures. We also acknowledge the support and encouragement of the other members of the Social Insect Lab at UNCG and the North Carolina Honey Bee Research Consortium.