Identification and Characterization of BmNPV m6A Sites and Their Possible Roles During Viral Infection

Bombyx mori nucleopolyhedrovirus (BmNPV) is one of the most serious pathogens and causes serious economic losses in sericulture. At present, there is no epigenetic modification of BmNPV transcripts, especially of m6A, and this modification mediates diverse cellular and viral functions. This study showed that m6A modifications are widespread in BmNPV transcripts in virally infected cells and the identified m6A peaks with a conserved RRACH sequence. m6A sites predominantly appear in the coding sequences (CDS) and the 3′-end of CDS. About 37% of viral genes with m6A sites deleted from the viral genome did not produce any infectious virions in KOV-transfected cells. Among the viral genes related to replication and proliferation, ie-1 mRNA was identified with a higher m6A level than other viral genes. The m6A sites in the ie-1 mRNA may be negatively related to the protein expression. Viral replication was markedly inhibited in cells overexpressed with BmYTHDF3 in a dose-dependent manner, and a contrary effect was found in si-BmYTHDF3-transfected cells. Collectively, the identification of putative m6A modification in BmNPV transcripts provides a foundation for comprehensively understanding the viral infection, replication, and pathobiology in silkworms.


INTRODUCTION
Bombyx mori nucleopolyhedrovirus (BmNPV) is one of the most serious pathogens and causes serious economic losses in sericulture (1). Its genome is circular double-strand DNA with 130 kb and encodes 141 putative proteins (2). BmNPV is a classical member of baculoviruses with two types of virions in the whole life cycle. The budded virus (BV) is responsible for the horizontal transmission among the cells, and the occlusion-derived virus (ODV) is responsible for the vertical transmission (3). BmNPV belongs to the nuclear-replicated virus and is inevitably modified by a series of host modification enzymes, including m6A modification. N 6 -Methyladenosine (m6A) is one of eukaryotic mRNA's most abundant internal modifications. It is also a dynamic reversible RNA modification, which plays an important role in the structure, positioning, and function of RNA (4,5). m6A is involved in various biological processes, including stress response, fertility, stem cell differentiation, circadian rhythm, microRNA biogenesis, and m6A modification has a positive regulatory effect in the life cycle of Simian virus 40 (SV40), and overexpression of YTHDF2 can accelerate virus replication (30). The m6A modification site located in the 3′UTRs and 3′ epsilon loop of the hepatitis B virus (HBV) genome is not conducive to the stability of HBV RNA. The m6A modification site located in the 5′ ϵ ring of the HBV genome has a positive regulatory effect on pregenomic RNA (pgRNA) reverse transcription (31). EBV virus (EBV) transcripts are many m6A modification sites. YTHDF1 can recognize m6A sites accelerating viral RNA uncapping and recruit RNA degradation complexes to inhibit viral infection and replication (32). m6A modification in adenovirus (ADV) transcripts involves regulating splicing and viral RNA expression (33). The porcine epidemic diarrhea virus (PEDV) genome contains multiple m6A modification sites and negatively regulates the virus infection (34). Our preliminary findings showed that the expression levels of writers (BmMETTL3 and BmMETTL3) and reader (BmYTHDF3) of m6A machinery in BmN cells were negative to the expression of viral structural protein VP39 in BmNPV infection (35). These studies indicate that both DNA and RNA viruses are regulated by the m6A modification system in host cells during viral infection, but the m6A modification system has different regulatory mechanisms for different viruses.
BmNPV belongs to the nuclear-replicated virus. Whether m6A modification is introduced into the viral transcripts in the viral proliferation process, and whether this apparent modification is beneficial for the virus to evade the detection of the host's innate immune system, is still unclear. In this study, the m6A modification was investigated by methylated RNA immunoprecipitation (MeRIP) sequencing from the BmNPVinfected midgut, and m6A sites were widespread in the viral transcripts. Furthermore, the m6A machinery negatively mediated the viral replication. The results will uncover the roles of BmNPV transcript m6A modification in the viral life cycle and provide a scientific base for anti-BmNPV.

Sample Preparation
The domesticated silkworm strain (Jingsong×Haoyue) was reared in our lab with fresh mulberry leaves at 25°C. The polyhedrin (1*10 6 ) of BmNPV stored in our lab at Soochow University was evenly coated on the mulberry leaves with 4*4 cm for the first day of 5 instar larvae. 48 h postinfection, the midguts were collected and extracted for total RNA. BmN cells were cultured in TC-100 medium with 10% fetal bovine serum (HyClone, Logan, UT, USA) in a 26°C incubator. The recombinant BmNPV-GFP was gifted by Professor Xiaofeng Wu from the Zhejiang University.
Methylated RNA Immunoprecipitation Sequencing 30 silkworm larvae were infected by BmNPV, and the midguts were used to extract the total RNA with the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The total RNA was fragmented with the corresponding buffer (10 mM ZnCl 2 , 10 mM Tris-HCl pH 7.0). The obtained fragments were captured with the m6A antibody, and the corresponding libraries were utilized for MeRIP sequencing (MeRIP-Seq), and further bioinformatics analysis was also conducted by Shanghai OE Biotech. Co., Ltd. (Shanghai, China). The sequencing data were deposited in the NCBI with accession numbers (SRR10141250, SRR10141249, SRR10141248, and SRR10141247) (35). BmNPV T3 (accession number: L33180.1) was used as the reference genome for mapping.

MeRIP-qPCR
The method for preparation of total RNA for MeRIP-qPCR was the same protocol as MeRIP sequencing, except that total RNA was not treated with the fragmentation buffer. Briefly, BmNPVinfected BmN cells were lysed with lysis buffer. The lysed cells were separated into parts. One part of lysate was incubated with an anti-m6A antibody, and another part was incubated with IgG control (Rabbit) conjugated to Protein A/G beads in 1 ml RNA immunoprecipitation (RIP) buffer. After incubation at 4°C overnight, the beads were washed with wash buffer 3 times, and the RIP samples were reverse transcribed into cDNA using the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China) and subjected to RT-qPCR. Translation initiation factor 4a (TIF-4A) was used as a housekeeping gene (36). The specific primers for each viral gene used in real-time PCR are listed in Supplementary Table 1. All experiments were conducted three times.

Plasmids and siRNA
The genome loci of ie-1 in BmNPV T3 (L33180.1) are 116994.118748 and were synthesized by the Sangon Biotechnology company (Shanghai, China). The DNA sequence was inserted into the Kpn I and Xba I sites in the pIZT-V5/his vector (Invitrogen, Frederick, MD, USA) for pIZT-V5his-ie. A mutant fragment of ie-1 was mutated A to T at the potential m6A sites. This mutant DNA sequence was also synthesized by the Sangon Biotechnology company (Shanghai, China) and was inserted into the sites of Kpn I and Xba I in the pIZT-V5/his vector for pIZT-V5his-ie-1mut. pIZT-BmYTHDF3 and siRNA (si-BmYTHDF3) specifically targeting BmYTHDF3 were reported in our previous study (35).

Overexpression and Depletion of BmYTHDF3 in BmN Cells and BmNPV Infection
BmN cells (1*10 5 cells/well) were transfected with a recombinant expression plasmid (pIZT-BmYTHDF3) with the transfection reagent. 48 h post-transfection, BmN cells were infected with BmNPV-GFP (MOI = 2) and harvested for extraction of total proteins at 12 h postinfection for Western blotting. Meanwhile, BmN cells (1*10 5 cells/well) were transfected with specific siRNA targeting BmYTHDF3 (si-BmYTHDF3) with the transfection reagent. 48 h post-transfection, BmN cells were infected with BmNPV-GFP (MOI = 2) and harvested for extraction of total proteins at 12 h postinfection for Western blotting. The expression of the Bm59 gene of BmNPV was used as the indicator of viral replication (37).
Bioinformatics m6A modification sites on the interested RNA sequences were carried out by the SRAMP prediction server (38). The motif discovery from the sets of sequences was investigated by HOMER software (39). The flexibility of protein local structures was studied through (i) the B-factor of the X-ray experiment and (ii) the fluctuation of residues during molecular dynamics simulations (40). The Flexibility Prediction was utilized to analyze the flexibility classes (0, 1, and 2 represent rigid, intermediate, and flexible, respectively) according to the observed dynamic properties from the target protein sequences.

m6A Modifications Are Widespread in BmNPV Transcripts in Virally Infected Cells
BmNPV-infected silkworm midguts were collected and sequenced with methylated RNA immunoprecipitation combined with high-throughput sequencing (MeRIP-Seq) to understand the m6A modification in mRNAs, including host mRNAs and viral mRNAs ( Figure 1A). Our previous study found that thousands of m6A peaks were identified from BmNPV-infected tissues, and the host m6A modification system was related to the expression of the viral gene. However, this modification is still a puzzle whether it appears in the BmNPV mRNAs. To better understand the posttranscription modification, especially m6A sites on the viral genes and the potential roles in the viral life cycle, we compared the MeRIP-Seq data with the BmNPV T3 (L33180.1) reference genome and found that 74 regions with high potential m6A methylation (m6A peak) were identified in 59 mRNA transcripts of BmNPV ( Figure 1B and Table 1). Among these m6A peaks, the mean size of the m6A peak is 262 bp, while the maximum and minimum peaks are 599 and 135 bp, respectively. HOMER motif enrichment analysis was utilized to analyze the RRACH motif conservation in 74 identified m6A peaks. The results showed that four conserved RRACH motifs were identified from the alignments between the de novo motif and known motifs and met the requirement of the motif domain of RRAmCH (R=G/A, G > A; H=A, C, U) ( Figure 1C). The positions of m6A sites were analyzed to understand the appearance of m6A sites in viral mRNAs. The results indicated that m6A modification mainly appeared in coding sequences (CDS) and the 3′-end of CDS ( Figure 1D). These results suggested a potential biofunctional significance for m6A modification in viral mRNA.

Viral Genes With m6A Sites Were Classified Into Four Phenotypes
The KOV (bacmid)s were subdivided into four phenotypes (A to D) according to the results from the knockout BmNPVs (KOVs) for each viral gene using the lambda red recombination system (41). The defined phenotypes of type-A and -B KOVs could produce infectious viruses, but types C and D did not produce any infectious virions in KOV-transfected cells (2). According to the appearance of m6A in the viral genes, we found 33, 4, 17, and 5 viral genes in types A, B, C, and D, respectively ( Table 2).

Validation of m6A Modification in Viral Genes
To know the existence of the m6A site in viral genes, multiple genes related to viral replication (lef-8, lef-5, lef-3, VP80, GP64/ 78, and ie-1) were selected for validation using m6A antibody immunoprecipitation combined with qPCR (MeRIP-qPCR). The results showed that all these viral genes could be enriched by the m6A antibody and detected by qPCR (real-time PCR) using the corresponding primers of each gene (Supplementary Table 1). Compared to the relative expression levels of viral genes in the complex of the m6A antibody and viral mRNAs with m6A modification, the results showed that ie-1 mRNA with a higher m6A level than other viral genes was identified from the complex ( Figure 2), indicating that enriched m6A modification in the ie-1 gene may be related to the transcription or translation.

The Expression Level of IE-1 Was Related to the m6A Sites
To further validate the function of m6A sites in the ie-1 mRNA, the wild-type and mutant ie-1 plasmids were constructed by cloning a wild-type 1,752-nucleotide (nt) length fragment (remove the TAA) encompassing the m6A peak region in the ie-1 mRNA and a mutant fragment with A > T mutations at the potential m6A sites, respectively ( Figure 3A). The plasmids were transfected into BmN cells for Western blotting assay, which exhibited a higher IE-1 expression level on the mutant plasmid than wild-type ie-1 ( Figure 3B), confirming that the presence of m6A sites in the ie-1 mRNA may be related to its expression. To exclude the influence of codon on protein expression, the flexibility of protein local structures was studied through (i) the B-factor of the X-ray experiment and (ii) the fluctuation of residues during molecular dynamics simulations (40). The prediction results showed that these 5 amino acid mutants in the IE-1 had not altered the protein flexibility (Figure 4).

m6A Reader (BmYTHDF3) Negatively Regulated the Viral Replication
Our previous study concluded that cytoplasmic YTH-domain family 3 (BmYTHDF3) was indicated as the m6A reader of the m6A machinery in silkworm, Bombyx mori, and the m6A modification system was predicted to play important roles in the BmNPV infection (35). The appearance of m6A modification in the ie-1 gene means that BmYTHDF3 could recognize these m6A sites. However, the exact roles of BmYTHDF3 in the replication of BmNPV are still not clear, and we used the recombinant expression plasmid (pIZT-BmYTHDF3) and

DISCUSSION
Increasing data indicated that m6A modification not only exists in the cellular mRNAs but also is widespread in the viral transcripts, including DNA and RNA viruses. Viral m6A modifications have been reported as pro-viral or antiviral functions depending on the virus species, cell type, and viral RNA location (20). We showed that m6A modification is widespread in BmNPV transcripts in virally infected cells. Viral genes with m6A sites were classified into four phenotypes according to the results from the knockout BmNPVs (KOVs) for each viral gene using the lambda red recombination system (41). We demonstrated that ie-1 mRNA with a higher m6A level than other viral genes was examined, and the expression level of IE-1 was related to the m6A sites in the viral ie-1 gene. Furthermore, m6A reader protein BmYTHDF3 suppressed the BmNPV infection in BmN cells, maybe recognizing the m6A sites in the viral gene. m6A machinery regulates cellular RNAs of viral RNAs by writers, erasers, and readers in cells (42). The appearance of m6A modification in viral RNAs and cellular m6A modification systems played vital roles in the viral life cycle. Our previous study found that the cellular mRNAs with m6A modification were altered following BmNPV infection and the expression levels of writers and readers negatively mediated BmNPV infection (35). However, the precise regulation of m6A modification on viral infection was vague. Therefore, we made assumptions about widespread m6A sites in the BmNPV transcripts like other viruses, and the viral RNA with m6A sites could be regulated by cellular m6A reader proteins. Presently, any reports showed that BmNPV transcripts could be modified by the m6A machinery. MeRIP-Seq data were aligned with the BmNPV T3 (L33180.1) reference genome, and 74 regions with high potential m6A methylation (m6A peak) were distributed in 59 viral mRNA transcripts. Among these A and -B KOVs could produce infectious viruses, but type-C and -D did not produce any infectious virions in KOV-transfected cells. m6A peaks, four conserved RRACH motifs were identified as other species. The enrichment analysis of m6A peaks on viral mRNAs showed that most were distributed on the CDS and 3′-CDS of viral mRNAs. It is known that m6A modification was found near the 3′ UTR and stop codon in mammals and mice (43). The methylation landscape across human adult tissues showed m6A sites preferentially around stop codons, essential for maintaining basic cellular function. Most of the m6A in the genes with tissue specificity were found in the 5′-UTR (44), suggesting that these viral genes with m6A sites may be essential for distinct regulatory effects in the aspect of host range, viral replication, assembly, and host-virus interaction. 22/59 viral genes were demonstrated to be required to produce infectious viruses, and the m6A site appearance in these genes may be related to the viral replication. Furthermore, several viral genes related to viral replication were selected to be validated by MeRIP-qPCR, and the results also showed that m6A modifications were enriched in these selected viral genes. Among them, the modification level of m6A on the ie-1 gene was higher than other genes, indicating that these m6A sites might regulate their transcription or translation.
In our previous study, we found that the expression level of VP39 was decreased in BmN cells overexpressed with BmMETTL3 or BmYTHDF3, while the reverse results were obtained in the BmN cell depletion of BmMETTL3 or BmYTHDF3 using corresponding specific siRNAs (35). These suggested that the m6A modification system plays an important role in BmNPV infection with an unknown mechanism. The role of the YTHDF readers on viral replication was extensively examined by depletion and overexpression studies in numerous DNA and RNA viruses (29,32,45). Increasing data showed that YTHDF readers regulate the mRNA fate of cellular or viral m6A RNAs. It is described to regulate splicing, nuclear export, cap-independent translation, and decay (46). To better understand the regulatory effects of m6A machinery on the expression of viral genes or viral replication, the expression level of IE-1 was detected from BmN cells transfected with the wild-type and mutant ie-1 plasmids, and Western blotting assay exhibited a higher IE-1 expression level on the mutant plasmid than wild-type ie-1, confirming that the presence of m6A sites in the ie-1 mRNA may be related to its expression. When these mutant sites were introduced into ie-1, any adverse effect on protein flexibility was detected. Using the BmNPV-GFP-infected BmN cells overexpressed with BmYTHDF3 and transfected with si-YTHDF3, we found that the EGFP and viral gene Bm59 were significantly decreased and increased, respectively. These results suggest that BmYTHDF3 suppressed the BmNPV infection in BmN cells via recognizing the m6A sites in viral genes. The m6A modification system has different regulatory mechanisms for different viruses (45,47). The appearance of m6A modification on viral RNA will be used as a molecular marker for evading host innate immune recognition of non-self RNA (28). Altogether, many BmNPV transcripts have been modified by the m6A machinery, including many genes related to viral infectivity. BmYTHDF3 suppressed the BmNPV infection by recognizing the m6A sites in viral genes, resulting in impairment of viral gene expression, avoiding viral overproliferation.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
XH and CG contributed to the conception and design of the study. XZ, YZ, and JP performed the detection and analysis.
XZ and XH wrote the first draft of the manuscript. CG wrote sections of the manuscript. All authors contributed to the manuscript revision and read and approved the submitted version.