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

Background 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. Results 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. Conclusion 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.

Background 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 viruse, 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 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 M. anisopliae caused the differential expression of Dnmt genes in the larvae of G. mellonella infected in a natural manner, suggesting that DNA methylation may play roles in the immune response of insects against parasitic fungi [17].Our previous study found that 27 genes were shown to have both differential expression and differential methylation in the midgut and fat body of B. 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 [23].
BmNPV, a circular double-stranded DNA (dsDNA) virus, belongs to the family Baculoviridae[24] and is a big challenge for the development and sustainability of the sericultural industry [25]. 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 [26]. Up to now, there has been scarce report on epigenetic function on BmNPV infection and immune response in B. mori.
In this study, we treated the B. 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 12h followed by MTT assay. The result showed that there were no significant differences in cell survival between 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.
In order to examine whether DNA methyltransferase inhibition affects BmNPV replication, BmN cells were then treated with 100 μM Zeb for 12h and then replaced with a new medium followed by infection of recombinant BmNPV BVs containing an egfp marker gene for 72h. By observing under the fluorescence microscope we observed that the fluorescence intensity of the treated group (Zeb/BmNPV) was significantly weaker 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 patterns of these two genes were very similar with the tendency of increase at 6h while decrease at 12h-48h in Zeb-treatedcells (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 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 (http://www.ncbi.nlm.nih.gov/sra/) under the accession number PRJNA541379. The number of significantly differentially expressed genes (DEGs) in infected midgut was 2171, with 920 genes up-regulated and 1251 genes down-regulated ( 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: Fig. S1).

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.29% of genomic cytosine is methylated including mCs at CG, CHG and CHH sites, and the average methylation level at CG location is about 0.95%. By using Circos software, which visualizes data in circular layout[27], the global difference methylation level between T group and C group was displayed in Fig. 3A. Relatively, the methylation level of scaffold 16 is highest than that of other scaffolds. On scaffold 14, 11 and 6, some methylated sites were shown hyper-methylation while on scaffold 15, 19 and 3, a lot of hypomethylated sites were observed. Furthermore, we comparatively analyzed the methylation level on different genomic functional regions between two groups and the results revealed that in silkworm DNA methylation is mostly targeted to gene bodies. It peaked in exon region and sharply decreased in intron region. Promoter region shows more methylation than the repeat region (Fig. 3B). Fractional methylation levels in exon region exhibited the tendency of hyper-methylation in control midguts than in BmNPV-infected midguts (Fig. 3B).

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).

Effects of DMRs on gene expression in silkworm challenged by BmNPV infection
To investigate whether the variances in CG methylation detected in silkworm following BmNPV infection changed 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: Fig. 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 B. mori [28]. Furthermore, by qRT-PCR and BS-PCR, we validated the up-regulation (Fig.5A) and hypermethylation of BGIBMGA014008 (Fig. 5B, 5C), encoding for structural maintenance of chromosomes protein in BmNPV-infected midgut.

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 changes 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 whilethat of stat has no significant change ( 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 [32,33]. In this study, we found 26 DEGs showing differential methylated levels in the midgut of BmNPV-infected silkworm larvae with 17 hypermethylated and 9 hypomethylated, 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 [34]. In addition, SMC4 expression was significantly associated with tumor size, dedifferentiation, advanced stages and vascular invasion of the primary liver cancers [35]. 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 [36]. It is reported that down-regulation of BmDhodh inhibits cell growth and the endomitotic DNA replication in silk gland cells [37]. In our study, we observed the increased mRNA level and hypomethylation in promoter region of BmDhodh (BGIBMGA011887). Generally, promoter methylation represses gene expression [38]. 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 that replicate in the nucleus. Regulation of this process is paramount to processes such as cell division and differentiation, but is also critically important for viral replication and pathogenesis [39].
Nucleocytoplasmic transport is mediated by nuclear pore complexes (NPCs), which are embedded in the nuclear envelope [40]. 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 [41,42]. 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, which could affect the replication of BmNPV.
Other interesting genes, such as Rad50 (BGIBMGA005449), which is associated with DNA repair functions [43][44][45] and the histidine triad protein gene (BGIBMGA011899), which plays roles in transcription, signal transduction and many peripheral and central nervous system diseases [46] 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 [47]. 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 [48,49]. Studies on virus-driven dysregulation of host immune-related gene expression through DNA methylation present a novel viral mechanism to inhibit immune responses [50]. Many DNA and RNA virus could utilize this mechanism to downregulates expression of host immune-related genes [6,51,52].
Previous reports revealed that BmNPV infection could activate the JAK/STAT signaling pathway and has slight effects on Imd and Toll signaling pathways [53][54][55]. The JAK signal transducer and STAT signaling pathway is an important regulator of cell proliferation, differentiation, survival, apoptosis and immune response [56]. 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 [57], displayed the opposite expression pattern with up-regulation of socs6 and down-regulation of socs2 and stat expression has no significant change (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 [58]. Inhibitor of apoptosis (IAP) protein family have demonstrated functions in both apoptosis and innate immunity [59]. Both BmIAP and Baculoviruses encode IAP had an inhibitory effect on apoptosis in insects [60,61]. Survivin, another apoptosis inhibitor has been confirmed to be an essential regulator of cell division [62]. In our study, we found the 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.

Conclusions
In summary, B. mori, one of the insect species with completely sequenced genomes is well proven to be a good model organism. The application of this model extends the molecular mechanism of host-

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 [63,64]. As to BmNPV gene lef-3 (GeneID: 1488686) and ie-1 (GeneID: 1488755), absolute qRT-PCR was carried out as in our previous report [63]. 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 is 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 [63]. 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). Then, the membrane was further incubated with HRP labeled goat anti-rabbit IgG (1:20000; Beyotime, China).
The same blots were re-probed with ɑ-Tubulin antibody (1:5000; Sigma, USA) to confirm equal loading of samples.
Transcriptome analysis RNA-Seq experiment 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 (http://www.genome.jp/kegg/) [65].

Whole genome bisulfite sequencing
For whole-genome bisulfite sequencing,a total amount of 5.2 mg genomic DNA were fragmented by sonication to 200-300bp, followed by end repair, adenylation and ligating methylated adaptors. Then these DNA fragments were treated twice with bisulfite using EZ DNA Methylation-Gold TM Kit (Zymo Research) before PCR amplification. Finally, sequencing was performed on an Illumina Hiseq 4000 platform and 125bp/150bp 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) [66] 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 [67]. 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 10kb. 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:

Differentially methylated analysis
Differentially methylated regions (DMRs) were identified using the DSS software(version DSS_2.12.0) [8,64,68]. 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 100bp, 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 ZymoTaq TN 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.        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.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.