Transcriptomic and Epigenomic Dynamics of Honey Bees in Response to Lethal Viral Infection

Honey bees (Apis mellifera L.) suffer from many brood pathogens, including viruses. Despite considerable research, the molecular responses and dynamics of honey bee pupae to viral pathogens remain poorly understood. Israeli Acute Paralysis Virus (IAPV) is emerging as a model virus since its association with severe colony losses. Using worker pupae, we studied the transcriptomic and methylomic consequences of IAPV infection over three distinct time points after inoculation. Contrasts of gene expression and 5 mC DNA methylation profiles between IAPV-infected and control individuals at these time points – corresponding to the pre-replicative (5 h), replicative (20 h), and terminal (48 h) phase of infection – indicate that profound immune responses and distinct manipulation of host molecular processes accompany the lethal progression of this virus. We identify the temporal dynamics of the transcriptomic response to with more genes differentially expressed in the replicative and terminal phases than in the pre-replicative phase. However, the number of differentially methylated regions decreased dramatically from the pre-replicative to the replicative and terminal phase. Several cellular pathways experienced hyper- and hypo-methylation in the pre-replicative phase and later dramatically increased in gene expression at the terminal phase, including the MAPK, Jak-STAT, Hippo, mTOR, TGF-beta signaling pathways, ubiquitin mediated proteolysis, and spliceosome. These affected biological functions suggest that adaptive host responses to combat the virus are mixed with viral manipulations of the host to increase its own reproduction, all of which are involved in anti-viral immune response, cell growth, and proliferation. Comparative genomic analyses with other studies of viral infections of honey bees and fruit flies indicated that similar immune pathways are shared. Our results further suggest that dynamic DNA methylation responds to viral infections quickly, regulating subsequent gene activities. Our study provides new insights of molecular mechanisms involved in epigenetic that can serve as foundation for the long-term goal to develop anti-viral strategies for honey bees, the most important commercial pollinator.


INTRODUCTION
Infectious disease results from the breakdown of an organism's normal physiological state because of the presence or proliferation of a pathogen. This disruption can be molecularly characterized by transcriptome profiling to permit a systemic understanding of host-pathogen interactions, particularly in combination with other system-level approaches (Arvey et al., 2012;Tanca et al., 2013;Gardy and Loman, 2018). Transcriptomic and epigenomic changes in response to pathogen infections are important to understand because they are interdependent but relate to different temporal dynamics.
DNA methylation of Cytosine (5 mC) is a central epigenetic regulator of gene expression and alternative splicing (Zemach et al., 2010;Li-Byarlay et al., 2013), affecting diverse aspects of organismal function and disease (Edwards and Myers, 2007). Epigenetic programming may also be the target of host manipulations by pathogens (De Monerri and Kim, 2014) and affect host defenses and susceptibility to infection (Merkling et al., 2015;Kuss-Duerkop et al., 2018). The virus induced changes in host immune response and gene expression via DNA methylation is still a new study field. It is understudied whether changing host DNA methylation such as gene expression of certain antiviral immune genes is a main mechanism utilized by DNA or RNA viruses to manipulate immune responses. Genome-wide analyses of methylome and transcriptome indicated that the infection of DNA tumor virus induced hypermethylation of host genes during viral persistence (Kuss-Duerkop et al., 2018). In plants, DNA hypomethylation was reported after 24 h of viral infection (Pavet et al., 2006).
Honey bees (Apis mellifera L.) were the first insects in which a fully functional 5 mC DNA methylation system was discovered (Wang et al., 2006;Kucharski et al., 2008). Social hymenopterans, such as bees, wasps, and ants, have low levels of DNA methylation (Hunt et al., 2010;Glastad et al., 2015;Bewick et al., 2017). The methylated 5 mC sites can be one of three different kinds (CpG dinucleotides, CHG, and CHH, H represents A, T or C) (Lyko et al., 2010). The CpG is the most common kind in honey bees. CpG DNA methylation in insects occurs predominantly in gene bodies, associated with histone modifications and alternative splicing (Lyko et al., 2010;Zemach et al., 2010;Hunt et al., 2013;Li-Byarlay et al., 2013;Glastad et al., 2014;Yan et al., 2014;Glastad et al., 2016). Its phenotypic consequences range from behavioral plasticity (Herb et al., 2012;Opachaloemphan et al., 2018) to alternative development trajectories (Kucharski et al., 2008;Glastad et al., 2019) and disease responses (Galbraith et al., 2015). Especially the study from Herb et al. (2012) showed gene expression differences related to differentially methylated regions (DMRs) and also suggests wider effect because many DMRs are related to transcription factors, in which the alternative splicing may affect their function. Despite this breath of effects, few studies have investigated the role of CpG DNA methylation in insect development, immune system, and disease (Mukherjee et al., 2015;Burggren, 2017;Vilcinskas, 2017;Grozinger and Flenniken, 2019).
Honey bees are the most important managed pollinator in agriculture and crop production (Calderone, 2012). However, recent declines in honey bee health have led to unsustainably high annual colony losses in the apicultural industry (Calderone, 2012;Kulhanek et al., 2017) for which multiple, potentially interacting factors are likely the cause (Cornman et al., 2012;Goulson et al., 2015). One major contributor to widespread colony mortality is pathogenic viral infection. Many of the most important viruses are either associated with or can be directly vectored by parasitic Varroa mites (Bailey et al., 1963;Chen and Siede, 2007;Di Prisco et al., 2011;Gisder and Genersch, 2017). Pathogenic viral stressors affect the morphology, physiology, behavior, and growth of honey bees at different developmental stages (Chen and Siede, 2007;Runckel et al., 2011;Cornman et al., 2012;McMenamin and Genersch, 2015), ultimately leading to reduced colony productivity and mortality. The developmental pupal stage is critical for many viral diseases because the Varroa mite feeds on the pupae for a prolonged time (Di Prisco et al., 2011).
Honey bees employ a combination of individual-and grouplevel defenses against pathogens (Wilson-Rich et al., 2009). The main innate immune pathways of insects are present in honey bees and respond to virus infections (Evans et al., 2006;Brutscher et al., 2015), including the Toll, Janus kinase and Signal Transducer and Activator of Transcription (JAK-STAT), Immune Deficiency (Imd), c-Jun NH2-Terminal Kinase (JNK) pathways, antimicrobial peptides (AMPs), endocytosis, and phagocytosis mechanisms (Brutscher et al., 2015). Activation of some immune pathways, such as the RNAi mechanism, can lead to increased virus resistance (Brutscher et al., 2017). Thus, we understand that some honey bee genes respond to virus infection but a systematic characterization of the temporal dynamics of the immune responses and the general transcriptome and methylome have not been investigated, although such studies can yield important information (Lemaitre and Hoffmann, 2007).
Israeli Acute Paralysis Virus (IAPV) has been associated with colony losses (Cox-Foster et al., 2007;Maori et al., 2009;Tantillo et al., 2015), particularly during the winter (Chen et al., 2014). IAPV belongs to the picornavirus-like family Dicistroviridae, that includes single-strand, positive sense RNA viruses that are pathogenic to a range of insects (Bonning and Miller, 2010). IAPV is relatively common in honey bees and affects all life history stages and castes (Chen et al., 2014). When cells are infected, the positive strand of RNA is released into the cytosol and translated by the host ribosomes. Then the assembled complexes for viral replication synthesize negativesense RNA from which more descendant positive strands of the RNA are produced. These are assembled with viral proteins into virions that are then released to infect other cells (Nagy and Pogany, 2012;Chen et al., 2014). Covert infections persist through vertical and horizontal transmission, but IAPV can also readily cause paralysis and death of infected individuals (Maori et al., 2007). Acute infections of IAPV profoundly affect cellular transcription at the pupal stage , and microarray analyses of larvae and adults indicate significant and stage-specific responses to IAPV infection (Chen et al., 2014). Parallel transcriptomic and epigenomic profiling of adult IAPV infection identified some additional molecular responses, including expression changes in immune genes and epigenetic pathways (Galbraith et al., 2015). However, little overlap of identified responses within and among studies leave the molecular characterization of IAPV infection in honey bees unresolved, and more detailed studies of dynamic responses in the time-course of IAPV infection are needed.
Here, we experimentally test the hypothesis that viral infections can affect not only transcriptomic profiles but also induce DNA methylomic changes in honey bee pupae. Specifically, we compared the complete transcriptome and methylome profiles across three distinct phases of a lethal IAPV infection. Furthermore, we tested if our results indicate significant overlap of immune genes with fruit flies as well as other stages of honey bees under IAPV viral infections.
Relative to the PBS control, IAPV-infected pupae had 198 genes significantly up-regulated and 234 genes significantly down-regulated after 5h; 2,738 genes were significantly upregulated and 3,175 genes significantly down-regulated after 20 h; and 3,021 genes significantly up-regulated and 2,963 genes significantly down-regulated after 48 h ( Supplementary  Figures 1-3). DEG overlap between all significantly up-or down-regulated genes was significant in all pairwise directional comparisons (up-and down-regulated genes separately) among the three time points (Figure 2; up-regulated 5 h vs. 20 h: p < 0.05; all other tests: p < 0.001). For the > 8-fold DEGs, significant overlap was found for all pairwise comparisons among the IAPV up-regulated gene sets (p < 0.001) and also for the down-regulated gene sets at 20 and 48 h (p < 0.001) (Supplementary Table 2).

Gene Enrichment and Functional Analyses
To identify IAPV-induced expression patterns across all genes regardless of significance threshold, GO Mann-Whitney U-tests and a Weighted Gene Co-expression Network Analysis were performed on the expression differences between IAPV infected and control samples at three time points. Post 5h of IAPV infection, multiple biological process GO terms were observed with upregulation of RNA processing (ncRNA, tRNA, and rRNA) as the most significant GO terms (p < 0.001). Only one term was significantly downregulated at this time point: "small GTPase mediated signal transduction." 20 h post IAPV infection, no GO terms passed the strictest significance threshold (p < 0.001) for upregulation, and GO term "translation" was the most significantly downregulated (p < 0.001). Finally, in post 48h IAPV infected samples, there was still no term passing the strictest significance in upregulation, and GO terms of "translation, " "homophilic cell adhesion, " and "aminoglycan metabolic process" were the most significantly downregulated (Figure 3). As to molecular function, significant GO terms associated with downregulation in post 20 and 48 h IAPV infected samples were "structural molecules, " "structural constituent of ribosome, " "structural constituent of cuticle, " and "chitin binding" (Supplementary Figure 4). As to cellular component, significant GO terms associated with downregulation in post 20 and 48 h IAPV infected samples were "ribosomal subunit, " "ribonucleoprotein complex, " and "intracellular non-membrane organelle" (Supplementary Figure 4). Although numerous other terms were less significantly enriched in up-or down-regulated genes, no other GO terms surpassed the most stringent significance criterion.  IAPV reads was observed in the IAPV-infected pupae, with the proportion being similarly high at 20 h (86% ave) and 48 h (90% ave). Additionally, five of nine control samples (all time points included) had a large proportion of reads mapping to DWV.

Comparative Transcriptomic Analyses of Immune Genes With Other Pathogen Infection in Honey Bees and Drosophila
A comparison with a previously published list of 186 immune genes from Evans et al. (2006) showed that 1 and 89 of them were represented in our 5 h, 20 h/48 h DEG lists respectively ( Figure 4A

Temporal Dynamics of Genes and Molecular Pathways Based on Transcriptomic Data
Analyses of the three infection time points revealed dramatic changes of gene expression levels from several key metabolic and immune pathways, including lipid metabolism, the Jak-STAT pathway, and AMPs. With regards to lipid metabolism, the expression of the gene Cyp6as5 was significantly increased in the infected samples from 5 to 20 h, and from 20 to 48 h, similarly to 9 more genes that displayed similar expression profiles (Supplementary Figure 7).
In the Jak-STAT immune pathway, the expression of gene

Changes of DNA Methylation After Viral Infections and Comparative Methylomic Analysis
Total methylated 5 mC sites in CpG, CHG, and CHH settings were 1,955,471; 878,838; and 1,433,905 respectively (Figure 7).
The ratio of methylated CpG sites versus total CpG sites was highest in 5 h PBS control (10%), followed by 5  We identified 523, 5, and 18 differentially methylated regions (DMRs) between IAPV and control groups in the 5, 20, and 48 h post-infection samples, respectively (Supplementary Table 7 for a detailed gene list). Significant overlap of DMRs existed between 5 and 20 h (p < 0.05) and 5 and 48 h (p < 0.0001) comparisons (Supplementary Table 4).
Our comparison of our DMR lists with the previously published immune genes of Evans et al. (2006)

Differential DNA Methylation Leads to Differential Gene Expression Later in Several Immune Pathways
To look into the potential association between gene expression and CpG methylation, we studied whether genes that showed early changes of CpG methylation also exhibited significant changes of gene expression later. A list of 38 immune genes that were hypo-or hyper-methylated at the first time point (5 h post-IAPV infection), were found differentially expressed later (48 h post IAPV infection) (Supplementary Table 8). These cellular pathways included the MAPK (mitogen-activated protein kinase) pathway, JAK/STAT (Janus kinase/signal transducer and activator of transcription) pathway, Hippo signaling pathway, mTOR (mammalian target of rapamycin) signaling pathway, TGF-beta (transforming growth factor-beta) signaling pathway, ubiquitin mediated proteolysis, and spliceosome. Interestingly, we detected 12 genes from this group with differentially spliced patterns with a significant statistical q value ( Table 2).

Motif
To reveal the potential regulation of CpG methylation, we examined the 6-mer motif of each CpG site, including the 2 nucleotides upstream and downstream of each CpG site. Analysis of 6-mer motifs from 5 h DMRs showed no common pattern of the motif (Supplementary Figure 11). Similarly, analysis of motifs from 20 h DMRs containing CpG sites revealed little signal beyond the selected central CpG sites in the middle of 6-mer motifs. However, DMRs from 48 h post-IAPV infection tend to have a T or C following the central CpG (Supplementary  Figure 11).

DISCUSSION
To test our hypothesis that viral infections not only affects gene expression profiles but also DNA methylomic marks in honey bee pupae, we compared the complete transcriptome and methylome responses across three distinct phases (5, 20, and 48 h) of a lethal IAPV infection. While methylation is most affected early during the infection, differential gene expression increases over time and severity of the infection. Many up-and down-regulated genes are involved in RNA processing and translations, proteolysis, and metabolic processes (Figure 3). A collection of genes have differential transcript usage or alternative splicing post IAPV infections (Figure 5). When compared to differentially expressed genes from adult honey bee and fruit fly post viral infection, we have discovered significant overlap of our IAPV-responsive genes with immune genes from these other studies (Evans et al., 2006).
The pupal stage is a critical phase for metamorphosis and an important temporal window for Varroa mite parasitism. Our study, which incorporates and confirms the separate observations of adult or larval honey bees (Cox-Foster et al., 2007;Cornman et al., 2013;Chen et al., 2014;Galbraith et al., 2015), includes new observations that provide a more complete understanding of how the pupae respond to lethal viral attack. Previous studies showed that honey bee pupae can respond to disease with increased immune gene expression to decrease the success of parasites and pathogens (Kuster et al., 2014;Brutscher et al., 2017). Our results reveal that rapid epigenetic changes in response to viral infection precede and presumably trigger profound and widespread transcriptome responses in honey bees. A combination of both transcriptomics and methylomics at multiple time points can reveal the role of and the relation between 5 mC methylation and gene expression profiles under viral infection of insects (Huang et al., 2019).
Alternative splicing in response to viral infection can affect subsequent gene expression through several mechanisms (De Maio et al., 2016). Our results of differential transcript usage or alternative splicing indicated the viral infection can cause large changes in alternative splicing of the bee hosts, which brings a new perspective on molecular mechanisms of IAPV virus pathogenesis. It is known that RNA virus can alter host mRNA splicing (Dubois et al., 2014;Chan et al., 2019). Previous research revealed that DNA methylation is linked to regulate    alternative splicing in honeybees via RNA interference knocking down DNMT3 (Li-Byarlay et al., 2013). The epigenetic changes we report here are interesting because the 5 mC methylome profile indicated a large amount of DMRs after 5 h IAPV infection. However, DMR number dramatically dropped after 20 h and 48 h of infection. This phenomenon may be explained by the cell death and apoptosis after severe infections after 20 and 48 h. We propose a possible explanation is that 5 mC changes as a primary host-defense mechanism responding to viral infection at the onset. During the early stage of viral infection, 5 mC methylation may be the first to react leading to molecular changes genome-wide, potentially interacting with the transcription process such as transcription factor binding. The mechanism on how to induce abnormal methylation at the molecular and cellular levels by viral infection is yet to be determined. A few potential explanations are that pathogenic infection may lead to cellular proliferation and inaccuracy in epigenetic processes, or molecular signaling pathways involved in the infection response affect these epigenetic processes. Previous studies showed that viral infection can change 5 mC methylation in human cancer (Kusano et al., 2006;Toyota and Yamamoto, 2011). Other studies also showed that changes in epigenetics marks can affect immune responses, cell-cycle checkpoints, cell death, and cell fate (Toyota and Yamamoto, 2011;Turner and Diaz-Munoz, 2018).
In honey bees, some evidence links 5 mC methylation marks to gene regulation (Kucharski et al., 2008;Li-Byarlay et al., 2013;Wang and Li-Byarlay, 2015;Li-Byarlay, 2016;Glastad et al., 2019)  but limited reports showed 5 mC methylation are associated with change of gene expressions (Herb et al., 2012). The data here indicates that 38 of the initially differentially methylated genes subsequently exhibit significant changes of gene expression (Table 1), and 12 of these genes displayed significant changes in their splicing events. In general, the epigenetic control of 5 mC methylation and gene expression has been reported for plants and mammals (Allis and Jenuwein, 2016;Saripalli et al., 2020), but the function of 5 mC methylation in the regulation of genes in insects is still not clear (Li-Byarlay et al., 2013;Wang and Li-Byarlay, 2015). One hypothesis is that DNA methylation is involved in cell signaling or involved in the regulation of co-transcription via transcription factors (Toyota et al., 2009;Shukla et al., 2011;Marina et al., 2016). DNA methyl-transferases 1 and 3 are critical enzymes for a functional methylation system. Our data show the IAPV infection reduced the gene expression of these two key enzymes over time. This finding is consistent with previous reports showing that viral infection can change the expression of DNMTs in animal and cell studies (Fang et al., 2012;Tian et al., 2013). The significant fold change of these enzymes is moderate but sufficient to cause dramatic changes in the global DNA methylome, as indicated in previous studies (Wang et al., 2006;Kucharski et al., 2008;Li-Byarlay et al., 2013). The DNA methylation changes in turn can have further consequences for widespread gene expression patterns. We found more DMR post 5 h IAPV infection, when the DNMTs are not yet differently expressed. Possible explanation is that the initial viral infection induced the methylation changes at the initial stage by other potential molecular mechanism such as histone modifications (Paschos and Allday, 2010;Dantas Machado et al., 2015). DNMTs may be regulated to restore homeostasis which may explain why fewer DMRs at the later stages (Matilainen et al., 2017).
At every time-point, we identified numerous differentially expressed genes that confirmed the profound consequences of IAPV infection . Such effects can certainly be caused by cell expansion or loss, rather than direct effects on gene activity. However, the systemic consequences for the health and function of the organism may be similar: the level of immunity might drop due to declining cell numbers or lowered gene activity per cell. Molecular pathways from differentially expressed and methylated genes indicated that the MAPK signaling pathway, Jak/STAT signaling pathway, Hippo signaling pathway, mTOR signaling pathway, TGF-beta signaling pathway, ubiquitin mediated proteolysis, and spliceosome are critical to respond to viral infections. Both MAPK and Jak/STAT signaling pathway are immune pathways with antiviral immune response in honey bees and dipterans (Dostert et al., 2005;Lemaitre and Hoffmann, 2007;Cirimotich et al., 2010;Chen et al., 2014;Galbraith et al., 2015;Brutscher et al., 2015;Grozinger and Flenniken, 2019). The Hippo signaling pathway is highly conserved from insects to mammals, which regulates cell death and differentiation (Munoz-Wolf and Lavelle, 2017). Previous study also showed differentially expressed genes regulated by IAPV in honey bee brood are involved in the mTOR, TGF-beta pathways, and ubiquitin mediated proteolysis (Chen et al., 2014).
In addition, our data showed that IAPV infection can induce Defensin-1, a key antimicrobial peptide responding to viral infections (Kuster et al., 2014;Brutscher et al., 2015). P450 genes are known as responding to oxidative stress in Drosophila during larval development (Li et al., 2008). In honey bees, P450s are involved in xenobiotic detoxification (Johnson et al., 2009;Berenbaum and Johnson, 2015). From a wider perspective, Chikungunya virus infection specifically activates the MAPK signaling pathways in humans (Varghese et al., 2016) and is a stress-responsive part of the intestinal innate immunity of C. elegans (Sakaguchi et al., 2004). Recent report showed that MAPK may play a role in regulating certain P450 genes (Yang et al., 2020). On the other hand, IAPV infection presumably causes major deregulation of the cellular homeostasis through the cooption of the cellular machinery to replicate itself . Thus, deregulation of ribosome biosynthesis (Figure 3) and major energetic pathways is not surprising. Interestingly, transcriptomic analyses of Colony Collapse Disorder, which had been associated with IAPV (Cox-Foster et al., 2007), have shown similar immune responses (Johnson et al., 2009;Aufauvre et al., 2014). The methylomic profile of pupae 48 h after IAPV infection indicated genes involved in anti-parasitoid immune response (Howell et al., 2012;Haddad et al., 2016) to be affected. One of these genes that overlaps with a previous genomic study of honey bee immune genes (Evans et al., 2006) is Peroxin 23, a protein involved in peroxisomal protein import and associated with autophagosome formation in Drosophila muscles and central nervous system (Hazelett et al., 2012;Zirin et al., 2015). Compared to a previous report of the 5 mC methylomic profiles of honey bee adults infected with IAPV (Galbraith et al., 2015), our analysis revealed only one overlapping gene (GB51998), which is related to ATP binding, circadian regulation of gene expression systems (Claridge-Chang et al., 2001), and muscle morphogenesis and function in Drosophila (Schnorrer et al., 2010).
The list of identified DMRs linked to genes associated with ATP binding, phagosome, fatty acid degradation, RNA transport, and RNA degradation. Lipid metabolism, signaling, and biosynthesis are related to virus-host interactions (Chukkapalli et al., 2012). Fatty acid degradation and metabolism also responded to cold stress in other insects and might indicate increased energetic demands of stressed organisms (Vermeulen et al., 2013). Alternatively, IAPV may remodel the host cells by hijacking the lipid signaling pathways and synthesis process for its own replication, similar to virus manipulation of host RNA replication and transport mechanisms (Nagy and Pogany, 2012;Boncristiani et al., 2013).
Our study only represents a first step in the understanding of the intricate interplay of viruses and their honey bee host's physiology. Future research needs to characterize the temporal transcriptome trends in more detail and tissue-specificity. Pupae are particularly relevant but also challenging because the ongoing metamorphosis entails drastic hormonal and presumably transcriptional changes even in the absence of disease. The discovered functional pathways need to be functionally investigated, particularly the putative 5 mC methylation control of immune genes.

Bee Samples and Virus Preparation
All honey bee samples were acquired from unselected experimental hives in the research apiary of the University of North Carolina Greensboro. A previously prepared and characterized IAPV solution in PBS was used for inoculations . White-eyed pupae were either injected with PBS solution as control treatment or with 10 4 genome copies of IAPV in 1.0 µl of the viral suspension per pupa. Pupae were maintained on folder filter papers in a sterile lab incubator at 32 • C and 60% R.H. After 5, 20, and 48 h of injections, each individual pupal sample was flash-frozen in liquid nitrogen and then stored in a −80 • C freezer. Figure 1 illustrates the experimental design.

RNA-Seq Library Preparation and Sequencing
Treatment control (PBS-injection) and IAPV-infected pupae were compared across three time points (5, 20, and 48 h). Three pupae were collected for each treatment and each time point for a total of 18 samples. Dual extraction of nucleic acid (RNA and DNA) from each pupa was carried out using the MasterPure dual extraction kit (Epicenter, MC85200). Briefly, whole pupae were individually homogenized and lysed first by using a micropestle. The lysed cells were mixed with extraction buffer and centrifuged to remove proteins and undesired debris of the cell. Subsequent RNA and DNA fractions were treated with DNase and RNase accordingly. RNA-seq library preparation and sequencing were performed by the Genomic Sciences Laboratory of the North Carolina State University. Total RNA from each pupa was used to generate tagged Illumina libraries, using the NEBNext R RNA Library Prep Kit (New England Biolab). All RNA-seq libraries were sequenced in two flow cells on the Illumina Next-Seq 500 (paired-end 150 bp length reads).

BS-Seq Library Preparation and Sequencing
Bisulfite (BS) conversion of the genomic DNA from the abovementioned dual extraction of each sample was carried out by using the EZ DNA Methylation-Lightning Kit (Zymo Research, D5030, Irvine, CA, United States). BS-converted DNA (100 ng) of each sample was used for preparing the BS-seq library via the TruSeq DNA Methylation Kit (Illumina, EGMK91324, San Diego, CA, United States). The quality of the libraries was checked using Qubit assays (Q32850, Life Technologies) and with the 2100 Bioanalyzer (Agilent Technologies, Inc.). The pooled 18 libraries were run on two replicate lanes of the Illumina NextSeq 500 (150 bp paired-end reads) by the Genomic Sciences Laboratory at NCSU.

Bioinformatic Analyses
Overall quality control analysis of all resulting data was carried out by using FastQC. 1 Raw data were trimmed by Trimmomatic (java -jar PE -phred33 ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 MINLEN:36). TopHat2 v.2.0.12 (Kim et al., 2013) was used to align trimmed RNA-seq reads to the Apis mellifera 4.5 reference genome. General alignment parameters were chosen to allow three mismatches per segment, 12 mismatches per read, and gaps of up to 3 bp (tophat -read-mismatches 12segment-mismatches 3 -read-gap-length 3 -read-edit-dist 15 -b2-sensitive) to minimize the possibility that small sequence differences of our samples with the reference genome would bias expression estimates. These settings resulted in the successful alignment of 35% of raw reads (117-143 million per sample). Reconstructed transcripts were aligned with Cufflinks (v2.2.1) to the Apis mellifera OGS 3.2 allowing gaps of up to 1 bp (cufflinks -overlap-radius 1 -library-type fr-firststrand). To annotate promoters and transcription starts, we employed the "cuffcompare" command using OGS3.2 gff file (cuffcompare -s -CG -r).
Cuffdiff2 (v.2.1.1) (Kim et al., 2013) was used to test for differential expression between IAPV and control groups (n = 3 pupae treated as biological replicates at each of the 3 different time points) with the multiple mapping correction and a false discovery rate threshold set at 0.05. All three IAPV groups (IAPV_5h, IAPV_20h, and IAPV_48h) were compared to their corresponding sham-control groups (PBS_5h, PBS_20h, PBS_48h) to identify differentially expressed genes (DEGs). The data were explored and visualized with CummeRbundv.2.0.0 2 and custom R scripts (Rsoftwarev.2.15.0). 3 Differentially expressed genes with similar expression temporal dynamics were analyzed by cufflinks in R. Total of output of gene numbers is set to be 10 (for example, consider gene GB49890 (cyp6as5), we can find nine other genes in the database with similar expression patterns, then plot the expression patterns in Supplementary Figure 7).
The reference genome and official gene set were the same as described before. The sequencing raw data are published at the NCBI database (Sequence Read Archive (SRA) submission: SUB3404557, BioProject: PRJNA429508). All codes for the bioinformatic analyses are online. 5 To determine overlap between DEG lists, either up-or down-regulated genes were compared between different time points. This directional overlap analysis avoids artifacts that can confound undirected DEG comparisons (Lawhorn et al., 2018). The analysis was compared focusing on DEGs with a more than 8-fold difference, and with DMRs (either hyperor hypo-methylated). All overlap analyses were performed with hypergeometric tests to calculate the p-value for overlap based on the cumulative distribution function (CDF) of the hypergeometric distribution as published previously (Li-Byarlay et al., 2013). GO Mann-Whitney U tests were carried out on the gene expression data using codes published in Github. 6 These analyses search for patterns of up-or down-regulation in genes associated with particular GO-terms based on their signed, uncorrected p-values regardless of significance thresholds (Wright et al., 2015).
For differentially transcript usage or spliced genes, transcriptlevel counts were quantified using Salmon1 v1.2.1 (Patro et al., 2017) against the A. mellifera HAv3.1 annotation. Viral RNAs were included in the transcriptome used for alignment but were omitted from isoform analysis. Salmon results were imported directly into the IsoformSwitchAnalyzeR R package (Vitting-Seerup and Sandelin, 2017; Vitting-Seerup and Sandelin, 2019), which normalizes transcript counts based on abundances and uses the TMM method to adjust effective library size. The version of IsoformSwitchAnalyzeR used was 1.11.3, obtained from GitHub commit 5a7ab4a due to a bug in version 1.11.2 on Bioconductor. Although CDS were annotated in the GFF file for A. mellifera, these caused problems in IsoformSwitchAnalyzeR due to inclusion of stop codons, as well as some transcripts having no annotated CDS or CDS that were obviously incorrect. Therefore, the analyzeORF function was used to predict open reading frames de novo from the transcript sequences, which in the vast majority of cases resulted in agreement with the published annotation. Nucleotide and amino acid sequences for all transcripts for significant genes were then exported for analysis with external tools: CPAT9 was used to determine whether a transcript was likely to be coding (Wang et al., 2013). It was run from the webserver using the Drosophila (dm3, BDGP release 5) assembly. A coding probability of 0.7 was used as the cutoff for classifying a transcript as coding or non-coding. Pfam10 (Punta et al., 2012) version 32.0 was used to predict protein domains, using pfam_scan v1.6 run on the Biocluster.

Statistical Analysis of Differential Transcript Usage
A design matrix was constructed in which the condition was the combination of treatment and time point (six conditions total), 5 https://www.researchgate.net/publication/343237015_RNAseq_and_BSseq_ codes_for_public 6 https://github.com/z0on/GO_MWU and DWV infection (determined by alignment of reads to the DWV mRNA) was included as a covariate. Sample degradation was excluded as a covariate as it introduced singularity into the model at 48 h. The contrasts tested were IAPV_5h -PBS_5h, IAPV_20h -PBS_20h, and IAPV_48h -PBS_48h. The preFilter function in IsoformSwitchAnalyzeR was used with default settings to filter transcripts and genes. Genes were removed if they only had one isoform or if they did not have at least 1 RPKM in at least one condition. Transcripts were removed if they had no expression in any condition. After filtering, 15,808 isoforms belonging to 4,667 genes remained.
Differential transcript usage (DTU) was assessed using DEXSeq from within IsoformSwitchAnalyzeR (Anders et al., 2012;Ritchie et al., 2015;Vitting-Seerup and Sandelin, 2019). Genes were only retained for further analysis if they had an FDR < 0.05 for DTU, and if isoform usage (proportion of counts from a gene belonging to a given isoform) changed by at least 0.1 for at least two isoforms in opposite directions (i.e., the isoforms "switched").

Enrichment of Overlapping Genes Among Different Studies
To test whether our DEGs and DMRs had significant overlap with known immune genes, these lists were compared to published immune genes (Evans et al., 2006) using hypergeometric tests. Additionally, a directional overlap analysis was performed on our DEGs (separated by up-or down-regulation) with a previous RNA-seq study of IAPV infection in adult honey bee workers (Galbraith et al., 2015).
To identify potential motifs that are targeted for epigenetic modification after IAPV infection, differentially methylated regions (DMR) were analyzed post 5, 20, and 48 h IAPV infection generating 6-mer sequence logos with the online tool weblogo. 7

DATA AVAILABILITY STATEMENT
The sequencing raw data are published at the NCBI database