Larvae of the domesticated silkworm strain Jingsong × Haoyue were reared in our laboratory with fresh mulberry leaves at 25°C. BmNPV was propagated in silkworms and the polyhedrons were collected by centrifugation at 8000 rpm. One million polyhedrons from BmNPV were evenly coated on the mulberry leaves with in an area of 4 × 4 cm (length ^∗ width). Fifth instar silkworm larvae (on 3rd day after molt) were fed on these leaves. Larvae at the same instar stage fed on untreated leaves were used as the control group.

Thirty normal healthy silkworms and 30 BmNPV-infected silkworms were collected at 72 h after feeding for midgut dissection. Total RNAs were extracted from normal midguts and BmNPV-infected midguts with TRIzol reagent (Invitrogen, Carlsbad, CA, United States) and mRNAs with polyA were enriched with Oligo-dT. The complete mRNA sequences were disrupted with fragmentation buffer (10 mM ZnCl₂, 10 mM Tris–HCl pH 7.0). The obtained fragments were divided into two parts: one was captured with m6A antibody (ab151230; Abcam, Cambridge, MA, United States) for enrichment of mRNA with m6A modification and the other was used as input for transcriptomic library construction. The fragments collected from the m6A antibody (ab151230; Abcam, Cambridge, MA, United States) were also used in parallel for transcriptomic library construction. The two constructed libraries were used for sequencing with Illumina Hiseq 4000 at Shanghai OE Biotech. Co. Ltd. The raw data were submitted to the Sequence Read Archive with accession numbers (SRR10141250, SRR10141249, SRR10141248, and SRR10141247).

Raw data (raw reads) in fastq format were first processed using Trimmomatic software (Bolger et al., 2014). Clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N, and low-quality reads from the raw data. Then, 250,000 paired reads were randomly extracted from the clean data for alignment with the NT database¹ using blastn software, and those reads with e < 1 × 10^–10 and >80% coverage were collected for further analysis. rRNA reads were removed using SortMeRNA software (Kopylova et al., 2012). After discarding the rRNA reads, the remaining clean reads were mapped to the reference genome (B. mori assembly ASM15162v1) using HISAT2 with default parameters (Kim et al., 2015). Unique reads with high mapping quality were retained, and potential polymerase chain reaction duplicates were marked with Picard².

m6A-seq data quality was evaluated with the Trumpet R package (Zhang et al., 2018). The m6A-enriched peaks in each m6A immunoprecipitation sample were identified by MeTDiffpeak calling software (Cui et al., 2018) (parameters: “FRAGMENT_LENGTH=200, PEAK_CUTOFF_PVALUE = 0.05”), with the corresponding input sample serving as a control. Called peaks were annotated by intersection with gene architecture using ChIPseeker software (Yu et al., 2015). The differential expression of transcripts with m6A methylome between BmNPV-infected and uninfected midguts was detected with MeTDiff with the following parameters: “FRAGMENT_LENGTH=200, PEAK_CUTOFF_PVALUE = 0.05.” The differential peaks were annotated by ChIPseeker software. GO enrichment and KEGG pathway enrichment analysis of peaks and differential peaks were performed using R based on the hypergeometric distribution.

The cellular m6A machinery includes proteins that act as writers, erasers, and readers of m6A (Gokhale et al., 2016). In mammals, RNA m6A methylation is catalyzed by a polyprotein complex composed of methyltransferase-like (METTL) enzymes METTL3, METTL14, and other factors. The cytoplasmic YTH-domain family (YTHDF)1, YTHDF2, and YTHDF3 proteins bind to m6A through their C-terminal YTH domain (Gokhale et al., 2016). METTL3 (NM_019852.4), METTL14 (NM_020961.3), and YTHDF2 (NM_001172828.1) of humans were blasted to the silkworm EST database to predict the writers and readers of the m6A machinery. The domain analysis of writers and readers from the silkworm was conducted with the SMART database³. The molecular weight was calculated with Compute pI/Mw⁴.

The pIZT-V5-his plasmid (Invitrogen) was used to overexpress BmMETTL3 (XM_012695409.2), BmMETTL14 (XM_004924348.3), and BmYTHDF3 (XM_004924548.3) genes in BmN cells. Sequences of the ORFs for the recombinant plasmids pIZT-V5-his-BmMETTL3 (EcoRI/SacII), pIZT-V5-his-BmMETTL14 (KpnI/SacII), and pIZT-V5-his-BmYTHDF3 (EcoRI/XbaI) were inserted into the indicated restriction sites of the pIZT-V5-his vector, respectively. BmN cells (1.5 × 10⁵ cells/well) were seeded in six-well plates and cultured overnight at 28°C. Transfection was performed using Lipofectamine^® 3000 Transfection Reagent (Invitrogen). Briefly, 2 μg recombinant plasmid was mixed with the transfection reagent gently and then added to a 3-cm dish with 2 ml culture medium without 10% fetal calf serum. Four hours later, the medium was replaced with complete medium and cultured at 28°C. Expression levels of BmMETTL3, BmMETTL14, and BmYTHDF3 genes in transfected cells were detected with real-time PCR assay.

According to the mRNA sequences of BmMETTL3, BmMETTL14, and BmYTHDF3, specific siRNA-1/2/3 of BmMETTL3, BmMETTL14, and BmYTHDF3 genes and non-specific siRNA-NC as a negative control were designed and synthesized by Shanghai Integrated Biotech Solutions (Table 1). BmN cells (1.5 × 10⁵ cells/well) were seeded in six-well plates and cultured overnight at 28°C. Transfection was performed using Lipofectamine^® 3000 Transfection Reagent (Invitrogen). Briefly, 2 μg siRNA was mixed with the transfection reagent gently and then added to a 3-cm dish with 2 ml culture medium without 10% fetal calf serum. Four hours later, the medium was replaced with complete medium and cultured at 28°C. The non-specific siRNA-NC (20 nM) was used as the control. BmN cells were subsequently treated using the same procedure as mentioned above for the transient expression experiments. Each transfection was repeated three times. Expression levels of BmMETTL3, BmMETTL14, and BmYTHDF3 genes were detected with real-time PCR to evaluate the extent of gene depletion following siRNA treatment.

BmN cells were cultured for 12 h in 96-well plates before treatment with different siRNAs targeting BmMETTL3, BmMETTL14, BmMETTL3+14, and BmYTHDF3. For a period of 60 h, the treated cells were collected for cell viability detection. Another group of BmN cells was transfected with the recombinant plasmids pIZT-V5-his-BmMETTL3, pIZT-V5-his-BmMETTL14, or pIZT-V5-his-BmYTHDF3. At Forty-eight hours post-transfection, the treated cells were collected for cell viability detection. Cell viability was estimated with the3-(4,5-dimethyl-2-thiazolyl)-2,5-diphenyl-2-H-tetrazolium bromide (MTT) assay. Absorbance was detected at 490 nm with a reference wavelength of 595 nm.

Budded virus containing EGFP-tagged BmNPV (BmNPV-GFP) was kindly gifted by Professor Xiaofeng Wu, Zhejiang University. At 48 h post-transfection, BmNPV-GFP [multiplicity of infection (MOI) = 5] was used to infect siRNA- and plasmid-transfected cells. At 48 h post-infection, the cells were harvested for western blotting with BmNPV structural protein VP39 antibody (gifted by Professor Zhongjian Guo, Jiangsu University) as the primary antibody. α-Tubulin was used as the reference protein, which was detected with anti-TUBB1 mouse monoclonal antibody (Abcam). Goat anti-mouse IgG H&L (Abcam) was used as the secondary antibody. Quantitative analysis of the visible bands from western blotting was performed using the ImageJ program. The integrated density values of the bands corresponding to the detected proteins were normalized to α-tubulin.

The expression levels of BmMETTL3, BmMETTL14, and BmYTHDF3 mRNAs were detected in the cells transfected with expression plasmids or siRNAs with real-time PCR. The procedure for real-time PCR assay was carried out according to the protocol supplied by the manufacturer (Bio-Rad, United States). The expression level of mRNAs was normalized to TIF-4A gene as the internal standard control and calculated using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). All experiments were repeated in triplicate. The primers are listed in Supplementary Table 1.

MeRIP-Seq was used to analyze the m6A modification in the transcriptome of the midguts in 6-day-old fifth instars silkworm larvaes. BmNPV infection was confirmed via observation of typical pathological symptoms, such as swelling of the larvae. In addition, viral structural protein VP39 was detected in the BmNPV-infected midguts but not in uninfected midguts (Figure 1). These findings indicate that B. mori was infected with BmNPV and could be used for further study.

Sequencing data were further analyzed by Hisat2 software (Kim et al., 2015), which was used to align the clean reads to the reference genome of B. mori (assembly ASM15162v1) and acquire genetic information for the m6A modification. A total of 49,934,564 and 62,659,272 reads were identified from uninfected and BmNPV-infected midguts, respectively. After filtering out low-quality data, 36,054,610 (72.20% of the total reads) and 49,054,967 (78.29% of the total reads) high-quality reads were mapped to the reference genome of B. mori (assembly ASM15162v1). Among the mapped reads, 32,838,534 (65.76%) and 43,853,738 (69.99%) were uniquely mapped to the genome and 13,188,800 (26.41%) and 17,221,827 (27.48%) were mapped to splice reads, respectively (Table 2).

To validate the efficiency of the RIP-seq data, R-package-Trumpet was used for quality control and to confirm the reliability of the experimental results (Zhang et al., 2018). The quality of m6A-seq data was shown to be sufficient. The enrichment analysis of m6A peaks on mRNAs showed that most peaks were distributed on the 3′-end of mRNAs (Figures 2A,B). Analysis of the sequencing data of immunoprecipitates with m6A antibody identified 7384 and 9144 peaks in samples from uninfected midguts and BmNPV-infected midguts, respectively (Figure 2C). The average length of the peaks was 869.77 and 656.94 bp in uninfected and BmNPV-infected midguts, respectively. The results suggest that the transcriptome of BmNPV-infected midgut, which had more A sites compared to the uninfected samples, had been modified by the related methyltransferase enzymes of the m6A machinery.

The length and the enrichment multiple distribution of m6A methylation peaks were analyzed using the RNA sequences with m6A sites. The length of mRNA with m6A sites from BmNPV-infected midgut was similar to that of uninfected midgut (Figure 2D), while the enrichment multiple distribution analysis showed that more reads were enriched in the peaks of the BmNPV-infected midgut than in those of the uninfected midgut (Figure 2E). These results also demonstrated that more transcripts in the BmNPV-infected midgut may have been modified with the related methyltransferase enzymes from the m6A machinery.

To understand the preferential location of the m6A modification in the transcripts, we investigated the localization of the m6A peaks in the gene functional elements of the entire transcriptome from both uninfected and BmNPV-infected midguts. The distribution of m6A peaks on gene functional elements showed that most of the peaks were distributed at the 3′-UTR and other exons in both BmNPV-infected and uninfected midguts (Figures 3A,B). Considering that the annotation results of gene elements corresponding to the peaks may overlap, we further analyzed the distribution of detected peaks on gene functional elements. The annotation results showed that most of the peaks appeared in the exons, introns, 5′-UTRs and 3′-UTRs (Figures 3C,D). These results indicated that the m6A modification is closely related to RNA expression and may have hotspot regions in the host transcripts.

To explore the potential signaling pathways related to m6A-containing transcripts, GO analysis was carried out by Tritrydb online. A total of 305 IDs were uploaded to Tritrydb and analyzed with the GO Enrichment tool. Default settings were used; however, the results were limited to GO slim terms. The pathways related to m6A-genes were enriched in GO terms referring to binding (209), catalytic activity (109), transporter activity (14), protein binding transcription factor activity (11), molecular transducer activity (10), and structural molecule activity (9) in the molecular function category. Many genes were enriched in cellular process (239), metabolic process (176), biological regulation (157), regulation of biological process (149), cellular component organization or biogenesis (126), and response to stimulus (104) in biological process (Figures 4A,B). The pathways related to the mRNAs modified by m6A were enriched in spliceosome, ribosome, and RNA transport in BmNPV-infected midguts, while proteasome, RNA transport, and spliceosome were enriched in uninfected midguts (Figures 4C,D). These results indicate that the m6A modification may be closely related to gene expression and regulation.

To uncover the differential m6A methylation between uninfected and BmNPV-infected midgut, statistics were applied with threshold p < 0.05. We identified 1221 significantly differentially expressed m6A transcripts, of which 788 and 433 were upregulated and downregulated, respectively, during BmNPV infection. The distribution of the differentially expressed peaks on transcripts was nearly identical between the two groups (Figures 5A,B). Based on GO analysis of these differentially expressed m6A transcripts, the major molecular functions in these transcripts were related to binding activity, catalytic activity, structural molecular activity, molecular transducer activity, transporter activity, and enzyme regulator activity (Figure 5C). Based on KEGG pathway analysis transcripts with differential m6A modification associated with signal transduction and translation (Figure 5D). These results indicate that m6A methylation may play a crucial role in numerous RNA biology events and that the appearance of peaks of m6A in transcripts may be related to BmNPV infection.

Proteins that act as writers and readers of m6A were identified in the silkworm genome and subjected to bioinformatics analysis. Human enzymes METTL3 (NM_019852.4) and METTL14 (NM_020961.3) were blasted to the silkworm EST database, which resulted in the identification of B. mori N6-adenosine-methyltransferase 70 kDa subunit (XM_012695409.2) and B. mori METTL protein 14 homolog (XM_004924348.3) as candidate writers of the silkworm m6A machinery, named as BmMETTL3 and BmMETTL14, respectively. The human cytoplasmic YTHDF2 (NM_001172828.1) was also blasted to the silkworm EST database, and B. mori YTH domain-containing family protein 3 (XM_004924548.3) was predicted to be a reader of m6A machinery, named as BmYTHDF3.

The in silico analysis showed that the coding sequences of BmMETTL3, BmMETTL14, and BmYTHDF3 were 1710, 1146, and 1959 bp corresponding to molecular weights of 62.98, 40.09, and 73.95 kDa, respectively. The specific domains in the m6A machinery components were predicted with the SMART database, which showed that BmMETTL3 and BmMETTL14 contained an MT-A70 domain, originally identified as the S-adenosylmethionine-binding subunit of human N6-adenosine-methyltransferase (MTase), an enzyme that sequence-specifically methylates adenines in pre-mRNAs. The YTH domain in BmYTHDF3 corresponds to an evolutionarily conserved m6A-dependent RNA binding domain (Figure 6). These domains from the three proteins in the m6A machinery were similar to those identified in other species, which confirmed that the m6A machinery is conserved across different species.

To detect the effect of the m6A modification of mRNAs on BmNPV infection, BmMETTL3 and BmMETTL14 and m6A reader protein BmYTHDF3 were depleted by siRNA in BmN cells and these cells were subsequently infected with BmNPV (MOI = 5). To understand whether the depletion of the m6A machinery genes had an adverse effect on cells, BmN cell viability was examined with the MTT assay and the results showed that depletion of BmMETTL3, BmMETTL14, BmMETTL3+14 (methyltransferase complex), and BmYTHDF3 by >50% did not impair the cells (Figures 7A,B). At 48 h post-infection, the infected cells were harvested and subjected to western blotting using VP39 antibody. Depletion of BmMETTL3, BmMETTL14, BmMETTL3+14, and BmYTHDF3 significantly increased expression of viral protein VP39, compared with cells treated with non-targeted control siRNA (Figure 7C).

To detect the effect of m6A on BmNPV infection, BmMETTL3 and BmMETTL14 and m6A reader protein BmYTHDF3 were overexpressed by recombinant plasmids in BmN cells and these cells were infected with BmNPV (MOI = 5). The results showed that the expression levels of BmMETTL3, BmMETTL14, and BmYTHDF3 were significantly increased by recombinant pIZT-V5 vectors following transfection (Figure 8A), and that overexpression of these genes in BmN cells significantly decreased expression of VP39 (Figure 8B). These results demonstrated that the m6A machinery may play an important role in the process of BmNPV infection.