Transcriptome-level assessment of the impact of deformed wing virus on honey bee larvae

Deformed wing virus (DWV) prevalence is high in honey bee (Apis mellifera) populations. The virus infects honey bees through vertical and horizontal transmission, leading to behavioural changes, wing deformity, and early mortality. To better understand the impacts of viral infection in the larval stage of honey bees, artificially reared honey bee larvae were infected with DWV (1.55 × 1010 copies/per larva). No significant mortality occurred in infected honey bee larvae, while the survival rates decreased significantly at the pupal stage. Examination of DWV replication revealed that viral replication began at 2 days post inoculation (d.p.i.), increased dramatically to 4 d.p.i., and then continuously increased in the pupal stage. To better understand the impact of DWV on the larval stage, DWV-infected and control groups were subjected to transcriptomic analysis at 4 d.p.i. Two hundred fifty-five differentially expressed genes (DEGs) (fold change ≥ 2 or ≤ -2) were identified. Of these DEGs, 168 genes were downregulated, and 87 genes were upregulated. Gene Ontology (GO) analysis showed that 141 DEGs (55.3%) were categorized into molecular functions, cellular components and biological processes. One hundred eleven genes (38 upregulated and 73 downregulated) were annotated by KO (KEGG Orthology) pathway mapping and involved metabolic pathways, biosynthesis of secondary metabolites and glycine, serine and threonine metabolism pathways. Validation of DEGs was performed, and the related gene expression levels showed a similar tendency to the DEG predictions at 4 d.p.i.; cell wall integrity and stress response component 1 (wsc1), cuticular protein and myo-inositol 2-dehydrogenase (iolG) were significantly upregulated, and small conductance calcium-activated potassium channel protein (SK) was significantly downregulated at 4 d.p.i. Related gene expression levels at different d.p.i. revealed that these DEGs were significantly regulated from the larval stage to the pupal stage, indicating the potential impacts of gene expression levels from the larval to the pupal stages. Taken together, DWV infection in the honey bee larval stage potentially influences the gene expression levels from larvae to pupae and reduces the survival rate of the pupal stage. This information emphasizes the consequences of DWV prevalence in honey bee larvae for apiculture.


Materials and methods
Artificial rearing of A. mellifera larvae. Honey bee larvae were collected from healthy honey bee (A. mellifera) colonies of NIU apiaries (National Ilan University, Taiwan). Healthy honey bee colonies were defined as containing nine frames of the hive with ca. 25,000 workers and a normal spawning queen. More than 300 1-day-old larvae were collected into 24-well culture plates (10 larvae/per well) containing a basic larval diet (BLD) 35 . The method for collecting 1-day-old larvae followed the protocol proposed by Ko et al. 36 . The artificial rearing method was modified based on methods proposed by Hanley et al. and Ko et al. 36,37 . The 1-day-old larvae were reared at 34-35 °C and 95% relative humidity in freshly prepared BLD. The BLD was changed daily for 3 days. The larvae were then transferred to new 24-well culture plates (5 larvae/per well) containing freshly prepared BLD for the DWV infection experiment.
Screening of DWV-carrying honey bee samples. All honey bee colonies used in this study were screened for DWV and six other common honey bee virus infections-ABPV, BQCV, CBPV, IAPV, SBV and Varroa destructor virus-1 (VDV-1)-by reverse transcription polymerase chain reaction (RT-PCR) as described below. Midgut tissue was collected from three honey bees in 500 µL of 0.1% phosphate buffered saline (PBS) solution as one sample and homogenized with a homogenizing pestle. Total RNA was extracted from 200 µL of homogenized midgut tissue using an RNA Isolation Kit (FairBiotech, TW), and the remaining part was stored at − 80℃. The extracted RNA (1 μg) was then treated with DNase I at 25 °C for 15 min and inactivated by EDTA at 65 °C for 10 min. DNase I-treated total RNA was reverse transcribed using a GScript RTase kit (GeneDi-reX, USA) following the manufacturer's instructions. The reaction was incubated at 42 °C for 1.5 h. and then terminated at 70 °C for 15 min. RT-PCR was performed using eight honey bee viral gene-specific primer sets (Supplementary Table 1). PCR amplification was performed as follows using a Primus 96 Plus Thermal Cycler (MWG-Biotech, DE): an initial preheating step at 94 °C for 5 min, and then 30 cycles at 95 °C for 30 s, 50 °C for 30 s, and 72 °C for 1 min, followed by a 5-min final extension at 72 °C and storage at 20 °C. The PCR product was analysed by electrophoresis on a 2% agarose gel in 1 × TAE buffer. Ten DWV-positive samples without other viral signals were collected and subjected to quantification of viral RNA copies.

Quantification of DWV and viral infection based on a plasmid standard curve. The selected
DWV-positive samples were subjected to real-time quantitative polymerase chain reaction (RT-qPCR) to calculate the viral copy number. For RT-qPCR, partial DWV gene fragments were amplified by PCR and then cloned using a T&A clone kit (RBC Bioscience, TW); the ligated plasmid DNAs were transformed into Escherichia coli www.nature.com/scientificreports/ DH5α (RBC Bioscience, TW). Bacterial colonies of the expected recombinant plasmids were chosen for further analysis by commercial sequencing. The sequence-confirmed plasmids were then used to generate a qPCR standard curve ( Supplementary Fig. 1). Undetermined DWV samples were diluted tenfold and tested using Applied Biosystems StepOne Real-time PCR systems with DyNAmo ColorFlash SYBR Green qPCR reagents (Thermo Scientific, USA). The qPCR programme was set as follows: 95 °C for 1 min and 40 cycles of 95 °C for 15 s and 60 °C for 1 min. The results for serially diluted standard plasmids were used to construct a standard curve. The linear standard equation for DWV quantification was generated by plotting the crossing point (Cp) versus the log 10 of the initial plasmid copy number as follows: y = -3.3831x + 45. 306, R 2 = 0.9992. The range of the assay was from 10 4 to 10 10 copies per sample. Therefore, the viral genome copy numbers of undetermined DWV samples were calculated based on the formula for the plasmid standard curve. The quantitated samples were used to infect honey bee larvae.
DWV infection of honey bee larvae. The larvae of the honey bee colonies were detected for the prevalence of 7 viruses before DWV infection as described above (Supplementary Table 1; Supplementary Fig. 2). Honey bee colonies without other identified viral signals were used for viral infection. One-day-old larvae were collected and divided into control (noninfected) and infected groups (approximately 150 larvae for each) and then reared as mentioned above. Each honey bee larva was inoculated in 60 µl of 0.1% PBS solution containing 1.55 × 10 10 DWV copies, while the control group was inoculated in 60 µl of 0.1% PBS solution only. After the larvae consumed the PBS solution for 12 h, fresh BLD was added to the larva. All experiments were performed in triplicate. The survival rates of the infected and control groups were observed and recorded daily until pupation. The survival analysis was performed using SPSS (IBM SPSS Statistics version 20).

Detection of viral copies after DWV infection.
To determine the viral copy number, three DWVinfected larvae were collected as one sample and in triplicate at 0 (the larvae infected with DWV at 3 h. post infection), 2, 4 and 9 d.p.i. (3-, 5-, 7-and 12-day-old larvae). After RNA extraction, conventional RT-PCR and quantitative analysis of DWV were performed as described above using the DWV primer set (Supplementary Table 1).
Preparation of RNA samples for next-generation sequencing. The honey bee larvae were collected at 4 d.p.i. (7-day-old larvae) and homogenized using Tissue Lyser II (QIAGEN, DE) at f = 30/s for 1 min, and this procedure was repeated three times. Total RNA was extracted using TRIzol reagent (Invitrogen Life Technologies, USA) following the manufacturer's instructions. The extracted RNA sample was then purified using the RNeasy Mini Kit (QIAGEN, DE). The RNA purity and integrity were checked using a Nanodrop 1000 Spectrophotometer V3.5 (Thermo Scientific, USA) and Qubit 2.0 Fluorometer (Invitrogen, USA). The total RNA samples from three DWV-infected larvae were pooled as one replicate, and two biological replicates were used for library construction. The size distribution of the RNA samples was determined using the Fragment Analyzer system (Agilent, USA). The RNA samples were then used to construct a library for transcriptome sequencing.
Next-generation sequencing and raw data processing. The mRNA from the DWV-infected and noninfected larvae was further purified for sequencing library construction using an TruSeq stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's protocol. Finally, the DNA fragments were selectively enriched by PCR and then qualified and quantified using the Fragment Analyzer system (Agilent, USA) and qPCR, respectively, for DNA library quality control analysis. These two samples were sequenced in parallel using an Illumina NovaSeq 6000 sequencer to generate high-throughput transcriptome sequencing data with a read length of 151 bp by paired-end (PE) technology. The sequencing adaptors were trimmed from raw PE reads using Trimmomatic 38 , and then read quality checking from trimmed PE reads was performed using PRINseq 39 . The PE reads were also mapped to the genomes of 7 viruses to screen for viral infections.
Differentially expressed gene (DEG) analysis. The quality PE reads were mapped to the honey bee genome index (A. mellifera Amel_HAv3.1) using HISAT2 40 with corresponding gene annotation information (GCF_003254395.2_Amel_HAv3.1_genomic.gff). The read counts of genes were computed using HTSeq 41 . To exclude genes with extremely low expression in the analysis, the most expressed replicate in either the control or DWV-infected sample must have had a minimum of 10 reads (after normalizing for sequencing depth). The read counts were then subjected to differentially expressed gene (DEG) analysis using edgeR 42 ( Supplementary  Fig. 3). Multiple testing errors were corrected using the false discovery rate (FDR). Statistical analysis of fold change was performed, and a significant difference in the expression of these two groups was defined as an FDR < 0.05. Additionally, we considered differentially expressed genes as having a log2 ratio ≥ 1 or ≤ -1 (fold change ≥ 2 or ≤ -2). To facilitate logarithmic transformation, a value of 1 was added to the CPM value of the remaining genes. The gene expression heat map was then generated using the R package 43 . DEGs were further subjected to GO enrichment analysis (http:// geneo ntolo gy. org/ page/ go-enric hment-analy sis) and KEGG Automatic Annotation Server (KAAS, http:// www. genome. jp/ tools/ kaas/) analysis for gene orthologue assignment and pathway mapping 44-46 . qRT-PCR validation. To Table 2). The RNA of the DWV-infected and control groups at 0 (larvae infected with DWV at 3 h. post infection), 2, 4, 6 and 9 d.p.i. was extracted as described above (three larvae were pooled as one sample). Before performing RT-qPCR, the RNA was treated with DNase I (Invitrogen Life Technologies, USA) following the manufacturer's instructions to reduce genomic DNA contamination. The DNase I-treated total RNA samples were reverse transcribed using the GScript RTase Kit (GeneDireX, USA) following the manufacturer's instructions. The reaction mixture was incubated at 42 °C for 1.5 h., and then the reaction was terminated at 70 °C for 15 min. Real-time qPCR was performed using a Thermo Scientific Verso SYBR Green 1-step qRT-PCR ROX Mix Kit (Thermo Scientific, USA) in a 96-well Bio-Rad CFX96 Real-Time PCR System (Bio-Rad, USA). All reactions were performed in five replicates. The relative gene expression levels were calculated using the 2 −△△Ct method 47 .

Results and discussion
DWV infection in honey bee larvae. Honey bee larvae were infected with 1.55 × 10 10 copies/per larva at 3 days old (0 d.p.i.). The number of surviving larvae from both the infected and noninfected groups was recorded daily from 0 to 22 days after infection and analysed using Kaplan-Meier survival analysis. A similar survival rate for the DWV-infected and noninfected groups was observed, but a significantly lower survival rate (chi-square = 9.785; 1 d.f.; P = 0.002) was identified in the DWV-infected group during the period between the late larval stage and prepupal stage (Fig. 1). The replication of DWV was detected in the infected group at the larval stage (0, 2 and 4 d.p.i.) and pupal stage (9 d.p.i.) (Fig. 2). The DWV copies gradually increased from 2 d.p.i. to 4 d.p.i. and dramatically increased to the pupal stage (9 d.p.i.). Based on the results, the DWV copies changed dramatically from the last phase of the larval stage (4 d.p.i.) to the pupal stage (9 d.p.i.), and the mortality was increased after the pupal stage (Figs. 1 and 2); therefore, the transcriptomic experiment focused on the last phase of the larval stage to investigate the gene expression profile in this period.
The survival rates for both groups were similar until 7 days of age (4 d.p.i.); however, in the 10-day-old samples (pupal stage), the survival rate was greatly decreased in the infected group, while DWV copies increased significantly during the pupal stage. Honey bee eggs were exposed to DWV predominantly by virus transovum 25 ; thus, honey bee eggs show a high risk of DWV infection once larvae hatch. Ten days after DWV infection, we observed that the survival rate was substantially decreased in the infected group, suggesting that DWV may have negative effects starting in the early life stage of honey bees, revealing the risk of infection to the honey bee population (Fig. 1).
Statistical analysis of next-generation sequencing. In total, 47,569,753 and 44,405,693 raw reads were generated from the DWV-infected and noninfected groups, respectively. Quality paired-end reads were obtained by adaptor trimming and quality checking. The quality reads were then mapped to the honey bee gene index and corresponding annotation database (GCF_000002195.4_Amel_4.5_genomic.gff). The average mapping rates to map RNA transcripts were 82.5% for the DWV-infected group and 91% for the noninfected group (Supplementary Table 3). The mapped reads were located in 9820 genes of honey bee mRNA and then were processed for further DEG analysis. The PE reads were further mapped to the genomes of common honey bee viruses and confirmed that the control and DWV infection groups did not have other honey bee viral infections (Supplementary Table 4). www.nature.com/scientificreports/

Detection of differentially expressed genes (DEGs).
To identify the significantly differentially expressed genes, expression profile comparisons between the infected and noninfected groups were performed. The DEG analysis results showed that most of the gene expression levels were similar between these groups; however, some of the genes showed different expression levels ( Fig. 3; Supplementary Fig. 6). Two hundred sixtyfive genes (0.26%) exhibited significant fold changes in the expression profile (fold change ≥ 2 or ≤ -2) between the DWV-infected and noninfected groups (  Tables 5 and 6). The DEGs with higher fold changes (fc ≥ 4 or ≤ -4) were summarized; of these DEGs, 22.99% (20 genes) of genes were more highly upregulated, and 17.26% (29 genes) were more highly downregulated (Tables 1 and 2). Thus, the transcriptomes of DWV-infected and noninfected honey bee larvae are moderate, but some physical mechanisms may be subjected to the effect of DWV infection.   www.nature.com/scientificreports/ upregulated with fold enrichments of 33.2, 25.9 and 25.9, respectively (FDR < 0.05; Table 3 Table 4). The GO enrichment results indicated that DWV infection may influence the cellular component and membrane structure and affect the binding affinity of immaterial or material entities with granularity of cells. Therefore, impacting the transmission of the signal from one side of the cell membrane to the other side initiates a change in cell activity. Rabies virus (RABV) infectivity is drastically decreased after metabotropic glutamate receptor subtype 2 (mGluR2) siRNA knockdown in cells, and mGluR2 modulates a functional cellular entry receptor for RABV 48 . GO analysis revealed a negative effect of the plasma membrane and cell surface, which may suppress DWV infection. Additionally, the catalytic activity and response to stimulus-related pathways were enriched, suggesting that honey bee larvae attempted to respond to DWV infection. Based on the survival data, honey bee larvae successfully escaped DWV infection; however, the protective mechanism may fail to work and lead to a lower survival rate during metamorphosis (Fig. 1).  Tables 9 and 10). Among these mapped DEGs, pathways of many metabolism-related enzymes possessed the highest numbers of hits, followed by "biosynthesis of secondary metabolites", "glycine, serine and threonine metabolism", "Rap1 signaling pathway", "cAMP signaling pathway", "calcium signaling pathway" and "inositol phosphate metabolism". Further analysis of the DEGs involved in these pathways revealed two upregulated DEGs [myo-inositol 2-dehydrogenase (iolG, LOC552024) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH, LOC413924)] were involved in 5 and 12 pathways, respectively (Table 1; Supplementary Table 9). From our data, IolG was www.nature.com/scientificreports/ mapped to 5 pathways, including microbial metabolism in diverse environments (ko01120) and inositol phosphate metabolism (ko00562). The iolG gene is involved in inositol catabolism in Aerobacter aerogenes 49 . Inositol monophosphatase mediates the synthesis of myo-inositol, which is likely involved in lipid metabolism 50 . A recent study also demonstrated that inositol phosphates enhances RNA virus (HIV-1) assembly and prevents viruses from cellular defence 51 , suggesting that upregulation of IolG after honey bee infection with DWV may facilitate viral assembly and spread. According to the analysis, GAPDH is involved in 12 different pathways, such as microbial metabolism in diverse environments (ko01120), carbon metabolism (ko01200), and glycolysis/gluconeogenesis (ko00010), indicating the importance of this gene in the biological function of the cells 52 .

KEGG analysis of differentially expressed genes.
GAPDH is a suitable housekeeping gene in honey bees; however, GAPDH is associated with and regulated after viral infection 53,54 . Therefore, hijacking host cellular resources by upregulating IolG and GAPDH after honey bee infection with DWV may facilitate viral infection.
Regarding the downregulated DEGs, multiple inositol polyphosphate phosphatase 1 (MINPP1, LOC409751) involved in 3 pathways (Table 2; Supplementary Table 10) was more relevant to viral infection. MINPP1 is an enzyme that hydrolyses abundant metabolites, such as inositol pentakisphosphate and inositol hexakisphosphate 55 . Inositol hexakisphosphate stimulates both immature and mature HIV-1 particle assembly to become an infectious form 51 . Suppressing MINPP1 expression may be a mechanism of honey bees against DWV infection.
qRT-PCR validation of DEGs. In total, 20 genes were upregulated and 29 genes were downregulated in DWV-infected honey bee larvae. To further validate the transcriptome data, four upregulated [wsc1 (LOC100577331), cuticular protein (LOC724464), iolG (LOC552024) and GAPDH (LOC413924)] and three   5 and 6). Regarding the upregulated DEGs, four genes increased the expression levels from 0 to 4 d.p.i. and reached a high peak at 4 d.p.i., while a declining pattern was observed at the pupal stage (Fig. 5). Most of the downregulated genes were consistent with expectations; the related gene expression levels were suppressed from 0 to 4 d.p.i., while CYP6A1 altered the related expression level from upregulation (2 d.p.i.) to downregulation (after 4 d.p.i.) (Fig. 6).
Proteins containing WSC domains have many biological functions, such as mediating intracellular responses to environmental stress 56,57 . Interestingly, the functions of wsc-related genes play a role in the heat shock response involved in the Pkc1-MPK1 signalling pathway 57 . In honey bees, the MAPK cascade transmits signals from the outer cell surface to the nucleus is an antiviral mechanism related to endocytosis 58 . The wsc1 gene was upregulated from 2 to 6 d.p.i. and then was suppressed at 9 d.p.i. (fc = − 7.61), indicating that DWV infection triggers the host defence system; however, downregulation of the wsc1 gene reflects a significant drop in the survival rate of DWV-infected honey bees in the pupal stage.
Cuticular protein is the protein plays an important role in the protection of insects. During DWV infection, the cuticular protein gene of honey bees was upregulated at 2 d.p.i. to 4 d.p.i. (fc = 2.37 and 12.63, respectively) but downregulated at 6 d.p.i. to 9 d.p.i. (fc = − 3.46 and − 5.42, respectively) (Fig. 5). Additionally, cuticular protein is downregulated when honey bees are infected with the microsporidian Nosema ceranae 59 , and the data are similar to the related expression level of cuticular protein in honey bee pupae infected with DWV, suggesting a similar outcome under this gene suppression. As mentioned above, inositol phosphates enhances RNA virus (HIV-1) assembly and prevents viruses from cellular defence 51 . Additionally, GAPDH is associated with and regulated after viral infection 53,54 . Furthermore, relative quantitative RT-PCR showed that the iolG gene was upregulated at 2 d.p.i., 4 d.p.i. and 9 d.p.i. (fc = 2.15, 9.29 and 3.26, respectively), and GAPDH was upregulated at 4 to 6 d.p.i (fc = 5.62 and 2.15, respectively). Therefore, upregulation of IolG and GAPDH after honey bee infection with DWV may facilitate the viral infection process.
Regarding the three downregulated genes, the related expression level showed a continuous declining pattern from the beginning (larval stage) of DWV infection to the pupal stage, except that CYP6A1 was upregulated at 2 d.p.i. (fc = 2.89) (Fig. 6). CYP6A1 is a member of the cytochrome P450 family and performs various enzyme reactions. Previous studies have indicated that Drosophila melanogaster P450 enzymes are related to developmental www.nature.com/scientificreports/ processes and are involved in the detoxification of foreign compounds 60 . In honey bees, CYP6A1 upregulation may be induced by DWV infection; however, downregulation is mostly observed during DWV infection, suggesting that insufficient detoxifying proteins, such as CYP6A1, may reduce the survival fitness of honey bees. The SK channel plays a fundamental role in all excitable cells, is potassium selective and is activated by increased intracellular calcium levels, such as during action potentials 61 . Additionally, dendritic cells (DCs), which are controlled by Ca 2+ signalling, play an important role in innate and adaptive immunity. Therefore, alterations of cytosolic Ca 2+ can trigger immune suppression or switch off DC activity 62 . From our data, SK was  www.nature.com/scientificreports/ downregulated at all time points, and significant downregulation was observed at 4 d.p.i. and 9 d.p.i. (fc = − 3.51 and − 10.64, respectively), indicating that the signals for the immune system may be dysregulate after DWV infection.
As mentioned above, MINPP1 hydrolyses inositol pentakisphosphate and inositol hexakisphosphate 55 . Inositol hexakisphosphate stimulates both immature and mature HIV-1 particle assembly to become an infectious form 51 . Based on the results of relative quantitative RT-PCR, MINPP1 expression was downregulated at all time points. Suppressed MINPP1 expression may be a mechanism of honey bees against DWV infection. Therefore, the gene expression pattern shows that continuous downregulation of MINPP1 in the pupal stage may affect the fitness of honey bees during pathogen infections.

Conclusions
In this study, the effect of DWV on honey bee larvae was evaluated based on transcriptomic level assessment. DEG identification showed that the transcriptomes of DWV-infected and noninfected honey bee larvae were moderate, while some physical mechanisms may be subjected to the effect of DWV infection. By validating these DEGs, wsc1, cuticular protein and iolG were confirmed to be significantly upregulated, and SK was significantly downregulated at 4 d.p.i. and continuously to the pupal stage. Additionally, these DEGs were significantly regulated in the pupal stage, indicating the potential effects of the gene expression levels from the larval to the pupal stages. In conclusion, DWV infection in the honey bee larval stage may influence the gene expression levels from larvae to pupae and reduce the survival rate of the pupal stage. This information highlights the consequences of DWV prevalence in honey bee larvae for apiculture. Therefore, detection and management strategies for DWVinfected honey bee colonies at an early stage must be developed to avoid colony losses.

Data availability
The datasets generated in this study are available in the Supplementary Information files (available online). The sequence data were submitted to the NCBI Sequence Read Archive (BioProject accession: PRJNA669279) or are available from the corresponding author (Dr. Yu-Shin Nai) by request.