The Role of Pathogen Dynamics and Immune Gene Expression in the Survival of Feral Honey Bees

Introduction

Feralization is the process by which previously domesticated organisms establish populations in the wild in the absence of anthropogenic influence (Gering et al., 2019a). The outcomes of feralization are typically studied in an evolutionary context, examining how environmental and genetic factors affect the fitness of feral organisms compared to their domesticated sources. It has been suggested that because of genetic bottlenecks and artificial selection during the domestication process, fitness decreases outside of captivity, making feral organisms exhibit reduced fitness upon reintroduction to the wild (Araki et al., 2009; Baskett and Waples, 2013; Meyer and Purugganan, 2013). However, feral organisms often thrive (e.g., cats, dogs, pigs) and do not always revert back to the ‘wild-type’ (Taylor et al., 1998; Bellard et al., 2017). Indeed, feral organisms often outnumber their wild counterparts, and can lead to shifts in community composition at the ecosystem level by increasing predation pressure on available prey and potentially spilling over pathogens to wild species (Leiser et al., 2013; Bevins et al., 2014; Maeda et al., 2019; Lepczyk et al., 2020). These ecological effects of feralization are gaining more attention due to their implications for conservation biology and ecosystem management.

Feral organisms frequently interact with both domesticated and wild species, and play a critical ecological role in the dynamics of pathogens shared among these closely related groups. The increased fecundity and expanded geographic ranges that result from the domestication processes, often result in large uncontrolled populations of feral organisms harboring infections and serving as a bridge between domesticated and wild hosts (Bevins et al., 2014). This is the case of the feral swine (Sus scrofa domestica), a species that has high rates of reproduction, high pathogen loads, and overlaps in range with domestic pigs (S. scrofa domestica) and wild boars (S. scrofa) (Taylor et al., 1998; Hill et al., 2014). Transmission of several pathogens, including viruses and bacteria, have been documented from feral swine to domestic pigs and wild boars, posing concerns for the role of feral organisms as reservoirs of pathogen transmission to both wild and domesticated populations (Cvetnic et al., 2003; Le Potier et al., 2006; Meng et al., 2009).

The ecological conditions of feralization may facilitate tolerance or resistance to disease. Pathogens of domesticated species are usually managed by humans to avoid the rapid spread of diseases among domesticated animals. On the contrary, pathogen transmission in feral populations is uncontrolled. Under these conditions, host–pathogen interactions in feral populations may therefore facilitate the rapid evolution of natural mechanisms of disease tolerance or resistance (LeConte et al., 2007; Locke, 2016). Thus, the maintenance of traits associated with disease resilience may be relaxed in domesticated species, while greater ability to mitigate the negative effects of pathogens is critical for the survival of feral organisms (Moreira et al., 2018).

Feral colonies of the honey bee (Apis mellifera L.) provide an ideal model to investigate the hypothesis that host–pathogen dynamics during feralization favor higher expression of defenses and disease tolerance in feral organisms. Apis mellifera is a eusocial bee species that has undergone extensive domestication efforts for traits such as increased honey production, decreased aggressive behavior and reduced frequency of swarming (Lecocq, 2018). Managed honey bee colonies frequently colonize wild environments and become feral because colonies reproduce through swarming (Winston, 1991). Both domesticated and feral honey bees face serious challenges due to a large number of pests and pathogens (Calderone, 2012; Mcmahon et al., 2016). One of the major drivers of disease and colony losses among honey bees is the ectoparasitic mite Varroa destructor, which acts as a vector for multiple bee-infecting RNA viruses that significantly weaken colonies and decrease their overwintering survival (Gisder et al., 2009; Martin et al., 2012). Varroa mites and associated viruses are considered major antagonists of honey bee health and, because of their strong effects on honey bee survival, managed colonies are often treated with chemical acaricides multiple times per year to decrease mite numbers. If untreated, most managed honey bee colonies die within the first year (Kraus and Page, 1995; LeConte et al., 2010). Nevertheless, feral honey bee colonies have been documented as surviving long-term in the wild in the absence of beekeeper management, where Varroa mites and viruses are not artificially controlled, and can therefore pose high selective pressure on colonies (Locke, 2016).

Previous studies have indicated that feral honey bee colonies may exhibit higher immune responses than managed colonies (Youngsteadt et al., 2015, but see Lowe et al., 2011). However, it is unclear how the expression of different immune phenotypes in managed and feral conditions is associated with colony survival and resistance or tolerance to parasites. Honey bees rely on both individual and social mechanisms of immunity to protect the colony from pests and pathogens, and management likely has an influence on both types of defenses (Neumann and Blacquière, 2017; Taric et al., 2020). While behavioral responses, such as hygienic behavior, play key roles in protection against pathogens (Simone-Finstrom, 2017), humoral immune responses in individual bees are also critical for pathogen defense and control of infections (Di Prisco et al., 2016; McMenamin et al., 2018). Several immune pathways are involved in the immune response against viruses, bacteria and fungi. For example, antimicrobial peptides (AMPs) produced from several immune pathways (e.g., Toll and Imd) have key general roles in insect immune systems (Yi et al., 2014; Brutscher et al., 2015). The RNA interference (RNAi) pathway is an antiviral defense pathway in insects that targets sequence-specific double-stranded RNA produced during RNA virus replication (Gammon and Mello, 2015).

Here, we investigate the role of pathogen infections and immune gene expression in the survival of feral and managed honey bees to answer the following questions: (1) are feral colonies reservoirs of pathogens with increased levels of pathogens compared to managed colonies?; (2) do increased pathogen levels lead to higher expression of immune genes in feral colonies than in managed colonies?; (3) is immune gene expression correlated with survival of honey bee colonies? Over a 2-year period, we sampled feral and managed colonies in the same landscapes with the participation of beekeepers who reported the location of colonies. We collected individuals from a total of 44 colonies during the spring and fall of the years 2017 and 2018 to compare the ecoimmunology and survival of colonies in apiaries and in the wild. The virus most directly linked to Varroa mite infestations (DWV) was found at significantly higher levels in feral than managed colonies, and the expression of most immune genes was also higher in feral colonies. We also identified two immune genes that are associated with colony survival and that can potentially be used as biomarkers of health in honey bee colonies.

Materials and Methods

Sample Collection

We established a collaborative community science project to locate feral honey bee colonies across the state of Pennsylvania, United States. Participants reported the locations of feral colonies by submitting GPS locations of colonies into a web portal or via email. To be included in the study, all feral colonies needed to survive at least one winter in wild, unmanaged conditions. We checked each reported colony in early spring to corroborate activity and record overwintering survival (Figure 1). We paired each feral colony with one managed colony located within a seven-mile radius to control for site variation between colonies located in geographic areas with different landscapes and climates. A total of eight pairs of feral and managed colonies (n = 16 colonies) in 2017 and 17 pairs (n = 34) in 2018 were included in the laboratory analyses (Figure 2). In the case where a managed colony was not able to be sampled a second time due to death or other reasons, the feral colony was paired with a different managed colony in the same location (between 2017 and 2018, n = 3; between spring and fall 2018, n = 1). This resulted in 20 unique feral colonies and 24 unique managed colonies being sampled over the course of this study (Supplementary Data Sheet 1). Additionally, we were unable to obtain data on the overwintering status of one managed colony in 2018, therefore it was omitted from the analyses of overwintering survival. Due to the death of either a managed or feral colony in a pair, only two pairs of colonies were sampled in both 2017 and 2018.

Figure 1. Cavities used by feral honey bee colonies included in the study. (A) Colony inhabiting a tree cavity with abundant propolis at the entrance of the nest in Saxonburg (PA, United States). (B) Feral colony nesting inside the wall of a house in New Bethlehem (PA, United States) (C) Entrance of a feral colony nesting in an abandoned shed in Harrison Valley (PA, United States).

Figure 2. Geographic locations of feral and managed colonies sampled in this study. (A) Map of North America indicating the location of Pennsylvania in the northeast of the United States. (B) Location of feral (circles) and managed (squares) colonies sampled in 2017. (C) Location of feral (circles) and managed (squares) colonies sampled in 2018. The color of each shape corresponds to relative DWV levels (2^–ΔΔCT values) recorded in the fall for that colony. Colony locations are approximate and have been adjusted slightly due to overlapping points. GPS coordinates of each colony can be found in Supplementary Data Sheet 1.

We sampled approximately 75 forager bees from the entrance of each colony in the spring (March–June) and fall (August–October) (dates in Supplementary Data Sheet 1). Individuals were sampled with aerial insect nets and transferred to 50 ml conical tubes that were placed on dry ice to preserve RNA quality before long-term storage at −80°C. All sampling sites were on private property and permission was obtained from the land owners. The specific locations and contact information of participants was kept confidential. No protected or endangered species were involved in these studies.

Selection of Pathogens and Immune Genes

To characterize disease dynamics in honey bee colonies, we quantified three pathogens that commonly infect honey bees and negatively impact colony health. Deformed wing virus (DWV) is an RNA virus that is considered the most detrimental honey bee pathogen for its ubiquity, global distribution, and role in overwintering losses of honey bees (Martin et al., 2012; Brutscher et al., 2016). It can be transmitted horizontally and vertically within a colony, but is also efficiently transmitted by Varroa mites which can increase the titer of this virus leading to overt, clinically symptomatic infections (Gisder et al., 2009; Möckel et al., 2011). We also quantified Black queen cell virus (BQCV), another RNA virus that is highly prevalent in adult honey bees globally, but predominantly affects immature bee stages (prepupae and pupae) (Mondet et al., 2014). BQCV is transmitted vertically and horizontally between adults and from adult bees to developing bees, but has not been shown to be transmitted by Varroa mites (Chen et al., 2006). BQCV may also have synergistic interactions with other pathogens due to its correlation with viruses and the fungal parasite Nosema ceranae Fries (D’Alvise et al., 2019). Nosema ceranae is a common microsporidian gut parasite of honey bees, contributing to decreased lifespans in infected bees (Higes et al., 2008; Goblirsch, 2018). We quantified this pathogen in all colonies sampled in 2018.

To characterize immune gene expression in feral and managed colonies, we quantified transcript expression of six genes (argonaute-2, vago, pgrp-s2, pgrp-lc, defensin-1 and hymenoptaecin) from several immune pathways. The genes argonaute-2 (ago2) and vago, from the RNAi pathway, have been shown to be upregulated after viral infection (Brutscher et al., 2015). The gene pgrp-s2 encodes for an upstream recognition receptor involved in activation of the Toll immune pathway, and pgrp-lc encodes a transmembrane protein activator of the Imd (Immune Deficiency) pathway. Both of these genes are upregulated in pathogen-infected honey bees (Evans et al., 2006; Brutscher et al., 2017). Additionally, we quantified genes corresponding to two antimicrobial peptides (AMPs), defensin-1 (Def1) and hymenoptaecin (Hym), produced by the Toll and Imd pathways. These AMPs have key roles in honey bee immune responses to viruses, bacterial and fungal pathogens (Yi et al., 2014; Brutscher et al., 2015).

RNA Extraction and Quantitative PCR

For total RNA extraction, we dissected abdomens from thirty bees per colony by removing the stingers, hindguts, and midguts. Ten abdomens per sample (three samples per colony) were pooled into 2.0 ml tubes with 2.0 mm BashingBead Lysis Tubes (Zymo Research, Irvine, CA, United States) and homogenized using a BeadBlaster^TM24 (Benchmark Scientific, Edison, NJ, United States) at 6.0 m/s for three 30 s intervals. We extracted RNA from homogenate using RNeasy spin columns (QIAGEN, Hilden, Germany), according to the manufacturer’s protocol and eluted into nuclease-free water. We assessed the quantity and quality of RNA using a SpectraMax iD3 Multi-Mode Microplate Reader (Molecular Devices, San Jose, CA, United States).

We quantified the three pathogens and the expression of six immune genes through quantitative reverse-transcriptase PCR (qRT-PCR) using previously developed primer sequences (Table 1). The three RNA extracts from pooled individuals per colony were individually used as templates to produce cDNA using random primers and MultiScribe RT, according to the manufacturer’s protocol (Applied Biosystems, Foster City, CA, United States). A total of 2 μg RNA was used for each cDNA synthesis. A total of 40 ng of cDNA was used for each qPCR reaction. We carried out reactions in 384-well plates using a QuantStudio 5 Real-Time PCR System (Applied Biosystems). Each well contained 5 μl Luna Universal qPCR Master Mix (New England Biolabs, Ipswich, MA, United States), 0.25 μl of each of the forward and reverse primers (10 μM), 2.5 μl nuclease-free H₂O, and 2 μl cDNA template. The following reaction conditions were used: 60 s at 95°C for initial denaturation, then 40 cycles of 15 s at 95°C for denaturation, and 30 s at 60°C for annealing, extension, and data collection followed by a melting curve analysis of 15 s at 95°C, 60 s at 60°C, and 1 s at 95°C to determine the specificity of amplification products. In each plate, we ran all reactions in triplicate and included negative controls of nuclease-free water for each set of primers. After surveying three reference genes (ef1-alpha, eIFS8, and GAPDH1), we determined that elongation factor 1-alpha (ef1-alpha) was a suitable reference gene due to its similar level of expression in all samples, and we used it as the reference gene for these experiments (Table 1).

Table 1. Sequences of forward and reverse primers used for the quantification of pathogens and immune genes through quantitative PCR.

We determined the Ct value for each sample by taking the mean of the three technical replicates. We used the Ct value for the reference gene and subtracted this from the Ct value for the target to generate ΔCt values for each sample. These ΔCt values were then normalized to the managed colony with the lowest relative abundance (highest ΔCt) of the target for that sampling period (spring or fall of the year sampled), generating ΔΔCt values. We then calculated the relative amounts of transcripts using the 2^–ΔΔCT method (Livak and Schmittgen, 2001). As each colony had three samples consisting of ten bees each, we took the average 2^–ΔΔCT value of the three biological replicates and used this for subsequent analyses.

For absolute quantification of DWV copy numbers in samples collected in the fall of each year, we used synthetic DNA corresponding to the sequence of DWV amplified by our primers, in the form of gBlocks gene fragments (Integrated DNA Technologies, Coralville, IA, United States), producing ten-fold dilutions (10¹–10⁷ copies of synthetic DNA) run in each qPCR plate. We calculated the copy number using the formula: copy number = DNA concentration (ng/mL) × 6.02 × 10²³ (copies/mol)/length (140 bp) × 6.6 × 10¹¹ (Wu et al., 2017). The log DNA copy numbers were then plotted with Ct values, producing linear standard curves for each qPCR plate, allowing for the estimation of DWV copy number in each sample.

Statistical Analysis

Due to non-normality, 2^–ΔΔCT values were log-transformed and analyzed through generalized linear models (GLM) with a Gaussian distribution using the glm function in R package ‘stats’ (R Core Team, 2019). The relative expression of each target was analyzed separately using management (feral or managed), sampling time (spring or fall of 2017 or 2018), and their interactions as fixed effects. Two-way ANOVA was used to test for overall effects of management and sampling time on target expression. The estimated marginal means were then calculated with the emmeans function in R package ‘emmeans’ (Lenth, 2020). These values were then used for post hoc tests with a Tukey adjustment for multiple comparisons to determine differences in relative target abundance between feral and managed colonies at each timepoint. To evaluate associations between pathogen levels and immune gene expression among all samples, we calculated the Spearman’s rank correlation coefficients using the cor.test function in R package ‘stats.’

To assess the role of immune gene expression, pathogen levels, and management in overwintering survival, we used a log linear generalized linear mixed model (GLMM) with a binomial distribution using the glmer function in R package ‘lme4’ (Bates et al., 2015). The original full model included data from the fall of both years with overwintering survival as the response variable, and the relative levels of DWV and BQCV, relative expression of all six immune genes, and management as fixed effects, and year of sampling as a random effect. We used a backward model selection to identify the fixed effects that contributed significantly to the model. We calculated the variance inflation factors (VIF) for the model using the vif function in R package ‘car’ (Fox and Weisberg, 2019). All VIF values were less than 3, and thus all fixed effects were kept in the final model. Ultimately, the model with the lowest AIC value was chosen. All analyses were conducted in R version 3.6.2.

Results

Pathogen Levels

Our results indicate that mean DWV levels were significantly higher in feral colonies compared to managed colonies in fall of 2017 (P < 0.05, z = −2.470), but not at other timepoints (P > 0.05) (Table 2 and Figure 3). Levels of BQCV and N. ceranae did not differ between groups (P > 0.05). All 44 colonies tested positive for the presence of DWV and BQCV at all timepoints, while N. ceranae presence was more variable over time (Supplementary Figure 1). Out of the 34 colonies tested, N. ceranae was detected in eight feral and 11 managed colonies in the spring (47.1% and 65%, respectively), and 12 feral and 10 managed colonies in the fall (76% and 59%, respectively).

Table 2. Results of two-way ANOVA on generalized linear models (GLM) assessing the effects of time, management, and their interaction on pathogen levels and immune gene expression.

Figure 3. Boxplots showing relative abundance of pathogens in feral (blue) and managed (red) honey bee colonies for each sampling period. Relative abundance = 2^–ΔΔCT values for each colony (calculations described in main text). In 2017, n = 8 pairs of colonies. In 2018, n = 17 pairs of colonies. Horizontal line of each box represents the median and circles represent outliers. Asterisk denotes statistical significance (*P < 0.05).

Immune Gene Expression

In the spring of 2017, feral colonies exhibited higher average expression of five out of the six immune genes tested (defensin-1, hymenoptaecin, pgrp-lc, pgrp-s2, and ago2), although pathogen levels were not significantly different between managed and feral colonies (Table 2 and Figure 4). In the fall of 2017, expression of defensin-1, hymenoptaecin, pgrp-s2 remained higher in feral colonies, while the gene vago had higher average expression in managed colonies. In 2018, immune gene expression was similar between feral and managed colonies, yet feral colonies had higher average expression of hymenoptaecin and pgrp-s2 in the spring, regardless of the similar levels of pathogens in feral and managed colonies. No significant differences in average gene expression between feral and managed colonies were observed in the fall of 2018. Interestingly, pathogen pressures were also similar at this time point.

Figure 4. Boxplots showing immune gene expression of feral (blue) and managed (red) honey bee colonies for 2017 (n = 8 pairs of colonies) and 2018 (n = 17 pairs of colonies). Relative abundance = 2^–ΔΔCT values for each colony (calculations described in main text). Horizontal line of each box represents the median and circles represent outliers. Asterisks denote statistical significance (*P < 0.05, ***P < 0.001).

Spearman’s correlations (ρ) revealed that DWV levels were positively correlated with the expression of hymenoptaecin in the spring (ρ = 0.32), but not significantly correlated with any other pathogen levels or immune genes, while BQCV levels were positively correlated with N. ceranae levels in the spring (ρ = 0.48). Nosema ceranae levels were also positively correlated with the expression of defensin-1 and pgrp-s2 in the fall (ρ = 0.39; ρ = 0.35, respectively). However, all correlation coefficients between pathogen levels and immune gene expression were low (ρ < 0.4) (Mukaka, 2012), suggesting factors other than pathogen levels likely contribute to immune gene expression (Supplementary Table 1).

Overwintering Survival

Total survival of colonies over the 2017–2018 winter was 63% for both feral and managed colonies. For the 2018–2019 winter, survival was 47% and 38% for feral and managed colonies, respectively. Of the five feral colonies that survived the 2017–2018 winter, two also survived the 2018–2019 winter. Two managed colonies were sampled in both years, and one of these also survived the 2018–2019 winter. Despite the similar overall survival between managed and feral colonies, more feral than managed colonies survived with high copy numbers of DWV (>10⁷) (Supplementary Figure 2).

The best predictive model for overwintering survival included management, DWV levels, and expression of hymenoptaecin, vago, and pgrp-s2 as fixed effects (AIC = 65.5). The level of DWV in a colony was negatively correlated with overwintering survival (DWV: P < 0.05, standard estimation error of coefficient = 0.06). Management and pgrp-s2 expression were also negatively correlated with overwintering survival, although neither were significant predictors (Management: P = 0.105, s.e. = 0.877; pgrp-s2: P = 0.062, s.e. = 0.51) (Figure 5). The expression of hymenoptaecin and vago were both significantly positively correlated with survival (hymenoptaecin: P < 0.05, s.e. = 0.265; vago: P < 0.05, s.e. = 0.757) (Figure 5).

Figure 5. Effect sizes of the generalized linear mixed model (GLMM) showing estimates of each variable along with 95% confidence intervals. Asterisks denote variables that were statistically significant predictors of overwintering survival (*P < 0.05). Estimates of fixed effects are: Management = –1.422; Log (pgrp-s2) = –0.953; Log (hymenoptaecin) = 0.641; Log (vago) = 1.845; Log (DWV) = –0.124.

Discussion

In this study, we investigated the ecoimmunology of feral and managed bees over a 2-year period. Our results show that overall feral honey bee colonies have higher levels of DWV despite yearly and seasonal variation. Alternatively, BQCV and N. ceranae, which were constant across seasons, did not differ with management and were not linked to decreased survival in this study. While we did not investigate the transmission efficiencies of pathogens from feral to managed honey bees, our results provide support for the ability of feral colonies to serve as reservoirs of DWV. We also found evidence of higher immune gene expression in feral colonies, even at timepoints when DWV levels were similar between managed and feral colonies. Further analysis of all colonies revealed that levels of DWV infection were positively correlated with the expression of hymenoptaecin in the spring and N. ceranae levels were correlated with defensin-1 and pgrp-s2 in the fall. The strength of correlations was low suggesting additional factors such as genetic background and environmental conditions play an important role in immune phenotypes. Last, we found significant associations between the expression of two immune genes (hymenoptaecin and vago) and survival in both feral and managed colonies. These genes have been previously identified as differentially expressed in virus-infected honey bees, but this is the first report of expression being correlated with reduced host mortality (Kuster et al., 2014; Ryabov et al., 2014). Feral and managed colonies also had similar probabilities of survival, despite higher DWV titers in feral than managed colonies.

The survival of feral colonies in the presence of high viral titers suggests that feralization may facilitate the expression of traits that confer virus tolerance. Virus tolerance to DWV has previously been observed among mite-resistant honey bees in several locations across many years suggesting that it is heritable (Locke et al., 2014; Russo et al., 2020). Additionally, even without miticidal treatment, several feral colonies had low levels of DWV, suggesting these colonies may exhibit additional traits such as hygienic behavior to reduce the number of Varroa mites, and thus viruses present (Supplementary Figure 2) (Wagoner et al., 2019). Natural selection is a primary driver of evolution in host–pathogen dynamics and allows organisms to potentially respond to many stressors (Papkou et al., 2019). Management practices of domesticated organisms change the fitness landscape and can alter their evolutionary trajectories (Wilkins et al., 2014; Nygren et al., 2015; Milla et al., 2018). Managed honey bees are typically under different environmental conditions compared to their feral counterparts, and while beekeeping practices often aim to limit pathogen abundance, these practices may relax the selective pressures that pathogens exert on managed colonies, potentially delaying host-parasite coevolution (Neumann and Blacquière, 2017). While low input beekeeping management reduces pathogen burdens and oxidative stress compared to commercially managed colonies, highly invasive beekeeping practices may cause added stresses (Taric et al., 2019). Even though the degree of disease management varies greatly among beekeepers (Underwood et al., 2019), Varroa mite management may inhibit adaptations for increased resilience to mites and viruses in managed honey bee colonies (Blacquière et al., 2019). In addition to virus tolerance, other traits such as small colony sizes, frequent swarming, and increased hygienic and grooming behavior appear to be critical for the survival of feral colonies in the presence of high pathogen pressure (Gramacho and Spivak, 2003; Seeley, 2007; Locke, 2016; Loftus et al., 2016; Russo et al., 2020). Future studies should focus on providing information about the molecular mechanisms of disease tolerance in these feral colonies to directly test for the role of feralization in altering selective pressures and leading to different immune phenotypes in honey bees.

Although mechanisms of virus tolerance are still unknown for insects, one possible explanation is the ability of highly infected bees to limit the overreaction of immune responses. Immune effectors, such as pro-phenoloxidase, produced by insects can have cytotoxic effects that work to limit infections, but also cause damage to host tissues (Sadd and Siva-Jothy, 2006; Hillyer, 2016). This self-damage can lead to increased aging and higher mortality as a result of increased inflammation (Alaux et al., 2010; Khan et al., 2017). Virus-tolerant colonies may therefore have mechanisms to limit inflammation-induced damage. A second possible mechanism includes transgenerational immune priming, the development of immune memory via vertical transmission of immunological experiences, and its effects have been demonstrated in several invertebrates (Tetreau et al., 2019). Vertical transmission of DWV and the fact that queens may preside over a colony for several years would favor immune-primed offspring in feral colonies, in contrast to beekeeping operations where queens are replaced frequently. Another potential mechanism includes changes in the virulence of DWV and Varroa mites in feral colonies. Genotypes of DWV are known to differ in virulence and Varroa mite transmission can favor certain strains of the virus (Ryabov et al., 2019). While our study did not assess Varroa mite levels or the genetic diversity of DWV, previous work showed that Varroa mites from managed colonies had increased population growth compared to mites from feral colonies (Dynes et al., 2020), providing evidence for the role of management in selecting for mites with greater reproductive rates. While the specific mechanisms are unclear, the absence of human management and the process of feralization may lead to changes in virus virulence and host tolerance in honey bees.

In addition to different selective pressures experienced by feral and managed colonies, feral colonies can also have different genetic backgrounds, contributing to differences in immune phenotypes and disease outcomes. Previously, López-Uribe et al. (2017) showed that feral and managed honey bees exhibit some genetic differentiation even at short geographic scales. Although we provide evidence for differential immune gene expression between feral and managed colonies, the role of ancestry and genetic diversity in this difference remains unclear. Immune gene expression is heritable in honey bees and virus tolerance may have a heritable basis as well, making it plausible that different genetic backgrounds play an important role in the different immune phenotypes observed between feral and managed colonies (Decanini et al., 2007; Thaduri et al., 2019). The underlying variability and heritability of immunological traits combined with different selection regimes may lead to different evolutionary trajectories for feral and managed honey bees.

We found higher expression of several immune genes in feral colonies compared to managed colonies, suggesting that feralization has led to increased pathogen defenses, although this was not true for all genes and differed with time of sampling. Specifically, we found that the expression of hymenoptaecin was higher in feral colonies, while the expression of vago was higher in managed colonies. However, both genes were associated with increased colony survival. Hymenoptaecin is a general AMP involved in responses to many pathogens including DWV and Varroa mites (Evans et al., 2006; Kuster et al., 2014; Brutscher et al., 2015; Wu et al., 2020). Other studies have shown that this gene is consistently upregulated in response to pathogens and during wounding events in honey bees and may be a potential biomarker to quantify honey bee health (Galbraith et al., 2015; Brutscher et al., 2017; Doublet et al., 2017; Zanni et al., 2017). We also identified the expression of vago as important for the overall survival of colonies. This transcript is expressed upon activation of the RNAi pathway, which recognizes the dsRNA of viruses and leads to the increased expression of vago or its orthologs in mosquitoes, fruit flies, bumble bees, and honey bees (Deddouche et al., 2008; Paradkar et al., 2012; Ryabov et al., 2014; Niu et al., 2016). In DWV-infected honey bees, vago expression has been shown to be significantly increased, providing evidence for its role in antiviral responses (Ryabov et al., 2014). To our knowledge this is the first report of vago expression being directly linked to increased survival in honey bees. These two genes could be considered biomarkers of honey bee health that can be used to predict the ability of a colony to survive the winter (López-Uribe et al., 2020). Additionally, genome-level studies looking for signatures of selection at regulatory regions of immunity in honey bees may also provide important information on mechanisms of pathogen resilience.

Feral organisms offer valuable systems to study potential negative consequences of domestication and anthropogenic influence by examining host–pathogen interactions of organisms that recently escaped managed conditions (Burdon and Thrall, 2008; Gering et al., 2019b). Previous studies of feral honey bees have examined levels of mite infestation, pathogen pressures, or the combination of pathogen pressures and immune gene expression, but the association of host–pathogen dynamics with colony survival was not previously investigated (Seeley, 2007; Thompson et al., 2014; Youngsteadt et al., 2015). Here, we quantified pathogen levels, immune gene expression, and linked this to overwintering survival in managed and feral honey bee colonies. This allowed for the identification of specific genes associated with overwintering survival of honey bees and evidence of virus tolerance in feral colonies, linking immunity, infection, and survival under natural conditions. Further identification of the genetic mechanisms of virus tolerance and biomarkers of bee health can help breeding efforts to focus on increasing these traits in selected honey bee stocks (e.g., Robertson et al., 2020), thus decreasing overall colony losses for the beekeeping industry. Future studies should assess the role of feralization on pathogen dynamics and ecoimmunology in other domesticated species.

Data Availability Statement

The raw data and R script to reproduce the statistical analyses of this study can be found online at the ScholarSphere repository: https://doi.org/10.26207/w2wv-ax07.

Author Contributions

CH and MML-U conceived the study. CH and KCE led data collection. KCE was responsible for correspondence with community members and locating honey bee colonies. CH conducted the data analysis with insight from CR and MML-U. CH and MML-U drafted the manuscript. All authors critically revised the manuscript and gave final approval for publication.

Funding

CH was supported by the Pennsylvania State University Graduate Training Program in Integrative Pollinator Ecology, funded by the Penn State College of Agricultural Sciences’ Strategic Networking Initiative Program. CR was funded through the USDA NIFA Hatch Appropriations under Project PEN04652 and Accession No. 1016243. MML-U was funded through the USDA NIFA Appropriations under Projects PEN04716 Accession No. 1020527 and PEN04620 Accession No. 1011873.

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.

Acknowledgments

We would like to thank John Albright, Ross and Peggy Bell, Ken Bowman, Charles Brudowsky, Carolyn Burgh, Lou Caputzal, Fran and Bob Cooley, Laurie Emerson, Stephen Finke, Bill Fisher, Mardi Frye, Daryln Hoffstot, Robert Hoppe, Rita and Roper Houston, Justin Jacobeen, Bonnie Kaushik, Jackie Kimmel, Thomas Kretchik, Bob Landis, Daryl Martin, Hank and Patty Jo McClenahan, Charles McGee, Josh Nearhoof, Stephen Repasky, Bruce Rodriguez, Amy Sharp, Avi Sommerville, Joshua Tyrone, Lynn Urban, Charlie Vorisek, John Wenzel, and the beekeepers of Pennsylvania (United States) for their participation in this study and the Pennsylvania State Beekeepers Association for their support throughout the project. We are also thankful to Ryan Ford for help with laboratory analyses, Kristen Brochu for help with statistical analyses, Catalina Mejia for help with figure design, Romina Russo for discussions about experimental design, and Brock Harpur, the López-Uribe lab and the reviewers for comments on previous versions of the manuscript.

Supplementary Material

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

Supplementary Figure 1 | Prevalence of pathogens in feral and managed colonies. Deformed wing virus (DWV) and Black queen cell virus (BQCV) were detected in all colonies at all times of sampling. Presence of Nosema ceranae was only tested in colonies sampled in 2018. Nosema ceranae was found in 47.1% and 76.47% of feral colonies in spring and fall, respectively, while 64.7% and 58.8% of managed colonies tested positive for N. ceranae in spring and fall, respectively.

Supplementary Figure 2 | DWV copy number in feral and managed honey bee colonies that died (No) and survived (Yes). The effects of management, overwintering survival, and their interaction were analyzed using GLM, followed by ANOVA to assess their influence on DWV copy number (management: df = 1, F = 2.895, P = 0.096; overwintering survival: df = 1, F-value = 4.281, P = 0.045).

Supplementary Table 1 | Spearman’s correlations of pathogens and immune genes for combined data from feral and managed colonies sampled in 2017 and 2018, separated by season of sampling. Asterisks denote statistical significance (^∗P < 0.1, ^∗∗P < 0.05, ^∗∗∗P < 0.001).

Supplementary Data Sheet 1 | Summary of data from honey bee colonies including: date of sampling, calculated 2^–ΔΔCT gene and pathogen expression values, overwintering survival, abbreviated GPS coordinates, and distance from paired colony.

References