Skip to main content

DNA methylomes and transcriptomes analysis reveal implication of host DNA methylation machinery in BmNPV proliferation in Bombyx mori



Bombyx mori nucleopolyhedrosis virus (BmNPV) is a major pathogen that threatens the sustainability of the sericultural industry. DNA methylation is a widespread gene regulation mode in epigenetics, which plays an important role in host immune response. Until now, little has been known about epigenetic regulation on virus diseases in insects. This study aims to explore the role of DNA methylation in BmNPV proliferation.


Inhibiting DNA methyltransferase (DNMT) activity of silkworm can suppress BmNPV replication. The integrated analysis of transcriptomes and DNA methylomes in silkworm midguts infected with or without BmNPV showed that both the expression pattern of transcriptome and DNA methylation pattern are changed significantly upon BmNPV infection. A total of 241 differentially methylated regions (DMRs) were observed in BmNPV infected midguts, among which, 126 DMRs were hyper-methylated and 115 DMRs were hypo-methylated. Significant differences in both mRNA transcript level and DNA methylated levels were found in 26 genes. BS-PCR validated the hypermethylation of BGIBMGA014008, a structural maintenance of chromosomes protein gene in the BmNPV-infected midgut. In addition, DNMT inhibition reduced the expression of inhibitor of apoptosis family genes, iap1 from BmNPV, Bmiap2, BmSurvivin1 and BmSurvivin2.


Our results indicate that DNA methylation plays positive roles in BmNPV proliferation and loss of DNMT activity could induce the apoptosis of infected cells to suppress BmNPV proliferation. Our results may provide a new idea and research direction for the molecular mechanism on insect-virus interaction.


DNA methylation, as one of the important epigenetic regulations, occurs in both eukaryotes and prokaryotes. Accumulating evidences have linked epigenetic effects to various biological processes including gene regulation, development, nutrigenomics, tumorigenesis, and DNA repair in mammals and plants [1,2,3]. Recently, the roles of DNA methylation in host immune responses have attracted increasing attention. Many studies have demonstrated that viruses, especially DNA tumor viruses, could induce aberrant DNA methylation of host genome to suppress the host immune responses and evade antiviral immunity [4,5,6]. In addition, viruses could modulate host DNA methyltransferase (DNMT) for epigenetic dysregulation of immune-related gene expression in host cells [7, 8]. Interestingly, demethylation treatment of cancer cells could activate the virus RNA transcription from dormant endogenous retroviruses and trigger antiviral signaling [9, 10].

Although comparing to mammals, the researches on DNA methylation in insects are relatively lagging behind, with the improvement and innovation of DNA methylated research methods, especially for the rapid progress in large-scale sequencing technology, some progresses and achievements on DNA methylation in insects have been made in recent years. Based on whole genome bisulfite sequencing (WGBS), insects such as Drosophila melanogaster, Tribolium castaneum, Bombyx mori and Nasonia Vitripennis [11,12,13,14] have proven that DNA methylation in insect genome does exist and the maintenance and regulation mechanism of DNA methylation system in insects may greatly differ from that of mammals and plants [15,16,17]. The function of DNA methylation generally associates with the evolution, aging, memory and regulation of caste determination in honey bees and in ants [18,19,20,21].

Until now, researches on the function of DNA methylation in insect’s immune responses are very limited. It is reported the genome wide-patterns of DNA methylation in the Aedes aegypti are disrupted with the infection of a virulent strain of Wolbachia [22]. Three different strains of Metarhizium anisopliae caused the differential expression of dnmt genes in the larvae of Galleria mellonella infected in a natural manner, suggesting that DNA methylation may play roles in the immune response of insects against parasitic fungi [17]. In Bombyx mori only two DNMT have been reported, DNMT1 and DNMT2. BmDNMT-1 retained the function as maintenance DNMT [23].. Our previous study found that 27 genes were shown to have both differential expression and differential methylation in the midgut and fat body of Bombyx mori cytoplasmic polyhedrosis virus (BmCPV) infected larvae, respectively, indicating that the BmCPV infection-induced expression changes of these genes could be mediated by variations in DNA methylation [24].

Bombyx mori nucleopolyhedrosis virus (BmNPV), a circular double-stranded DNA (dsDNA) virus, belongs to the family Baculoviridae [25] and is a big challenge for the development and sustainability of the sericultural industry [26]. Early study has found that in Autographa calfomnica nuclear polyhedrosis virus (AcNPV), the activity of the p10 gene promoter was severely affected when the 5′-CCGG-3′ site in this promoter was methylated, suggesting that methylation of specific promoter sequences in an insect virus can affect viral gene expression [27]. Up to now, there has been scarce report on epigenetic function on BmNPV infection and immune response in Bombyx mori.

In this study, we treated the Bombyx mori cells with Zebularine (Zeb), a kind of DNMT inhibitor and found that BmNPV replication is significantly suppressed. Furthermore, genome-wide methylome and transcriptome analyses were carried out to explore the possible role of DNA methylation in silkworm immune response and BmNPV-host interaction.


DNA methyltransferase inhibition suppresses the replication of BmNPV

Zeb has been well known as an inhibitor of DNMT. We first tested the cytotoxicity of Zeb for BmN cells. BmN cells were treated with different concentrations of Zeb for 12 h followed by MTT assay. The result showed that there were no significant differences in cell survival among 20 μM, 50 μM, 100 μM Zeb-treated cells and control cells (Fig. 1a). The viability of 150 μM treated cells only accounted for 80.17% in control cells. As the Zeb concentration increased to 200 μM, the viability was less than 50% compared to that of control cells. This suggested that less than 100 μM Zeb has no significant harmful effects on BmN cell survival.

Fig. 1

Determination of the cytotoxicity of zebularine and the chemicals effect on BmNPV proliferation. a BmN cells were treated with different concentrations of zebularine (20, 50, 100, 150, 200, 250 μM) for 12 h. Cytotoxicity was measured by MTT assay. Cellular cytotoxicity was determined in duplicate and each experiment was repeated three times. ** p < 0.01.b BmN cells were treated with 100 μM zebularine and infected with recombinant BmNPV BVs containing an egfp marker gene (Zeb/BmNPV). DMSO-treated BmN cells with recombinant BmNPV infection were used as the control (DMSO/BmNPV). Fluorescence intensity was observed under fluorescence microscopeat 72 hpi. c RNAs were extracted from BmNPV-infected BmN cells treated with zebularine and DMSO, respectively. Absolute qRT-PCR was carried out to analyze the copies of ie-1 and lef-3 at different time points. The data were represented as mean ± SD. Three independent experiments with three technical replicates were performed. d Western blotting analysis of VP39. Protein samples were from BmNPV-infected BmN cells with zebularine and DMSO, respectively. Western blot analysis for detection of VP39 was performed by an anti-VP39 antibody. a-Tubulin was served as the loading control

In order to examine whether DNA methyltransferase inhibition affects BmNPV replication, BmN cells were then treated with 100 μM Zeb for 12 h and then replaced with a new medium followed by infection of recombinant BmNPV BVs containing an egfp marker gene for 72 h. By observing under the fluorescence microscope we observed that both the amount of cells with fluorescence and the fluorescence intensity of the treated group (Zeb/BmNPV) were significantly less than that of the control group (DMSO/BmNPV) (Fig. 1b). Subsequently, qRT-PCR was carried out to detect the expression of BmNPV ie-1 gene and lef-3 gene at different time points. Comparing the results of Zeb-treated cells to control cells, we found that the expression level of these two genes were very similar with the tendency of increase at 6 h while decrease at 12 h–48 h in Zeb-treated cells (Fig. 1c). Finally, western blot was performed to assess the Zeb impact on BmNPV replication (Fig. 1d). All of these results suggested that inhibition of DNMT could suppress BmNPV proliferation.

BmNPV infection alters the transcriptional profiles of silkworm midgut

By RNA-Seq technology, we analyzed the transcriptional profiles of normal and infected silkworm midgut. More than 43,000,000 clean reads were obtained from each sample and the summary of the sequenced data is shown in Additional file 1: Table S1. The sequence data are deposited in the NCBI Sequence Read Archive database ( under the accession number PRJNA541379. The number of significantly differentially expressed genes (DEGs) in infected midgut was 2171, with 920 up-regulated genes and 1251 down-regulated genes (Fig. 2a, Additional file 1: Table S2). GO analysis revealed that the down-regulated genes were enriched in ATP metabolic process, cytoplasm, catalytic activity, small molecule metabolic process, etc. Up-regulated genes were enriched in DNA binding, chromosome, DNA replication and so on (Fig. 2b). KEGG pathway analysis showed that up-regulated genes were involved in DNA replication, base excision repair, RNA transport, spliceosome and so on while down-regulated genes were related to valine, leucine and isoleucine degradation, fatty acid degradation, metabolic pathways, and protein processing in endoplasmic reticulum, etc. (Additional file 2: Figure S1).

Fig. 2

Differentially expressed genes associated with BmNPV infection. a The volcano plot of differentially expressed genes. Red points represent up-regulated genes; Green points represent down-regulated genes. Blue points represent non-changed genes. b.GO anotation was performed by mapping genes to GO terms in the GO database ( GO enrichment analysis was conducted by GOseq R package with corrected p-value ≤0.05 as a threshold. Gene numbers are listed for each category

Overview of whole genome bisulfite sequencing

To evaluate epigenetic changes in insect cells caused by BmNPV infection, genome-wide DNA methylation profiles of BmNPV-infected midgut (T group) and normal midgut (C group) were performed at single-base resolution. Based on more than 99.8% bisulfite conversation rate and more than 10× mean sequencing coverage on cytosine site, we observed that about 0.185% of genomic cytosine was methylated, which was higher than 0.11% in the silk glands of silkworm (Xiang, 2010). The average methylation level at CG location was about 0.875%. The global difference methylation level between T group and C group was displayed by Circos software, which can visualize data in circular layout [28](Fig. 3a). Relatively, on scaffold 4 and 11, some methylated sites were shown hyper-methylation while on scaffold 5 and 19; a lot of hypo-methylated sites were observed (T vs. C). Furthermore, we comparatively analyzed the methylation level on different genomic functional regions between these two groups and the results revealed that in silkworm, DNA methylation was mostly targeted to gene bodies. It peaked in exon region and sharply decreased in intron region. Promoter region showed more methylation than the repeat region (Fig. 3b). Fractional methylation levels in exon region exhibited the tendency of hypo-methylation in BmNPV-infected midguts than in control midguts (Fig. 3b).

Fig. 3

The global difference methylation level between T and C groups on CG site. a Circos images were used to show the different methylation levels between C group (normal midgut) and T group (BmNPV-infected midgut) at the genome-wide scale. Each circle from the outside to the inside depicts: (a) 20 scaffolds of Bombyx mori from the largest to the smallest size distributed in clockwise; (b) the methylation level of T group; (c) difference methylation level between T groups and C group; (d) the methyaltion level of C group. Color scale from green to blue represents DNA methylation level; Color scale from blue to red represents DNA methylation difference between two groups. b Distribution of DNA methylation levels of genes in T and Cgroup. The gene structure is defined by seven different features, denoted by the x-axis. The length of each feature was normalized and divided into equal numbers of bins

Differential DNA methylation regions (DMRs) upon BmNPV infection

To explore the potential function of DNA methylation in response to BmNPV infection, we compared the DMR of the whole genome between BmNPV infected midgut and the control tissues. The results revealed that 241 DMRs were obtained with 177 located in the gene region. Among these, 126 DMRs were hyper-methylated and 115 were hypo-methylated (Additional file 1: Table S3). Further study revealed that DMRs are preference for the intron region (Fig. 4a). KEGG pathway analysis of 177 differential methylation genes (DMGs) revealed that DMGs were involved in pathways of spliceosome, RNA transport, protein processing in endoplasmic reticulum, etc. (Fig. 4b).

Fig. 4

Differentially methylated genes associated with BmNPV infection. a The number of differentially methylated genes in different regions of the genome. b KEGG pathway enrichment of differentially methylated genes. RichFactor is the ratio of the number of differentially expressed genes in this pathway term to the number of all genes in this pathway term. Q-value is corrected for ranging from 0 to 1. Only the top 20 of the enriched pathway terms are displayed here

Effects of DMRs on gene expression in silkworm challenged by BmNPV infection

To investigate whether the variances in CG methylation upon BmNPV infection could alter gene expression, we comprehensively analyzed the data of DEGs and DMRs. As a result, we found 26 DMGs showing differential mRNA levels in the midgut of the infected silkworm (Table 1), suggesting that the expression levels of these genes involved in BmNPV infection may altered via variation in DNA methylation. Analysis of the distribution of DNA methylation on these genes revealed that CG methylation of 25 genes presented in the gene body, predominantly in the introns of genes, accounted for about 61.54% (Additional file 3: Figure S2). Interestingly, there were 2 genes (BGIBMGA002246 and BGIBMGA004695), whose intron-exon boundaries were DNA methylated, indicating that DNA methylation may be associated with alternative splicing in Bombyx mori [29]. Furthermore, by qRT-PCR and BS-PCR, we validated the up-regulation (Fig. 5a) and hyper-methylation of BGIBMGA014008 (Fig. 5b, c), encoding for structural maintenance of chromosomes protein in BmNPV-infected midgut.

Table 1 Gene informatin of both DEGs and DMRs in silkworm midguts with BmNPV infection
Fig. 5

Validation of both different expression and different methylation of BGIBMGA014008 in response to BmNPV infection. a BGIBMGA014008 was identified to be differentially expressed upon BmNPV infection by qRT-PCR. RNA was extracted from midguts with or without BmNPV infection. Data represent the relative transcript levels of genes in the infected midguts compared to the control tissues from three independent samples. The error bars indicate standard deviations (*p < 0.05). C group represents samples from normal midguts and T group representssamples from BmNPV-infected midguts. b The visualized image of differential methylation of BGIBMGA014008was generated by Integrative Genomics Viewer (IGV). The location of blue bars in the images indicates the methylated sites and the height of blue bars represents the methylation levels. c Bisulfite sequencing validation of differential methylation of BGIBMGA014008in normal verse infected midgut. Targeted bisulfite sequencing validation of a region overlapping the DMR is indicated by the rectangle over the genome browser tracks. Dark circles indicate methylated and open circles mean unmethylated cytosines. Each row consists of a single sequenced clone

DNMT inhibition decreases the expression of antiapoptosis-related genes in Bombyx mori

To explore the potential mechanism of DNA methylation on regulating host immune response to BmNPV infection, we detected the expressed level of 7 genes, which are associated with JAK/STAT pathway (stat, socs6 and socs2) and anti-apoptosis (iap1 from BmNPV, Bmiap2, BmSurvivin1, BmSurvivin2) between Zeb-treated cells and control cells. The results showed that the transcript level of socs6 is up-regulated and socs2 is down-regulated while that of stat has no significant change (Fig. 6). Interestingly, the expression levels of all the four genes related to anti-apoptosis are significantly down-regulated (Fig. 6), indicating that DNMT inhibition may promote the apoptosis of infected cells by suppressing the expression of anti-apoptosis genes both in virus and in Bombyx mori.

Fig. 6

DNMT inhibition alters the expression of genes involving in antiapoptosis. RNAs were extracted from BmNPV-infected BmN cells treated with Zeb and DMSO, respectively and then transcribed into cDNA as template for RT-qPCR. The data were normalized using BmGAPDH and are represented as mean ± SD from three independent experiments. Relative expression levels were calculated using the 2−ΔΔCt method(**p < 0.01)


The role of epigenetic regulation of gene expression has been widely accepted in animal and plants with the change of environmental stressors, such as the invasion of pathogens. Recent studies have shown that inhibition of DNA methylation by chemicals such as 5-aza-2′-deoxycytidine significantly induced antitumor immune responses in colon and ovarian cancers [9, 10] and resulted in cell cycle arrest, p53-dependent apoptosis, and IFN signaling activation of HPV-positive HNC cells [30]. In this study, we found that inhibition of DNMT by Zeb could suppress BmNPV replication.

To explore whether the BmNPV infection can induce the variance on genomic DNA methylation of silkworm and whether DNA methylation could regulate the expression of genes associated with BmNPV infection, we performed transcriptome analyses and whole-genome bisulfite sequencing between BmNPV-infected midguts and control midguts of silkworm larvae. As a result, the characteristics of genomic DNA methylation in Bombyx mori are similar to previous studies on other insects and distinct from that of animal and plant. For insects, the general methylation level is much lower than that in mammals. At current, the cytosine methylation level of most insects is about 0 ~ 1%, while in mammals and birds’ ranges from 3 to 10%, fish and amphibians around 10%, and some plants up to 50% [28]. In most mammals, about 60 to 90% of all CpG dinucleotides are methylated. DNA methylation in vertebrates distributes in the entire genome [31] while due to the low level of methylation in insects, DNA methylation in insects is scattered throughout the genome and is mainly concentrated in gene body (introns and exons). Promoters and repeat region DNAs in insects are hypo-methylated [32]. In this study, as to Bombyx mori, it exhibits low levels of DNA methylation with about 0.185% and CG methylation is substantially enriched in gene bodies.

The function of DNA methylation on gene expression is well established, but the underlying mechanisms of this function are poorly understood. It is speculated that cytosine methylation can facilitate or suppress interaction with DNA-binding protein, change the stability of nucleosomes, affect the local chromatin structure and access of transcript factors to genomic DNA, and thus result in changes of gene expression [33, 34]. In this study, we found 26 DEGs showing differential methylated levels in the midgut of BmNPV-infected silkworm larvae with 17 hyper-methylated and 9 hypo-methylated, which demonstrated that variance in DNA methylation may lead to expression changes of those genes associated with BmNPV infection.

BGIBMGA014008 is a structural maintenance of chromosomes protein gene. Studies showed that structural maintenance of chromosomes 4 (SMC-4) is a chromosomal ATPase which plays an important role in regulating chromosome assembly and segregation. SMC-4 expression was significantly higher in colorectal cancer and was associated with tumorigenesis. Knockdown of SMC-4 significantly suppressed the proliferation of cancer cells and degraded its malignant degree [35]. In addition, SMC4 expression was significantly associated with tumor size, dedifferentiation, advanced stages and vascular invasion of the primary liver cancers [36]. In our study, BGIBMGA014008 was identified up-regulation and hypermethylation in the midgut after BmNPV infection (Fig. 5), indicating that BGIBMGA014008 may involve in silkworm cell multiplication and play roles in BmNPV infection.

Dihydroorotate dehydrogenase (Dhodh), catalyzing the oxidation of dihydroorotate to orotate, is the fourth enzyme of pyrimidine synthesis in the de novo pyrimidine biosynthetic pathway [37]. It is reported that down-regulation of BmDhodh inhibits cell growth and the endomitotic DNA replication in silk gland cells [38]. In our study, we observed the increased mRNA level and hypo-methylation in promoter region of BmDhodh (BGIBMGA011887). Generally, promoter methylation represses gene expression [39]. We speculated that hypomethylation in the promoter region of BmDhodh activates BmDhodh expression, which is benefit for silkworm cell replication, and thus facilitates BmNPV proliferation.

Nuclear import and export of viral RNA and proteins are critical to the replication cycle of viruses. Regulation of this process is paramount to processes such as cell division and differentiation, but is also critically important for viral replication and pathogenesis [40]. Nucleocytoplasmic transport is mediated by nuclear pore complexes (NPCs), which are embedded in the nuclear envelope [41]. Nup188, one of the nucleoporins can regulate chromosome segregation, promote chromosome alignment, confine the passage of membrane proteins and is thus crucial for the homeostasis of the different nuclear membranes [42, 43]. In this study, both the mRNA transcript level of Nup188 (BGIBMGA002694) and DNA methylation level in the intron of Nup188 were significantly enhanced in BmNPV infected midgut, suggesting that DNA methylation may have impact on nucleoporins expression and thus affect the replication of BmNPV.

Other interesting genes, such as Rad50 (BGIBMGA005449), which is associated with DNA repair functions [44,45,46] and the histidine triad protein gene (BGIBMGA011899), which plays roles in transcription, signal transduction and many peripheral and central nervous system diseases [47] are also shown different mRNA transcript level and DNA methylated level upon BmNPV infection (Table 1).

Recent studies have revealed that gene expression regulation by DNA methylation may play a critical role in arms races between viruses and their hosts [48]. To evade detection and restriction by the host immune response, viruses also employ various mechanisms to control gene expression related to immunity, including hijacking epigenetic machinery [49, 50]. Studies on virus-driven dysregulation of host immune-related gene expression through DNA methylation present a novel viral mechanism to inhibit immune responses [51]. Many DNA and RNA virus could utilize this mechanism to down-regulates expression of host immune-related genes [6, 52, 53].

Previous reports revealed that BmNPV infection could activate the JAK/STAT signaling pathway and has slight effects on Imd and Toll signaling pathways [54,55,56]. The JAK signal transducer and STAT signaling pathway is an important regulator of cell proliferation, differentiation, survival, apoptosis and immune response [57]. Therefore, to explore whether DNA methylation could affect the JAK/STAT pathway, we detected the expression of three key genes involved in this pathway and we found that after inhibition of DNMT, two genes, which are negative regulators of the JAK/STAT signaling pathway [58], displayed the opposite expression pattern with up-regulation of socs6 and down-regulation of socs2. Another key gene, stat has no significant expression changes (Fig. 6), indicating that DNA methylation may have no obvious effects on JAK/STAT pathway.

Apoptosis as an important immune response plays a critical role in the antiviral defense in insects [59]. Inhibitor of apoptosis (IAP) protein family have demonstrated functions in both apoptosis and innate immunity [60]. Both BmIAP and Baculoviruses encode IAP had an inhibitory effect on apoptosis in insects [61, 62]. Survivin, another apoptosis inhibitor has been confirmed to be an essential regulator of cell division [63]. In our study, we found four genes belong to inhibitor of apoptosis family were all decreased their expression (Fig. 6), suggesting that DNMT1 inhibition may lead to promote the apoptosis of infected cells and consequently suppress the BmNPV proliferation.


In summary, based on parallel analyses of global gene expression and the cellular methylome upon BmNPV infection, our results suggest that alterations in DNA methylation status may have effects on the expression of genes associated with BmNPV infection. Inhibition of DNMT in silkworm cells could induce the apoptosis of infected cells. Our results may provide a new idea and research direction for the molecular mechanism of the interaction between silkworm and viruses.


Cell line, silkworm strain and virus inoculation

The Bombyx mori cell line BmN, which was kindly gifted by Dr. Xudong Tang from Silkworm Pathology Laboratory in Jiangsu University of Science and Technology, was cultured at 27 °C in TC-100 medium (United States Biological, USA) supplemented with 10% (V/V) fetal bovine serum (FBS) (Gibco, USA) in 6-well plates to a confluency of about 70%. BmN cells per well were infected with recombinant BmNPV BVs containing an egfp marker gene at a MOI of 5. The P50 silkworm strain was reared at room temperature and under a photoperiod of 12 h light and 12 h dark until the fourth molt. For viral inoculation, 2 μL BmNPV viral stock including recombinant BmNPV BVs was injected into each larva through intersegmental membrane. The control uninfected larvae were injected with 2 μL 0.9% NaC1 solution. The infected group and the uninfected group had 3 replicates respectively and each replicate was pooled with 5 larvae.

DNMT inhibitor and cellular cytotoxicity assay (MTT)

BmN cells were seeded in 96 well plates at a density of about 5000 cells per well and were treated with different concentrations of DNMT inhibitor, Zeb (MedChemExpress) for 12 h. 10 μl MTT (5 mg/ml, Sigma) was added to each well at 37 °C in 5% CO2 for 4 h. Later, the medium was removed carefully without disturbing the formazan crystals and 150 ml of dimethyl sulfoxide (DMSO) was added followed by incubation at 37 °C for 15 min. The absorbance of the suspension was measured at 490 nm. The percentages of metabolically active cells were compared with the percentage of control cells of the same culture plate. Cellular cytotoxicity was determined in duplicate and each experiment was repeated three times independently.

Quantitative real-time PCR analysis

Relative quantitative real-time PCR (qRT-PCR) as descried in our previous study was performed to detect the expression of silkworm gene [64, 65]. As to BmNPV gene lef-3 (GeneID: 1488686) and ie-1 (GeneID: 1488755), absolute qRT-PCR was carried out as in our previous report [64]. PCR reactions were run with a thermal cycling at 95 °C for 30 s followed by 40 cycles of 95 °C for 5 s, 60 °C for 31 s. Following the amplification, melting curves were constructed. Three independent experiments with three technical replicates were analyzed. All data were represented by the mean ± SD. The unpaired Student’s t test was used to compare the difference in means. P value < 0.05 was set for statistically significant. The sequences of the primers were listed in Additional file 1: Table S4.

Western blot assay

Western blotting assay was performed as descried in our previous study [64]. Briefly, a total of 30 μg protein sample from cells was loaded on each lane and separated on SDS-PAGE before transferred onto a nitrocellulose membrane. The membrane was incubated with rabbit VP39 (1:2000). The same blots were re-probed with ɑ-Tubulin ployclonal antibody (1:2000; Sigma, USA) to confirm equal loading of samples. The secondary antibody of VP39 and ɑ-tubulin is HRP labeled goat anti-rabbit IgG (1:1000; Beyotime, China).

Transcriptome analysis

RNA-Seq experiment with three biological replicates was performed by Novogene Bioinformatics Technology Co., LTD (Beijing, China). The process is described briefly as follows: total RNA was extracted by using the Trizol reagent (Invitrogen, USA) and RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer’s recommendations. The library fragments of preferentially 250~300 bp in length were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then, the fragments were ligated to sequencing adaptors and enriched by PCR amplification. The libraries were sequenced on an Illumina Hiseq platform after library quality was assessed on the Agilent Bioanalyzer 2100 system.

Clean reads, which were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data were mapped to silkworm genomic database(B. mori assembly ASM15162v1) using Hisat2 v2.0.4. Differentially expressed genes were identified by the DESeq R package (1.18.0). The resulting P-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 and fold change ≥2 were assigned as differentially expressed. Gene Ontology (GO) enrichment analysis of differentially expressed genes was performed by the GOseq R package. GO terms with corrected P-value less than 0.05 were considered significantly enriched by differential expressed genes. KOBAS software was used to test the statistical enrichment of differential expression genes in KEGG pathways ( [66].

Whole genome bisulfite sequencing

For whole-genome bisulfite sequencing, three biological replicates were performed. A total amount of 5.2 mg genomic DNA were fragmented by sonication to 200-300 bp, followed by end repair, adenylation and ligating methylated adaptors. Then these DNA fragments were treated twice with bisulfite using EZ DNA Methylation-Gold™ Kit (Zymo Research) before PCR amplification. Finally, sequencing was performed on an Illumina Hiseq 4000 platform and 125 bp/150 bp paired-end reads (raw reads) were generated. The clean reads were obtained from raw reads after pre-processed through Trimmomatic (Trimmomatic-0.36) software and all filtering steps.

Bismark software (version 0.16.3) [67] was used to perform alignments of bisulfite-treated reads to silkworm genome (B. mori assembly ASM15162v1). Firstly, silkworm genome was changed into bisulfite-converted version (C-to-T and G-to-A) followed by indexing using bowtie2 [68]. Then, sequence reads were changed into bisulfite-converted versions before aligned to converted genome in a directional manner. The unique best alignment reads are then compared to the normal genomic sequence and the methylation level of all cytosine in the reads is evaluated.

The results of methylation extractor were transformed into bigWig format for visualization using IGV browser. In order to calculate the methylation level of the sequence, we divided the sequence into multiple bins with bin size is 10 kb. The sum of methylated and unmethylated read counts in each window was calculated. Methylation level (ML) for each C site shows the fraction of methylated Cs, and is defined as:


Calculated ML was further corrected based on bisulfite non-conversion rate and sequencing coverage depth on a certain site. To obtain the accurate mC sites, threshold values were set up as follows: sequencing depth is≥5, q-value is≤0.01.

Differentially methylated analysis

Differentially methylated regions (DMRs) were identified using the DSS software (version DSS_2.12.0) [8, 65, 69] . The core of DSS is a new dispersion shrinkage method for estimating the dispersion parameter from Gamma-Poisson or Beta-Binomial distributions. Putative DMRs were identified based on the following standards and parameters: (1) in the smoothing process, the smoothing distance is 200 bp (parameter: smoothing = TRUE, smoothing.span = 200, delta = 0); (2) the proportion of the difference sites with P value less than 1e-05 was greater than 50% of the region (parameter: p.threshold = 1e-05; pct.sig = 0.5); (3) the number of the sites was greater than 3, and the length was greater than 50 (parameter: minlen = 50; minCG = 3); (4) when the distance between two DMR is less than 100 bp, the two regions are merged (parameter: dis.merge = 100). GO and KEGG pathway analysis of genes containing DMRs were performed as described in the “Transcriptome analysis” section.

Bisulfite-PCR validation of DMRs

Genomic DNA from BmNPV-infected and control tissues was extracted and treated with bisulfite as described above. Bisulfite converted DNA was uses for PCR amplification with ZymoTaqTN Preix DNA Polymerase (ZYMO RESEARCH, America). Nested primers for BS-PCR were designed using the online MethPrimer program (Additional file 1: Table S4). PCR products were purified and cloned into the pMD19-T vector (TaKaRa, Japan). Ten different clones for each group were sent to Sangon Biotech Co., Ltd. (Shanghai, China) for sequencing. DNA methylation of individual CG sites were analyzed using Quma software (

Availability of data and materials

The raw data in this study are deposited in the NCBI Sequence Read Archive database ( under the accession number PRJNA541379.



Bombyx mori nucleopolyhedrosis virus


Differentially expressed genes


Differential methylation genes


Differential DNA methylation regions


Dimethyl sulfoxide


DNA methyltransferase


Inhibitor of apoptosis


Whole genome bisulfite sequencing




  1. 1.

    Park J, Peng Z, Zeng J, Elango N, Park T, Wheeler D, Werren JH, Yi SV. Comparative analyses of DNA methylation and sequence evolution using Nasonia genomes. Mol Biol Evol. 2011;28(12):3345–54.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  2. 2.

    Baylin SB, Jones PA. A decade of exploring the cancer epigenome - biological and translational implications. Nat Rev Cancer. 2011;11(10):726–34.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  3. 3.

    Riccio A, Alvania RS, Lonze BE, Ramanan N, Kim T, Huang Y, Dawson TM, Snyder SH, Ginty DD. A nitric oxide signaling pathway controls CREB-mediated gene expression in neurons. Mol Cell. 2009;21(2):283–94.

    Article  CAS  Google Scholar 

  4. 4.

    Cicchini L, Blumhagen RZ, Westrich JA, Myers ME, Warren CJ, Siska C, Raben D, Kechris KJ, Pyeon D. High-risk human papillomavirus E7 alters host DNA Methylome and represses HLA-E expression in human keratinocytes. Sci Rep. 2017;7(1):3633.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  5. 5.

    Dong SM, Lee HG, Cho SG, Kwon SH, Yoon H, Kwon HJ, Lee JH, Kim H, Park PG, Kim H, et al. Hypermethylation of the interferon regulatory factor 5 promoter in Epstein-Barr virus-associated gastric carcinoma. J Microbiol. 2015;53(1):70–6.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  6. 6.

    Di Bartolo DL, Cannon M, Liu YF, Renne R, Chadburn A, Boshoff C, Cesarman E. KSHV LANA inhibits TGF-beta signaling through epigenetic silencing of the TGF-beta type II receptor. Blood. 2008;111(9):4731–40.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  7. 7.

    Kang MS, Kieff E. Epstein-Barr virus latent genes. Exp Mol Med. 2015;47:e131.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  8. 8.

    Wu J, Xu Y, Mo D, Huang P, Sun R, Huang L, Pan S, Xu J. Kaposi's sarcoma-associated herpesvirus (KSHV) vIL-6 promotes cell proliferation and migration by upregulating DNMT1 via STAT3 activation. PLoS One. 2014;9(3):e93478.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  9. 9.

    Roulois D, Loo Yau H, Singhania R, Wang Y, Danesh A, Shen SY, Han H, Liang G, Jones PA, Pugh TJ, et al. DNA-Demethylating agents target colorectal Cancer cells by inducing viral mimicry by endogenous transcripts. Cell. 2015;162(5):961–73.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  10. 10.

    Chiappinelli KB, Strissel PL, Desrichard A, Li H, Henke C, Akman B, Hein A, Rote NS, Cope LM, Snyder A, et al. Inhibiting DNA methylation causes an interferon response in Cancer via dsRNA including endogenous retroviruses. Cell. 2015;162(5):974–86.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. 11.

    Capuano F, Mulleder M, Kok R, Blom HJ, Ralser M. Cytosine DNA methylation is found in Drosophila melanogaster but absent in Saccharomyces cerevisiae, Schizosaccharomyces pombe, and other yeast species. Anal Chem. 2014;86(8):3697–702.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  12. 12.

    Xiang H, Zhu J, Chen Q, Dai F, Li X, Li M, Zhang H, Zhang G, Li D, Dong Y, et al. Single base-resolution methylome of the silkworm reveals a sparse epigenomic map. Nat Biotechnol. 2010;28(5):516–20.

    CAS  PubMed  Article  Google Scholar 

  13. 13.

    Feliciello I, Parazajder J, Akrap I, Ugarkovic D. First evidence of DNA methylation in insect Tribolium castaneum: environmental regulation of DNA methylation within heterochromatin. Epigenetics. 2013;8(5):534–41.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Beeler SM, Wong GT, Zheng JM, Bush EC, Remnant EJ, Oldroyd BP, Drewell RA. Whole-genome DNA methylation profile of the jewel wasp (Nasonia vitripennis). G3 (Bethesda). 2014;4(3):383–8.

    CAS  Article  Google Scholar 

  15. 15.

    Glastad KM, Hunt BG, Yi SV, Goodisman MA. DNA methylation in insects: on the brink of the epigenomic era. Insect Mol Biol. 2011;20(5):553–65.

    CAS  PubMed  Article  Google Scholar 

  16. 16.

    Standage DS, Berens AJ, Glastad KM, Severin AJ, Brendel VP, Toth AL. Genome, transcriptome and methylome sequencing of a primitively eusocial wasp reveal a greatly reduced DNA methylation system in a social insect. Mol Ecol. 2016;25(8):1769–84.

    CAS  PubMed  Article  Google Scholar 

  17. 17.

    Wu P, Jie W, Shang Q, Annan E, Jiang X, Hou C, Chen T, Guo X. DNA methylation in silkworm genome may provide insights into epigenetic regulation of response to Bombyx mori cypovirus infection. Sci Rep. 2017;7(1):16013.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  18. 18.

    Bonasio R, Li Q, Lian J, Mutti NS, Jin L, Zhao H, Zhang P, Wen P, Xiang H, Ding Y, et al. Genome-wide and caste-specific DNA methylomes of the ants Camponotus floridanus and Harpegnathos saltator. Curr Biol. 2012;22(19):1755–64.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  19. 19.

    Yan H, Bonasio R, Simola DF, Liebig J, Berger SL, Reinberg D. DNA methylation in social insects: how epigenetics can control behavior and longevity. Annu Rev Entomol. 2015;60(60):435–52.

    CAS  PubMed  Article  Google Scholar 

  20. 20.

    Biergans SD, Giovanni GC, Judith R, Charles C. DnmtsandTettarget memory-associated genes after appetitive olfactory training in honey bees. Sci Rep. 2015;5(1):21656.

    Article  CAS  Google Scholar 

  21. 21.

    Herb BR, Wolschin F, Hansen KD, Aryee MJ, Langmead B, Irizarry R, Amdam GV, Feinberg AP. Reversible switching between epigenetic states in honeybee behavioral subcastes. Nat Neurosci. 2012;15(10):1371–3.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. 22.

    Ye YH, Woolfit M, Huttley GA, Rances E, Caragata EP, Popovici J, O'Neill SL, McGraw EA. Infection with a virulent strain of Wolbachia disrupts genome wide-patterns of cytosine methylation in the mosquito Aedes aegypti. PLoS One. 2013;8(6):e66482.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  23. 23.

    Marhold J, Rothe N, Pauli A, Mund C, Kuehle K, Brueckner B, Lyko F. Conservation of DNA methylation in dipteran insects. Insect Mol Biol. 2004;13(2):117–23.

    CAS  PubMed  Article  Google Scholar 

  24. 24.

    Vilcinkas A. The role of epigenetics in host–parasite coevolution: lessons from the model host insects galleria mellonella and Tribolium castaneum. Zoology. 2016;119(4):273–80.

    Article  Google Scholar 

  25. 25.

    Rahman MM, Gopinathan KP. Systemic and in vitro infection process of Bombyx mori nucleopolyhedrovirus. Virus Res. 2004;101(2):109–18.

    CAS  PubMed  Article  Google Scholar 

  26. 26.

    Jiang L, Xia Q. The progress and future of enhancing antiviral capacity by transgenic technology in the silkworm Bombyx mori. Insect Biochem Mol Biol. 2014;48(1):1–7.

    CAS  PubMed  Article  Google Scholar 

  27. 27.

    Knebel D, Lubbert H, Doerfler W. The promoter of the late p10 gene in the insect nuclear polyhedrosis virus Autographa californica: activation by viral gene products and sensitivity to DNA methylation. EMBO J. 1985;4(5):1301–6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  28. 28.

    Krauss V, Eisenhardt C, Unger T. The genome of the stick insect Medauroidea extradentata is strongly methylated within genes and repetitive DNA. PLoS One. 2009;4(9):e7223.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  29. 29.

    Regulski M, Lu Z, Kendall J, Donoghue MT, Reinders J, Llaca V, Deschamps S, Smith A, Levy D, McCombie WR, et al. The maize methylome influences mRNA splice sites and reveals widespread paramutation-like switches guided by small RNA. Genome Res. 2013;23(10):1651–62.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  30. 30.

    Biktasova A, Hajek M, Sewell A, Gary C, Bellinger G, Deshpande HA, Bhatia A, Burtness B, Judson B, Mehra S, et al. Demethylation therapy as a targeted treatment for human papillomavirus-associated head and neck Cancer. Clin Cancer Res. 2017;23(23):7276–87.

    CAS  PubMed  Article  Google Scholar 

  31. 31.

    Okamura K, Matsumoto KA, Nakai K. Gradual transition from mosaic to global DNA methylation patterns during deuterostome evolution. BMC Bioinformatics. 2010;11(Suppl 7):S2.

    PubMed  PubMed Central  Google Scholar 

  32. 32.

    Zemach A, McDaniel IE, Silva P, Zilberman D. Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science. 2010;328(5980):916–9.

    CAS  PubMed  Article  Google Scholar 

  33. 33.

    Gebhard C, Benner C, Ehrich M, Schwarzfischer L, Schilling E, Klug M, Dietmaier W, Thiede C, Holler E, Andreesen R, et al. General transcription factor binding at CpG islands in normal cells correlates with resistance to de novo DNA methylation in cancer cells. Cancer Res. 2010;70(4):1398–407.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  34. 34.

    Dantas Machado AC, Zhou T, Rao S, Rastogi C, Bussemaker HJ, Rohs R. 22 evolving insights on how cytosine methylation affects protein-DNA binding. Brief Funct Genomics. 2015;14(1):61–73.

    CAS  PubMed  Article  Google Scholar 

  35. 35.

    Feng XD, Song Q, Li CW, Chen J, Tang HM, Peng ZH, Wang XC. Structural maintenance of chromosomes 4 is a predictor of survival and a novel therapeutic target in colorectal cancer. Asian Pac J Cancer Prev. 2014;15(21):9459–65.

    PubMed  Article  PubMed Central  Google Scholar 

  36. 36.

    Zhou B, Yuan T, Liu M, Liu H, Xie J, Shen Y, Chen P. Overexpression of the structural maintenance of chromosome 4 protein is associated with tumor de-differentiation, advanced stage and vascular invasion of primary liver cancer. Oncol Rep. 2012;28(4):1263–8.

    CAS  PubMed  Article  Google Scholar 

  37. 37.

    Hansen M, Le JN, Johansson E, Antal T, Ullrich A, Löffler M, Larsen S. Inhibitor binding in a class 2 dihydroorotate dehydrogenase causes variations in the membrane-associated N-terminal domain. Protein Sci. 2010;13(4):1031–42.

    Article  CAS  Google Scholar 

  38. 38.

    Zhao E, Jiang X, Cui H. Bombyx mori Dihydroorotate Dehydrogenase: Knockdown Inhibits Cell Growth and Proliferation via Inducing Cell Cycle Arrest. Int J Mol Sci. 2018;19(9).

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  39. 39.

    Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13(7):484–92.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  40. 40.

    Tomoyuki H, Keizo T. Nucleocytoplasmic shuttling of viral proteins in Borna disease virus infection. Viruses. 2013;5(8):1978–90.

    Article  CAS  Google Scholar 

  41. 41.

    Andersen KR, Onischenko E, Tang JH, Kumar P, Chen JZ, Ulrich A, Liphardt JT, Weis K, Schwartz TU. Scaffold nucleoporins Nup188 and Nup192 share structural and functional properties with nuclear transport receptors. Elife. 2013;2(11):e00745.

    PubMed  PubMed Central  Article  Google Scholar 

  42. 42.

    Itoh G, Sugino S, Ikeda M, Mizuguchi M, Kanno S, Amin MA, Iemura K, Yasui A, Hirota T, Tanaka K. Nucleoporin Nup188 is required for chromosome alignment in mitosis. Cancer Sci. 2013;104(7):871–9.

    CAS  PubMed  Article  Google Scholar 

  43. 43.

    Theerthagiri G, Eisenhardt N, Schwarz H, Antonin W. The nucleoporin Nup188 controls passage of membrane proteins across the nuclear pore complex. J Cell Biol. 2010;189(7):1129–42.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. 44.

    Boswell ZK, Rahman S, Canny MD, Latham MP. A dynamic allosteric pathway underlies Rad50 ABC ATPase function in DNA repair. Sci Rep. 2018;8(1):1639.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  45. 45.

    Myler LR, Gallardo IF, Soniat MM, Deshpande RA, Gonzalez XB, Kim Y, Paull TT, Finkelstein IJ. Single-molecule imaging reveals how Mre11-Rad50-Nbs1 initiates DNA break repair. Mol Cell. 2017;67(5):891–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  46. 46.

    Henrich SM, Usadel C, Werwein E, Burdova K, Janscak P, Ferrari S, Hess D, Klempnauer KH. Interplay with the Mre11-Rad50-Nbs1 complex and phosphorylation by GSK3β implicate human B-Myb in DNA-damage signaling. Sci Rep. 2017;7:41663.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  47. 47.

    Dang YH, Liu ZW, Liu P, Wang JB. Emerging roles of Histidine triad nucleotide binding protein 1 in neuropsychiatric diseases. Zhongguo Yi Xue Ke Xue Yuan Xue Bao. 2017;39(5):705–14.

    PubMed  Google Scholar 

  48. 48.

    Warren CJ, Van Doorslaer K, Pandey A, Espinosa JM, Pyeon D. Role of the host restriction factor APOBEC3 on papillomavirus evolution. Virus Evol. 2015;1(1):vev015.

    PubMed  PubMed Central  Article  Google Scholar 

  49. 49.

    Adhya D, Basu A. Epigenetic modulation of host: new insights into immune evasion by viruses. J Biosci. 2010;35(4):647–63.

    CAS  PubMed  Article  Google Scholar 

  50. 50.

    Paschos K, Allday MJ. Epigenetic reprogramming of host genes in viral and microbial pathogenesis. Trends Microbiol. 2010;18(10):439–47.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  51. 51.

    Kussduerkop SK, Westrich JA, Pyeon D. DNA Tumor Virus Regulation of Host DNA Methylation and Its Implications for Immune Evasion and Oncogenesis. Viruses. 2018;10(2):82.

    Article  CAS  Google Scholar 

  52. 52.

    Zhang X, Justice AC, Hu Y, Wang Z, Zhao H, Wang G, Johnson EO, Emu B, Sutton RE, Krystal JH, et al. Epigenome-wide differential DNA methylation between HIV-infected and uninfected individuals. Epigenetics. 2016;11(10):750–60.

    PubMed  PubMed Central  Article  Google Scholar 

  53. 53.

    Nakayamahosoya K, Ishida T, Youngblood B, Nakamura H, Hosoya N, Koga M, Koibuchi T, Iwamoto A, Kawanatachikawa A. Epigenetic repression of interleukin-2 expression in senescent CD4+ T cells during chronic human immunodeficiency virus type-1 infection. J Infect Dis. 2015;211(1):28.

    CAS  Article  Google Scholar 

  54. 54.

    Cheng T, Lin P, Huang L, Wu Y, Jin S, Liu C, Xia Q. Genome-wide analysis of host responses to four different types of microorganisms in Bombyx Mori (Lepidoptera: Bombycidae). J Insect Sci. 2016;16(1):69.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  55. 55.

    Zhang X, Guo R, Kumar D, Ma H, Liu J, Hu X, Cao G, Xue R, Gong C. Identification, gene expression and immune function of the novel Bm-STAT gene in virus-infected Bombyx mori. Gene. 2016;577(1):82–8.

    CAS  PubMed  Article  Google Scholar 

  56. 56.

    Liu W, Liu J, Lu Y, Gong Y, Zhu M, Chen F, Liang Z, Zhu L, Kuang S, Hu X, et al. Immune signaling pathways activated in response to different pathogenic micro-organisms in Bombyx mori. Mol Immunol. 2015;65(2):391–7.

    CAS  PubMed  Article  Google Scholar 

  57. 57.

    Truong AD, Rengaraj D, Hong Y, Hoang CT, Hong YH, Lillehoj HS. Differentially expressed JAK-STAT signaling pathway genes and target microRNAs in the spleen of necrotic enteritis-afflicted chicken lines. Res Vet Sci. 2017;115:235–43.

    CAS  PubMed  Article  Google Scholar 

  58. 58.

    Seif F, Khoshmirsafa M, Aazami H, Mohsenzadegan M, Sedighi G, Bahar M. The role of JAK-STAT signaling pathway and its regulators in the fate of T helper cells. Cell Commun Signal. 2017;15(1):23.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  59. 59.

    Huipeng YAO, Xiaofeng WU, Gokulamma K. Antiviral activity in the mulberry silkworm , Bombyx mori L. Journal of Zhejiang University Science A. 2006;7(2):350–6.

    Google Scholar 

  60. 60.

    Leu JH, Kuo YC, Kou GH, Lo CF. Molecular cloning and characterization of an inhibitor of apoptosis protein (IAP) from the tiger shrimp. Dev Comp Immunol. 2008;32(2):121–33.

    CAS  PubMed  Article  Google Scholar 

  61. 61.

    Ikeda M, Yamada H, Ito H, Kobayashi M. Baculovirus IAP1 induces caspase-dependent apoptosis in insect cells. J Gen Virol. 2011;92(Pt 11):2654–63.

    CAS  PubMed  Article  Google Scholar 

  62. 62.

    Hamajima R, Iwamoto A, Tomizaki M, Suganuma I, Kitaguchi K, Kobayashi M, Yamada H, Ikeda M. Functional analysis of inhibitor of apoptosis 1 of the silkworm Bombyx mori. Insect Biochem Mol Biol. 2016;79:97–107.

    CAS  PubMed  Article  Google Scholar 

  63. 63.

    Fortugno P, Beltrami E, Plescia J, Fontana J, Pradhan D, Marchisio PC, Sessa WC, Altieri DC. Regulation of survivin function by Hsp90. Proc Natl Acad Sci U S A. 2003;100(24):13791–6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  64. 64.

    Wu P, Shang Q, Dweteh OA, Huang H, Zhang S, Zhong J, Hou Q, Guo X. Over expression of bmo-miR-2819 suppresses BmNPV replication by regulating the BmNPV ie-1 gene in Bombyx mori. Mol Immunol. 2019;109:134–9.

    CAS  PubMed  Article  Google Scholar 

  65. 65.

    Hao W, Xu T, Hao F, Li C, Li B, Bing Y, Qin Z, Peng J, Conneely KN. Detection of differentially methylated regions from whole-genome bisulfite sequencing data without replicates. Nucleic Acids Res. 2015;43(21):e141.

    Google Scholar 

  66. 66.

    Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.

    CAS  PubMed  Article  Google Scholar 

  67. 67.

    Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  68. 68.

    Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  69. 69.

    Park Y, Wu H. Differential methylation analysis for BS-seq data under general experimental design. Bioinformatics. 2016;32(10):1446–53.

    CAS  PubMed  Article  Google Scholar 

Download references


We greatly thank Dr. Xudong Tang for providing us BmN cell line.


This work was financially supported by the National Nature Science Foundation of China (No.31502016) for the design and execution of the study; Grant from the Six Talent Peak Project of Jiangsu Province, China (No.NY-079) awarded to P.W. funded for cell culture and virus purification; Grant from the First Level Young Scholar Program of Jiangsu University, China awarded to P.W. funded for data analysis and writing the manuscript.

Author information




PW and XJG designed the research; HLH, SLZ and QS performed the experiments; HTY, JBZ and QRH prepared the experimental materials; HLH, PW and XJG analyzed the data and wrote the manuscript. All authors have read and approved the final manuscript.

Corresponding authors

Correspondence to Ping Wu or Xijie Guo.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Summary of the sequenced reads. Table S2. Differentially expressed genes of silkworm midguts upon BmNPV infection. Table S3. Differential DNA methylation regions of silkworm midguts upon BmNPV infection. Table S4. Sequences of primers for qRT-PCR and BS-PCR.

Additional file 2: Figure S1.

KEGG pathway enrichment of differentially expressed genes following BmNPV infection.

Additional file 3: Figure S2.

The number of both differentially methylated genes and differentially expressed genes in different regions of genome.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Huang, H., Wu, P., Zhang, S. et al. DNA methylomes and transcriptomes analysis reveal implication of host DNA methylation machinery in BmNPV proliferation in Bombyx mori. BMC Genomics 20, 736 (2019).

Download citation


  • Bombyx mori
  • Epigenetic regulation
  • BmNPV
  • DNA methylation
  • Host-virus interaction
  • Insect immune response, inhibitor of apoptosis