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℃ 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 12h light and 12h 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 12h. 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 (http://www.genome.jp/kegg/)[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-300bp, followed by end repair, adenylation and ligating methylated adaptors. Then these DNA fragments were treated twice with bisulfite using EZ DNA Methylation-GoldTM 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) [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 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:
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 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 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 (http://quma.cdb.riken.jp/).