Skip to main content

The dynamics of N6-methyladenine RNA modification in interactions between rice and plant viruses

Abstract

Background

N6-methyladenosine (m6A) is the most common RNA modification in eukaryotes and has been implicated as a novel epigenetic marker that is involved in various biological processes. The pattern and functional dissection of m6A in the regulation of several major human viral diseases have already been reported. However, the patterns and functions of m6A distribution in plant disease bursting remain largely unknown.

Results

We analyse the high-quality m6A methylomes in rice plants infected with two devastating viruses. We find that the m6A methylation is mainly associated with genes that are not actively expressed in virus-infected rice plants. We also detect different m6A peak distributions on the same gene, which may contribute to different antiviral modes between rice stripe virus or rice black-stripe dwarf virus infection. Interestingly, we observe increased levels of m6A methylation in rice plant response to virus infection. Several antiviral pathway-related genes, such as RNA silencing-, resistance-, and fundamental antiviral phytohormone metabolic-related genes, are also m6A methylated. The level of m6A methylation is tightly associated with its relative expression levels.

Conclusions

We revealed the dynamics of m6A modification during the interaction between rice and viruses, which may act as a main regulatory strategy in gene expression. Our investigations highlight the significance of m6A modifications in interactions between plant and viruses, especially in regulating the expression of genes involved in key pathways.

Background

N6-methyladenine (m6A) RNA methylation is one of the most common RNA modifications in prokaryotes and eukaryotes [1]. m6A methylation and its biological functions in prokaryotes and eukaryotes have been the focus of research, recently [2,3,4,5,6]. Modification of m6A methylation was first reported 40 years ago [7, 8]. The modification contributes to the generation, localisation, and functionality of messenger RNA (mRNA) through the regulation of stability and translation [2,3,4,5]. Our understanding of the essential roles of m6A in viruses, fungi, animals, and plants is increasing [9,10,11]. Most studies have focused on the influence of m6A on development, evolution, and physiology, especially in plants [12,13,14,15,16]. However, little is known about the precise functions of m6A in the interactions between plants and parasites, and whether these functions are involved in physiological and pathological changes.

Methylation of adenosine is catalysed by “WRITER”, a large molecular weight (> 1 MDa) RNA methyltransferase complex in plants. It is composed of two MTA-70 family proteins (MTA and MTB), FK506-binding protein 12 (FKBP12) interacting protein 37 (FIP37), VIRILIZER (VIR), and ubiquitin ligase HAKAI [17,18,19]. The m6A methylation is reversible and can be removed by “ERASER”, which includes two main components in Arabidopsis: ALKBH9B and ALKBH10B. Recognition of the m6A methylated genes involves “READER”, which recognises YTH domain-containing proteins, such as evolutionarily conserved C-terminal regions (ETC2), ETC3, and ETC4 [12, 20]. m6A methylation is performed by adding the methyl group to the N6 position of adenosine. S-adenosylmethionine (SAM) often serves as a methyl “DONOR” for almost all cellular methylation reactions [21, 22]. SAM generation from methionine and adenosine triphosphate (ATP) involves SAM synthetases [22]. The “WRITER”, “READER”, “ERASER”, and “DONOR” proteins are tightly associated with multiple biological processes in plants. Previous studies showed that loss-of-function mutants of “WRITER” genes (FIP37 and OsFIP) resulted in serious developmental malformations, and even embryonic lethality [14, 17]. The “ERASER” genes, which include AtALKBH9B and AtALKBH10B, affect the infectivity of the alfalfa mosaic virus (AMV) to Arabidopsis because of the interaction of the viral capsid protein (CP) and the eraser protein; these genes also affect the floral transition of Arabidopsis due to the altered stability of mRNAs targeting key flowering time genes [20, 23]. “READER” genes are involved in leaf formation and trichome morphology. The binding of ECT2 to the RNA “UGUA” m6A motif could regulate the transcript stability of trichome morphogenesis-related genes [12, 24]. The collective findings concerning m6A biological functions strongly indicate that m6A methylation has important roles in the regulation of gene expression in plants, such as the regulation of the mRNA stability of target genes.

Although the precise molecular functions of m6A dynamics are not fully understood, several studies have revealed relationships with RNA stability [25], triggering RNA structure switches [26], translation [27], splicing [28], microRNA (miRNA) processing [29], and RNA export [30]. The recent development of the high-throughput methylated RNA immunoprecipitation sequencing (MeRIP-seq) technology enabled the identification of transcriptome-wide m6A modification. The findings demonstrated the enrichment of m6A around the start codon, stop codon, and 3′-untranslated regions (3′-UTRs) with an “RRACH” consensus motif in all eukaryotes that have been analysed [15, 31, 32]. These findings strongly suggest a conserved mechanism of m6A deposition in eukaryotic mRNA and stimulated many hypotheses about their roles.

The dynamics of m6A modifications include the exact sites, frequency of methylation, and percentage of methylated genes. These dynamics may vary in plants exposed to biotic and abiotic stress, especially under parasitic infection. m6A modification is important in the regulation of viral replication and the viral life cycle in animal systems [6, 11, 33,34,35]. However, in plants, most viruses are RNA viruses. m6A modification, as a widespread modification in plants, may have profound potential roles in modulating virus infection. In Arabidopsis, the genomic RNA accumulation of AMV was reportedly decreased in the T-DNA insertion mutant of Atalkbh9b, and the virus infectivity was impaired, whereas the infectivity of cucumber mosaic virus (CMV) was not altered [23]. This may be because AtALKBH9B can interact with the CP of AMV, but not with the CP of CMV [23]. Infection of tobacco plants with tobacco mosaic virus (TMV) significantly reduced the overall m6A modification levels [36]. In addition, a conserved domain-containing ALKB has been identified in the genomic RNA of several single-stranded plant RNA viruses [37, 38]. These results imply that the m6A modification may act as a fine modulation mechanism for plants responding to viral infection. Some plant viruses have evolved strategies to defend against the host m6A modulation system. However, the m6A dynamics in interactions between plant and virus remain unclear, as do the molecular functions of the modification and the relationship between the expression levels of host main disease resistance pathway-related genes and the m6A modification region on the gene bodies.

Presently, combined MeRIP-seq and transcriptome analyses revealed the activation of the overall m6A modification levels during plant virus infection. Further, the distribution of m6A peaks in both viral and rice genomes were mapped for the first time. The m6A modification was tightly associated with genes that were not actively expressed in rice infected with viruses. Gene ontology (GO) analyses showed that RNA binding activity apparently influenced the molecular functions. Moreover, the most common consensus was analysed in rice with and without viruses’ infection.

We also found m6A levels were associated with the expression of the key genes involved in jasmonate (JA)-mediated RNA silencing. Our findings also revealed the involvement of the m6A modification in the relative expression of the main antiviral pathway-related genes in plants, such as the genes of main m6A methylation machinery, RNA silencing, and phytohormone metabolism. These data provide evidence that m6A modification participate in and alter the physiological and pathological status of rice plants during interactions with viruses.

Results

Confirmation of infection with RSV and RBSDV in rice seedlings

To confirm the infection of rice plants with the two plant viruses, disease development was observed. SBPH-infested rice plants harbouring RBSDV began to exhibit plant growth abnormalities, such as dwarfism and leaf darkening. At 20 days-post-transplantation (dpt), these plants showed more serious developmental abnormalities. At 60 dpt, SBPH-infested rice plants infected with RSV showed leaf yellowing, stripe, chlorosis, and slower plant growth (Fig. 1A). Leaves were collected from these symptomatic rice plants for RT-PCR and western blotting analyses to detect the target viruses (Fig. 1B, C). The specific pairs of primers corresponding to RSV and RBSDV are shown in Additional file 2: Table S1. RT-PCR showed a specific band with the expected size appeared in RSV- and RBSDV-infected rice plants compared with mock-treated plants, respectively. These results and those of western blots using RSV NS3 and RBSDV p10 antisera (Fig. 1B, C) confirmed the independent infection of the plants by the viruses following inoculation.

Fig. 1
figure1

Flow chart of the investigation of m6A methylation during infection of rice by RBSDV or RSV. A Symptoms of RBSDV- and RSV-infected rice plants in plastic buckets at 60 days post infection (dpt). The left plant is a mock-treated plant, the middle two plants are infected with RBSDV, and the plant on the right is infected with RSV. B RT-PCR and western blot (WB) detection of RSV using a specific pair of primers corresponding to RdRp and anti-NS3 specific antiserum. Total proteins were stained with Coomassie brilliant blue (CBB), which was treated as the loading control. C RT-PCR and WB were performed to detect RBSDV in rice using a specific pair of primers corresponding to CP. Anti-p10 specific antiserum was carried out for the WB. D Experimental flow chart of m6A-IP-seq and RNA-seq using RBSDV- and RSV-infected plants. NGS, next-generation sequencing. MeRIP, methylated RNA immunoprecipitation

Transcriptome-wide mapping of m6A in rice

To obtain a transcriptome m6A methylation modification map in rice, the mock-, RSV-, and RBSDV-infected rice samples were used for further analyses. The m6A analysis procedures were described, and the samples that were used for input (non-IP control), m6A-based RNA Immunoprecipitation (m6A-IP), and RNA-seq were clearly indicated (Fig. 1D). A series of IP, input, and mRNA libraries were constructed and sequenced, respectively (Additional file 2: Table S2). Samples of these series libraries came from mock-, RSV-, and RBSDV-infected rice leaves at 60 dpt. Each treatment was performed with two biological replicates. Raw sequencing data were further processed for adaptor and low-quality base removal. The obtained clear reads were aligned to the rice reference genome (Oryza sativa. IRGSP-1.0). Read distribution analysis showed that all the m6A-IP samples were highly enriched around the stop codon and within 3′-untranslated regions (3′-UTRs), which are in line with the previous reports in HIV-infected T cells and maize and suggested the m6A-IP sequencing data are reliable and with a high authenticity [6, 15] (Additional file 1: Fig. S1). The m6A-IP-seq analyses detected more than 26,000 m6A peaks in each individual treatment and biological replicate (Fig. 4A, Additional file 1: Fig. S2). For each treatment (one individual virus infection), high-confidence peaks were identified (Additional file 2: Table S3). Briefly, the regions that overlapped in at least one of the two replicates were designated high-confidence m6A peak regions. Confident peaks from different experimental conditions were further integrated into a unique m6A peak map. Consequently, a total of 26,390, 27,038, and 26,675 unique m6A peaks with high confidence (p < 0.05, fold change > 1.5) for mock-, RBSDV-, and RSV-infected samples were detected, respectively (Fig. 4A, Additional file 1: Fig. S2). After comparison with mock treatment, the RBSDV- and RSV-infected samples displayed 8011 and 6603 different regulated peaks, accounting for an average of approximately 1 m6A peak within transcription units from each gene. Among these differential m6A peaks, there are 3897 and 2900 new peaks appeared upon RBSDV and RSV infection, and 4113, and 3702 common peaks in RBSDV- and RSV-infected sample compared with mock-treated rice, and 1503 differential m6A peaks were both appeared in RBSDV- and RSV-infected sample (Additional file 2: Table S4). These results suggested that 48.7% and 43.9% differential m6A peaks newly appeared upon RBSDV- and RSV infection of rice, respectively, which also indicated that m6A methylation was tightly associated with viruses’ infection of the plant. At the genomic level, these unique m6A-methylated peaks for the three treatments were unevenly distributed across each rice chromosome. The common peak density was also mapped. The gene density according to the previously reported data is presented in Fig. 2, and the m6A peak distribution density was highly consistent with the corresponding gene density on the same chromosome position in the mock sample.

Fig. 2
figure2

Circos plots of the m6A methylome in rice plants infected with RBSDV or RSV. The six rings from the outside to the inside show the genomic positions (1st), gene density (2nd), peak density of mock-treated rice plants (3rd), peak density of RSV-treated rice plants (4th), peak density of RBSDV-treated rice plants (5th), and the common peak density of RBSDV- and RSV-infected rice plants (6th)

Widespread m6A methylation of RSV and RBSDV genomic RNA

The m6A-IP experiment was performed twice. The peak calling method used was stringent (false discovery rate < 0.01). The two replications of the next generation sequencing (NGS) data revealed high correlations with the bound RNAs (0.990), which indicated the high replicability of the sequencing results. The clear reads obtained by NGS were also aligned to the reference RBSDV and RSV genomic RNAs, and the m6A peaks spanning the full sequences of different segments of viruses were mapped (Fig. 3, Additional file 2: Table S5). In particular, clusters of m6A peaks were clearly observed in the 5′ terminal of RBSDV genomic S1, S2, S3, S4, S5, S6, S9, and S10, and some discrete peaks appeared in S4, and S7 (Fig. 3A, red arrows). In RSV-infected sample, the main m6A peak clusters were located on the genomic RNA2, RNA3, and RNA4, and compared to the input, several clearly m6A peaks were in the RNA1 to RNA4, two m6A peaks located to the 3′ terminal of RNA1 (Fig. 3B, red arrows). We also exhibited the fine m6A peaks that distributed to the each viral genomic RNAs (Additional file 2: Table S5), and the viral-specific m6A peaks that distributed on RBSDV S5, S6, and S9 (Fig. 3C) and RSV RNA1, RNA2, and RNA3 (Fig. 3D) were selected and exhibited. Our results suggested that the m6A modification often occurred in the 5′-terminal of the genomic RNAs of RBSDV, while they were random distributed on the RSV genome. These maybe resulted from the characteristic of the two completely different plant viruses (RSV is a single-strand RNA virus, while RBSDV is a double-stranded RNA virus). Taken together, in addition to the mRNA of the host plant, viral mRNAs could also be N6-methyladenosine methylated under interactions between virus and rice, and the m6A distribution pattern on viral genomic RNA was specific and novel.

Fig. 3
figure3

Circos plots of the m6A methylome in RBSDV and RSV genomic RNAs. A Distribution of m6A methylated reads on the ten RBSDV genomic RNAs. Six rings from outside to inside show genomic positions (1st), reads distribution of RBSDV_1_Input (2nd), reads distribution of RBSDV_1_IP (3rd), reads distribution of RBSDV_2_Input (4th), reads distribution of RBSDV_2_IP (5th), and the GC content of the genomic RNA (6th). B Distribution of m6A methylated reads on the four RSV genomic RNAs. Six outer rings were indicated similarly to RBSDV above. C m6A methylation peaks in the full-length RBSDV segment 5 (upper panel), 6 (middle panel), and 9 (bottom panel). The detail peaks regions and viral gene annotation are shown in the Additional file 2: Table S6. Top numbers show the full length of the analysed RNA segments, and bp is the short name of base-pair. Blue colour marked line shows the m6A peak region on viral genome, and number 1 and 2 mean the two replicate of the m6A-IP-sequencing. D Distribution of m6A peaks on the RSV genomic RNA1 (upper panel), RNA2 (middle panel), and RNA4 (bottom panel). The detail peak regions and viral gene annotation are shown in the Additional file 2: Table S6. Top numbers show the full-length of the analysed genomic RNAs, and the nt means nucleotide. Other marks are similar with Fig. 3C

Activation of rice m6A RNA methylation levels upon virus infection

To explore the changes of m6A RNA methylation levels in virus-infected rice, the sequencing data were analysed. Collectively, there were 15,977, 16,854, and 16,267 m6A methylated genes that corresponding to mock-, RBSDV-, and RSV-infected rice, respectively (Additional file 2: Table S6). The findings clearly indicated that rice m6A RNA methylation was enriched under RSV and RBSDV infection. In terms of differential m6A peak number, the enriched peaks were also increased in rice infected with viruses (Fig. 4A, Additional file 2: Table S6). Meanwhile, the confidence of the peak significances was determined by calculating the correlations of the peak numbers and the p value (Fig. 4B). The horizontal ordinate dimension (-log10[q-value]) of most of the m6A peaks ranged from 4 to 10. The findings indicated that most of the peaks were highly confident and credible. The differential m6A peaks were selected and compared to those obtained for mock-treated rice plants. There were 8,010 and 6,602 unique and different m6A peaks for RBSDV- and RSV-infected rice plants compared to the mock-treated sample, respectively (Fig. 4C). For the digital exhibition of the m6A signal intensity, a violin plot was used to show the fold enrichment distribution (Fig. 4D). The median number of the peak enrichment-fold was approximately 40 and 76 (Additional file 2: Table S7). To explore the confidence of the different peaks deposited in RBSDV- and RSV-infected samples, the significance distribution of different m6A peaks was analysed by correlation with the p- value of each peak and peak number (Fig. 4E). The main differential regulated m6A peaks were deposited in 3 on the horizontal axis. The finding indicated that the different m6A peaks of rice mRNA under viral infection was confident and reliable. Taken together, the m6A modification levels of rice mRNAs were enriched under infection by plant viruses.

Fig. 4
figure4

Activation of rice m6A methylation by virus infection. A Histograms show the number of unique peak in mock-treated and RBSDV- and RSV-infected rice plants. The Y-axis represents the peak number, and the X-axis represents the treatments. B Significant distribution analysis of the different peaks in mock-treated and RBSDV- and RSV-infected rice samples. The Y-axis represents the peak number, and the X-axis represents the negative value of the logarithm of the p- value base-10. C Histogram showing the different regulated numbers of m6A methylated genes in RBSDV- and RSV-infected rice compared with mock-treated rice samples. D Comparisons of the fold change of the different regulated peaks in RBSDV- and RSV-infected rice plants. The Y-axis represents the logarithmic of peak folded-enrichment base-2. E Significant distribution analysis of the different peaks in RBSDV- and RSV-infected rice samples. The Y-axis represents the numbers of different regulated peaks, and the X-axis represents the negative value of the logarithmic value of the p- value base-10

Association of m6A with genes that were not actively expressed in virus-infected rice plants

RNA deep sequencing (input samples) of the 60 dpt rice plants of mock-, RBSDV-, and RSV-infected samples with two biological replicates were performed, and each replicate was used to examine the correlation between m6A modification and gene expression in rice. The euclidian distance (ED) of the two replicates’ reads was small, and the square colour of the two was close (Fig. 5A). The findings suggested that the experiments were reproducible between the two replicates. Based on the confident sequencing results, the m6A peaks of each treatment were annotated to the genes. Rice genes that could be annotated by the m6A peaks were called m6A gene, and that couldn’t be annotated were called non-m6A gene in the next analyses. According to the Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) (Additional file 2: Table S6), these genes were divided into three groups (FPKM < 1, 1 < FPKM < 5, and FPKM > 5), and the ratio of gene number of each category were calculated and showed by the heatmap (Fig. 5B). Either in the mock-treated sample or the viruses’ infected rice, most of the analysed genes (m6A or non-m6A genes) were mainly distributed in the low expressed category (FRKM < 1) (Fig. 5B). Compared to the mock-treated sample, the ratio of m6A methylated genes were increased in low expressed category of RSV- and RBSDV-infected samples (Fig. 5B). When these genes (m6A and non-m6A) were divided into highly expressed genes (FRKM ≥ 1) and genes that were not highly expressed (FRKM < 1), the expression of genes that was high or low were showed in mock-, RBSDV-, and RSV-infected samples (Fig. 5C), and we found that most genes were non-methylated and distributed in low expression category in mock-, RBSDV-, and RSV-infected rice samples (Fig. 5C). The m6A methylated genes were slightly enriched in the low expression category upon viruses’ infection of rice compared to the mock-treated sample (Fig. 5C). Either in mock-treated or in RSV-infected and RBSDV-infected samples, the m6A-methylated genes were faintly expressed at a higher level than those of non-m6A methylated genes (Fig. 5D), and the highly expressed genes displayed a lower m6A occupancy (Fig. 5B). Hence, the m6A methylation mainly occurred in low expressed genes either in mock-treated or in RSV- and RBSDV-infected rice samples and the m6A modified gene number were slightly enriched in the low expressed category upon viruses’ infection.

Fig. 5
figure5

Analyses of the relationship between m6A methylation with expression levels of the target gene in rice under plant virus infection. A Euclidian distance (ED) coefficients among gene expression profiles generated by RNA-seq analysis of the two biological replicates of the three treatments. RNA-seq was performed simultaneously with m6A-IP-seq with total RNA extracted at 60 dpt. A lower value means a closer ED of the two compared objects, and with higher reproducibility of the two replicates. B The percentage of rice m6A methylated and un-methylated genes at a defined FPKM levels (< 1, 1–5, and > 5). Different colour densities indicate different percentages of the corresponding gene. C Comparisons of number of non-m6A methylated genes and number of m6A methylated genes in their gene bodies with high (FPKM > 1) and low (FPKM > 1) expression levels in the three treatments. Relationships between gene expression and number were tested using the chi-square test. “*” indicate p- value < 0.05, and “**” means p- value < 0.01. D Box plot comparing FPKM expression levels between non-m6A methylated genes and m6A methylated genes in the three treatments. A two-tailed unpaired Student’s t -test was performed to calculate the p- values of these three treatments

Different regulated m6A-methylated genes in virus-infected rice plants

To investigate which pathway-related genes were m6A methylated, the obtained different m6A peaks were mapped to the rice reference genome. The gene was divided into 5′-UTR, 3′-UTR, first exon region, and exon excluding the 5′-UTR and 3′-UTR in the protein-coding region. The sequenced clear m6A data were aligned to these functional elements. We found that 38%, 42%, and 47% of the m6A modification sites were located on the 3′-UTR of the genes in RBSDV-infected, RSV-infected, and common RBSDV and RSV peaks, respectively (Fig. 6A). Furthermore, 23%, 23%, and 24% of the m6A modification sites were deposited on the 5′-UTR regions of the aforementioned genes (Fig. 6A). These results indicated that most of the m6A sites were located in the non-coding region of genes. Gene ontology (GO) analyses were performed to provide insight into the biological functions of different m6A modifications in rice. The methylated genes were involved in multiple molecular functions, particularly in binding and catalytic activity (Fig. 6B). Thus, m6A modifications may act as epigenetic markers that mediate molecular interactions between rice and plant viruses. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the different m6A modification-related genes were performed. The five enriched pathways were carbohydrate metabolism, translation, protein folding/sorting/degradation, amino acid metabolism, and signal transduction (Fig. 6C). these results suggested that m6A modification was widely involved in and closely related to the intracellular carbohydrate metabolism, amino acid metabolism, and protein translation/folding/sorting/degradation in RSV- and RBSDV- infected rice sample.

Fig. 6
figure6

GO and KEGG analyses of different m6A methylated genes in RBSDV- and RSV-infected rice plants. A Percentages of the different gene bodies (5′-UTR, 3′-UTR, 1st exon, and other exons) in different m6A methylated genes in rice infected with RBSDV alone, RSV alone, or with both. B GO analysis of the different m6A methylated genes in rice infected with RBSDV alone, RSV alone, or with both compared with mock-treated rice plant. C KEGG analysis of the different m6A methylated genes in rice infected with RBSDV alone, RSV alone, or with both

Analysis of consensus motifs related to m6A modifications

To investigate whether there are conserved motifs in clear sequencing reads of mock-, RBSDV-, and RSV-infected rice samples, the MEME suite software was used for de novo scanning of the enriched motifs. In total, 172,801, 152,136, and 168,666 raw m6A peaks appeared in the sequencing data corresponding to the mock-, RBSDV-, and RSV-infected rice samples (Additional file 2: Table S8). When the parameter was set as the fold enrichment (FC) > 1.5 (p < 0.05), 26,389, 27,038, and 26,675 significantly enriched m6A peaks were identified in mock-, RBSDV-, and RSV-infected rice samples, respectively. These significantly enriched peaks were performed to do de novo searching of the common consensus by software. The “CGBCGKC” (B = C, G, or T; K = G or T), “CGRCGVC” (R = A or G; V = A, G or C), and “CGHCGDCG” (H = A, C or T; D = A, G or T) motifs were enriched in mock-, RBSDV-, and RSV-infected rice m6A methylated reads, respectively, using the DREME suite (Fig. 7A, C, and E). MEME suite scanning revealed the enrichment of the longer conserved motifs “CSYCBCCGCCSYCGCCGCSSY”, “GCGGCGGCGRCGGCG”, and “YCGCCGCCGBCGCCG” in m6A methylated reads of mock-, RBSDV-, and RSV-infected rice (Fig. 7B, D, and F). Besides, the enriched top five consensus in mock-, RBSDV-, and RSV-infected rice were selected and exhibited (Additional file 2: Table S9), and results showed the dynamic of the ranking with viruses’ infection. The top four widely studied motifs in other species were selected and analyses using the DREME and MEME suites in the mock-, RBSDV-, and RSV-infected rice samples (Additional file 2: Table S10), and the results showed that the most common motifs still were “RRACH” (68%), and “URUAY” (28%). Taken together, these results implied that the common recognition mechanism of the m6A methylated sites during the plant-virus interaction, and viruses’ infection caused the dynamics of the top five consensus ranking.

Fig. 7
figure7

Identification of predominant consensus motifs containing m6A methylation sites in mock-treated and RBSDV- and RSV-infected rice plant using DREME and MEME suites. A Sequences logo representations of the consensus motifs containing m6A sites in mock-treated rice samples. B The most enriched consensus motif in RBSDV-infected samples. C The most enriched consensus motif in RSV-infected samples

Rice m6A methylation modifications display sensitive and dynamic responses to infection of rice plants by viruses

Abiotic and biotic stress, such as heat, cold, salt, and drought, and a variety of fungal, bacterial, viral, and nematode plant pathogens often affect the growth and development of rice plants, which poses a significant threat to the yield and restricts the global distribution of rice plants. In addition to the visual phenotype changes to these stresses, there are sophisticated and fine gene expression and regulation networks behind. Hence, we wondered whether mRNA m6A modification was involved in the interactions between rice and plant viruses. LC-MS/MS was used to track the m6A methylation dynamics of rice mRNA in plants infected with RBSDV or RSV. The difference in the potential epigenetic basis of the virus infection on rice plants was also investigated. Compared to the mock treatment, the N6-methyladenosine level was significantly increased in RBSDV- and RSV-infected rice compared to the mock-treated sample just at the one-day-post-infection (dpi), and the FC of the m6A content was further increased with the infection time extended from 2 days to 16 days (Fig. 8A). These results demonstrated that the m6A methylation was sensitive and dynamic upon viruses’ infection. The FC of the m6A content in the RBSDV-infected samples was approximately 5 times higher than that in the mock-treated sample at 16 days, and 3 times higher in the RSV-infected sample (Fig. 8A). Total RNAs were extracted, then the mRNAs were isolated. Nucleic acid-based dot-blot analyses were performed. The blotting results showed that with the increase in mRNA loading, the colour reaction of virus-infected rice plants became stronger compared with the mock-treated samples (Fig. 8B). These findings indicated that the infection of plants with viruses markedly increased the N6-methyladenosine level of plant mRNAs, and the N6-methyladenosine level was tightly correlated with virus infection of rice plants.

Fig. 8
figure8

Rice m6A methylation levels are positively associated with the expression of key genes involved in antiviral RNA silencing pathways and plant hormone signals. A Comparisons of m6A methylation levels of the respect mock-treated and RBSDV- and RSV-infected rice plants at 0, 1, 2, 4, 8, and 16 dpi by LC-MS/MS. Error bars indicate mean ± SD, with three biological replicates. B Dot-blot analysis of m6A levels in extracted total RNA from samples at 16 dpi using the specific anti-m6A antibodies. The left side of the membrane depicts the amount of loaded mRNA from mock-treated and RBSDV- and RSV-infected rice plants, respectively. C qRT-PCR analysis of the relative expression of OsAGO18 in mock-treated and RBSDV- and RSV-infected rice plants at 0, 1, 2, 4, 8, and 16 dpi. D Relative expression of the OsSLRL1 in the three treatments at 0, 1, 2, 4, 8, and 16 dpi. E Analysis of m6A methylation levels on different fragments of OsAGO18 by m6A-IP-qPCR. The upper panel indicates the gene structures of OsAGO18 labelled with fragments amplified in the m6A-IP-qPCR assay. The results of positions 1 and 12 were chosen for figure exhibition. F m6A-IP-qPCR assay of the m6A methylation levels of different fragments on OsSLRL1. Similarly, the upper panel represents the gene structures labelled with fragments amplified in the m6A-IP-qPCR analyses. Results of positions 2, 3, and 4 were selected for the display in the figure. Error bars denote mean ± SD, n = 3 biological replicates in all qRT-PCR assays

Correlation of m6A levels with expression of key genes involved in plant antiviral RNA silencing pathway and hormone metabolism

The m6A methylation level of rice total mRNAs was strongly correlated with virus infection of the plants (Fig. 8B). This prompted us to explore the potential molecular mechanisms. We focused on the virus infection of rice and explored whether m6A modifications were involved in the regulation of the key genes’ expression that related to virus infection. Based on the m6A-IP sequencing results (Additional file 2: Table S11), the m6A methylated candidate genes, ARGONAUTE 18 (OsAGO18), and SLENDER RICE LIKE 1 (OsSLRL1), were selected. OsAGO18 was reported to participate in anti-virus RNA silencing pathways in rice by binding small interfering RNAs (siRNAs) [39]. The expression of OsAGO18 can be activated by jasmonate (JA) signalling [40]. OsSLRL1 is a DELLA family protein that regulates plant hormone metabolism and mediates the growth and development of plants [41,42,43]. qRT-PCR was performed and the relative expression levels of the selected candidate genes were determined in plants infected by the virus from 0 to 16 days. OsAGO18 was upregulated by 14-fold and 4-fold in RBSDV- and RSV-infected rice plants, respectively, compared with mock-treated plants, while the OsSLRL1 was respectively downregulated by 12-fold and 25-fold in RBSDV- and RSV-infected rice plants, respectively (Fig. 8C, D). Based on these results, we concluded that the antiviral RNA silencing pathway was activated, and the synthesis and degradation of plant hormones controlling the growth, development, and basal resistance to pathogens were regulated in virus-infected rice plants.

To validate whether the m6A modifications were correlated with the expression of key genes, m6A methylation-based RNA immunoprecipitation and qRT-PCR technology (m6A-IP-qPCR) was performed. The m6A levels in position 12 of OsAGO18, not position 1 or other regions, was increased by approximately 13-fold and 11-fold in RBSDV- and RSV-infected samples, respectively, compared with mock-infected rice plants (Fig. 8E). The m6A levels in positions 2 and 4 of OsSLRL1 were increased by 16-fold and 12-fold in RBSDV-infected rice plants, respectively, but not in position 3 (approximately 6-fold and not responsive to virus infection) (Fig. 8F). For RSV-infected samples, m6A levels in positions 2 and 4 of OsSLRL1 were not significantly changed (approximately 6- and 4-fold, respectively) under virus infection, whereas the m6A level was markedly up-regulated at position 3 (Fig. 8F). The collective findings indicated that changes in the relative expression of OsAGO18 and OsSLRL1 tightly correlated with the changes in their m6A modification level and that the m6A methylation of different gene regions may have different effects on its expression. These findings could contribute to a deeper understanding of the molecular basis of interactions between rice and viruses.

Involvement of rice m6A modification in the regulation of the expression of the main m6A modification machinery components in virus-infected plants

RNA m6A modifications occur in many eukaryotes, and it was done by the m6A methylation machinery. The WRITER, READER, ERASER, and methyl synthetases that produce the DONOR (methyl) are the four main components of the m6A methylation machinery. To investigate whether the gene expression of the main components of m6A methylation machinery can be affected by m6A modifications, qRT-PCR was performed. The relative gene expression levels of five WRITER genes (OsMAT1, OsMAT2, OsMAT3, OsMAT4, and OsFIP), five ERASER genes (OsALKBH10B, OsALKBH9B-1, Os05g0401500, OsALKBH9B-2, and Os03g0238800), 12 READER genes (OsYTH1–12), and five S-adenosyl-I-methionine synthetases (OsSAM1, OsSAM1L, OsSAM2, OsSAM2L, and OsSAM3) were determined in mock and virus-infected rice samples. The relationship between the relative expression level of target genes and m6A modification sites was counted by combination analyses (Additional file 2: Table S12), and no any fixed regular pattern has been found.

Compared with mock-treated samples, the OsMAT3 and OsMAT4 WRITER genes were significantly increased in RSV-infected samples, whereas they did not change in RBSDV-infected samples (Fig. 9A). For ERASER genes, the expression of OsALKBH10 was significantly suppressed in both the RBSDV- and RSV-infected samples, both Os05g0401500 and OsALKBH9B-2 were up-regulated in RSV-infected samples, whereas the expression of OsALKBH9B-2 was increased in RRSDV-infected samples (Fig. 9B). For READER genes, the expression levels of OsYTH1, OsYTH3, OsYTH5, and OsYTH7 were up-regulated in RSV-infected samples, whereas OsYTH8 was significantly decreased. Only the expression level of OsYTH6 was increased in RBSDV-infected samples (Fig. 9C). For the methyl donor producer, the OsSAM1 and OsSAM2 expression levels were significantly down-regulated in both RBSDV- and RSV-infected samples, whereas the expression of OsSAM2L and OsSAM3 in RSV-infected samples were up-regulated, and the expression of OsSAM3 was also markedly increased in RBSDV-infected samples (Fig. 9D). Based on the relative expression results, we integrated the m6A-IP sequencing data with the m6A methylation machinery in rice (Additional file 2: Table S12). Modification of m6A occurs in genes of m6A methylation machinery under plant viruses’ infection. For the donor producer, OsSAM2 was m6A methylated in both RBSDV- and RSV-infected samples. For the WRITER genes, OsMTA3 and OsMTA4 were m6A methylated in RBSDV- and RSV-infected samples, respectively. For ERASER genes, the m6A modification occurred on OsALKBH10B and OsALKBH9B in RSV-infected samples, whereas no ERASER genes were m6A methylated in RBSDV-infected samples. For READER genes, OsYTH01, OsYTH10, OsYTH11, and OsYTH12 in RBSDV-infected samples, and OsYTH05 and OsYTH08 in RSV-infected samples were m6A methylated (Fig. 9E). In summary, the m6A modifications occurred in the genes that encoding the plant m6A methylation machinery, and probably regulated the dynamics of the target gene expression, which may act as a main post-translational gene expression regulation strategy under plant virus infection.

Fig. 9
figure9

Integrated analyses of the main components of m6A methylation machinery in rice with m6A methylation modifications and gene expression in plants infected with viruses. A Relative expression levels of five “WRITER” components in rice plants by qRT-PCR analyses. B qRT-PCR analysis of the relative expression of five “ERASER” components in rice plants. C Relative expression levels of twelve “READER” component genes were determined by qRT-PCR. D Relative expression levels of five methyl “DONER” synthesis genes were analysed by qRT-PCR. E Rice m6A methylation pathways and related m6A methylated genes under plant virus infections. Blue coloured letters indicate the m6A methylated genes of certain treatments. For instance, RBSDV: OsMTA3, means the OsMTA3 gene was methylated in RBSDV-infected sample. All qRT-PCR assays were performed with three biological replicates, and the error bars denote the mean ± SD

Involvement of rice m6A modification in regulating the expression of the main antiviral RNA silencing components in virus-infected rice plants

To investigate whether m6A methylation modification was involved in the regulation of the main host antiviral RNA silencing pathways genes and the reported resistance genes corresponding to these two viruses. Genes, which included nine DCLs, five RNA-dependent RNA polymerases (RDR), 17 ARGONAUTE genes (OsAGOs), seven RBSDV resistance genes (OsGDI homologous), and two RSV resistance genes (OsSOT1 and OsStvb-i), were selected for relative expression analyses. Integrated analyses of the antiviral genes with m6A modification sites on the gene body and relative expression level were determined in virus-infected rice plants.

qRT-PCR revealed that both the OsDCL2b-1 and OsDCL2b-2 were significantly up-regulated in RBSDV- and RSV-infected samples than in mock-treated samples (Fig. 10A). In addition, the expression of OsDCL1b was greatly suppressed in RBSDV-infected samples, whereas it was increased in RSV-infected samples (Fig. 10A). Compared to the mock-treated samples, the expression of OsRDR1 and OsRDR3 in both RBSDV- and RSV-infected samples was markedly increased, and OsRDR2 was upregulated only in RSV-infected samples (Fig. 10B). In RBSDV-infected samples, only the expression of OsAGO18 was up-regulated significantly (Fig. 8C), whereas in the RSV-infected samples, OsAGO5c, OsAGO12, OsAGO13, OsAGO14, and OsAGO18 were markedly increased (Figs. 8C and 10C). Concerning resistance genes, seven ZmRabGDI (Rab GTPase dissociation inhibiter) homologous genes in rice, which were defined as RBSDV resistance genes on maize [44], were screened out. The RSV resistance genes OsSOT1 (Sulfotransferase 1) [45] and OsStvb-i (Stripe disease resistance i) [46], were also chosen. A total of nine genes were selected for further gene expression analyses. OsGDI-1-1, OsGDI-1-2, OsGDI-2, and OsSOT1 were significantly decreased both in RBSDV- and RSV-infected rice plants compared to that in the mock-treated sample (Fig. 10D). In contrast, OsGDIα was markedly upregulated in RBSDV- and RSV-infected rice plants (Fig. 10D). The collective findings indicated that OsDCL2b-1, OsDCL2b-2, OsRDR1, OsRDR3, OsAGO12, OsAGO13, and OsAGO18 were the main antiviral genes in the innate immune RNA silencing pathways, and that OsGDIα, OsGDI-1-1, OsGDI-1-2, OsGDI-2, and OsSOT1 may act as resistance genes for both RBSDV and RSV.

Fig. 10
figure10

Integrated analyses of main antiviral RNA silencing pathway-related genes with m6A methylation modifications and gene expression levels in rice infected with viruses. A Relative expression levels of nine OsDCL genes in mock-treated and RBSDV- and RSV-infected samples. B Relative expression levels of five OsRDR genes in mock-treated and RBSDV- and RSV-infected samples using qRT-PCR analyses. C Relative expression levels of 17 OsAGO genes were determined with respect to mock-treated and RBSDV- and RSV-infected samples using qRT-PCR analyses. D Relative expression levels of nine resistance genes, including seven OsGDI genes, one OsSOT1 gene, and one OsStvb-i were determined in mock-treated and RBSDV- and RSV-infected samples using qRT-PCR analyses. E Rice antiviral RNA silencing pathways and the related m6A methylated genes in virus-infected plants. Blue coloured letters depict the m6A methylated genes of a certain treatment, as detailed in Fig. 9E. All qRT-PCR assays were performed with three biological replicates. The error bars denote the mean ± SD

To analyse the potential roles of m6A methylation in the expression of 57 selected genes in RNA silencing pathways and viral resistance genes, we systematically analysed their m6A modification status according to the results of m6A-IP-seq (Additional file 2: Table S13). For the main RNA silencing pathway genes, OsAGO1c, OsAGO2, OsAGO18, and OsRDR1 were m6A methylated in RBSDV-infected rice. In RSV-infected samples, OsDCL1a, OsDCL3b, OsDCL4, OsAGO2, OsAGO12, OsAGO5c, OsAGO17, OsAGO18, OsRDR1, and OsRDR3 were significantly methylated (Fig. 10E). In contrast, no viral resistance genes were m6A methylated (Fig. 10E). These results indicated the possible involvement of m6A methylation in the main antiviral RNA silencing pathways. The methylation may act as a fine regulator to mediate the spatial and temporal expression of target genes in arm race of plant-virus interaction.

Integrated analyses of relative expression profile and m6A modification of seven phytohormone metabolism-related genes

Viruses have various strategies to reprogram the host’s cellular status to one that is prone to viral replication and spread. In plants, the specific environment created by viruses usually refers to phytohormone regulations of nearly all aspects of plant physiology, including development, growth, defence, and reproduction [47]. These phytohormones, which include jasmine acid (JA), salicylic acid (SA), abscisic acid (ABA), auxin, cytokinin (CTK), ethylene (ET), and brassinosteroids (BR), have significant roles in plant development and physiological regulations and are also involved in defence against pathogens [48, 49]. Previous findings suggest that phytohormones are strongly associated with virus infection and symptom development. We investigated whether these phytohormones are regulated by m6A modification and whether the expression profiles of their metabolic genes are altered by virus infection of rice plants. qRT-PCR was performed to determine the relative expression of the seven main phytohormone metabolism-related genes (Additional file 1: Fig. S3–S9). Further, integrated analyses of their relative expression and m6A modification sites were performed (Additional file 1: Fig. S10, Additional file 2: Table S14).

For the JA pathway, the biosynthesis genes OsLOX8 and OsLOX9 were markedly up-regulated in both RSV- and RBSDV-infected plants. OsJMT1-1, OsAOS2, and OsLOX2 were markedly decreased, whereas OsAOS1 and OsLOX5 were not changed (Additional file 1: Fig. S3A). In addition, OsJMT1, OsLOX1, and OsHPL3 were decreased in RSV-infected samples, whereas OsJMT1 was increased in RBSDV-infected plants (Additional file 1: Fig. S3A). The expression profiles of 21 JA responsive genes were further investigated. OsPR1a, OsWRKY28, OsPR2, and OsPR5-4 were significantly increased in both RSV- and RBSDV-infected plants, OsPR5-3, OsMYB2, OsMYB55/56-1, OsWRKY10, OsRbohA, OsRbohB, and OsRbohC were markedly decreased, and OsPR1b, OsPR5-2, OsbZIP52, and OsRbohE were not changed. OsPR5-1 and OsPR1 were increased in RBSDV- and RSV-infected rice, respectively (Additional file 1: Fig. S3B and S3C). These results indicated that the activated JA pathways mainly depend on the markedly increased expression of biosynthesis-related OsLOX8 and OsLOX9, and the up-regulated OsPR1a, OsWRKY28, OsPR2, and OsPR5-4 responsive genes.

For the SA pathway, we investigated the biosynthesis-related genes OsICS1, OsPAL, OsPAL1, OsAIM1, OsCM, and OsEDS1) and the response genes (OsPR1-101, OsWRKY45-1, OsWRKY45-2, OsSGT1, OsPR1b, OsPR1a, OsPR1-12, OsPR1-21, OsPR1-22, OsPR1-51, and OsPR1-121). OsICS1, OsPAL, and OsPAL1 were significantly suppressed in both the RBSDV- and RSV-infected plants, OsCM and OsEDS1 were not altered, and only OsAIM was slightly up-regulated in RSV-infected plants (Additional file 1: Fig. S4A). These results indicated that SA biosynthesis was suppressed upon virus infection. OsWRKY45-1, OsWRKY45-2, and OsRR1-51 were significantly increased in both the RBSDV- and RSV-infected plants, whereas the OsSGT1, OsPR1-21, OsPR1-22, and OsPR1-121 were dramatically decreased. OsPR1-101 and OsPR1a were markedly up-regulated in RSV-infected plants but significantly down-regulated in RBSDV-infected plants (Additional file 1: Fig. S4B and S4C). These results suggested that the SA pathway was suppressed by virus infection through down-regulation of the main biosynthesis genes.

For the ABA pathway, qRT-PCR revealed that OsNCED3 was significantly up-regulated in RBSDV-infected samples, whereas in RSV-infected plants, OsNCED1-1 was up-regulated. OsABA was suppressed in both RBSDV- and RSV-infected plants. Other biosynthesis-related genes, including OsABA1, OsbZIP72, and OsbZIP23, were unchanged in RBSDV- and RSV-infected plants (Additional file 1: Fig. S5A). ABA deactivation genes OsABA8OX1, OsABA8OX2, and OsABA8OX3 were significantly up-regulated in RBSDV-infected plants. In RSV-infected plants, OsABA8OX2 and OsABA8OX3 were slightly decreased (Additional file 1: Fig. S5B). The collective results indicated that the ABA pathway was activated by RBSDV infection, but was suppressed in RSV-infected plants.

For the Auxin pathway, the relative expression levels of auxin-related metabolic genes were affected differently by virus infection. Of the analysed biosynthesis genes, OsYUCCA1, OsYUCCA5, OsYUCCA6, OsYUCCA9, OsAO-2, OsAO3-L, and OsAOO3-2 were significantly down-regulated upon RBSDV and RSV infection, whereas OsYUCCA8, OsAO-1, OsAAO, and OsAAO3-1 were up-regulated (Additional file 1: Fig. S6A). The OsGH3.8 and OsGH3.2 auxin transformation genes were markedly increased in both RBSDV- and RSV-infected plants (Additional file 1: Fig. S6B). The OsPIN1A, OsPIN1B, OsPIN2, OsPILS7a, and OsPILS7b auxin transportation genes were markedly down-regulated in RBSDV- and RSV-infected plants (Additional file 1: Fig. S6C), while only OsPIN3A was increased (Additional file 1: Fig. S6D). The OsIAA3 and OsIAA7 signalling genes were markedly decreased in both RBSDV- and RSV-infected plants (Additional file 1: Fig. S6E). The other signalling gene (OsIAA20) was significantly up-regulated in RBSDV-infected samples, whereas it was severely suppressed in RSV-infected plants (Additional file 1: Fig. S6E). These results showed that the auxin pathway was suppressed in both RBSDV- and RSV-infected plants.

For the CTK pathway, genes responsible for CTK biosynthesis, including OsIPT3 and OsIPT7, were dramatically decreased in both RBSDV- and RSV-infected plants (Additional file 1: Fig. S7A). In contrast, the CTK transformation genes or oxidase genes (OsCKX4 and OscZOGT1) were significantly up-regulated in both RBSDV- and RSV-infected plants (Additional file 1: Fig. S7B). The OsPR8, OsPR10a-1, OsPR1a-1, OsPR1b, and OsPR10a-2 CTK responsive genes were markedly increased in both RBSDV- and RSV-infected plants (Additional file 1: Fig. S7C & S7D). Most of the analysed genes (OsPR2-1, OsPR2-2, OsPR3, OsPR4, OsPR5-1, and OsPR6) were suppressed (Additional file 1: Fig. S7D). Other responsive genes, which included OsPR1a-2, were dramatically upregulated in RBSDV-infected samples, but not in RSV-infected samples. The expression profile of OsPR9 was completely opposite (Additional file 1: Fig. S7D). These results showed that the CTK pathway in rice was suppressed by viral infection.

For the ET pathway, the OsACS2 biosynthesis gene was significantly up-regulated in both RBSDV- and RSV-infected plants, as was OsACS1 in RBSDV-infected plants (Additional file 1: Fig. S8). Most of the oxidase genes (OsACO7, OsACO1, and OsACO2) were dramatically decreased in both the RBSDV- and RSV-infected plants (Additional file 1: Fig. S8). The collective results indicated that the ET pathway was activated in rice upon virus infection.

For the BR pathway, 24 genes, including 14 biosynthesis-related genes and 10 signalling genes, were analysed. Most of the biosynthesis genes (OsD11, OsDS11-L, OsCPD1, OsGSK2, OsRAV1, OsRAV2, OsRAVL1, and OsBZR1) were markedly suppressed in both RBSDV- and RSV-infected plants, whereas OsDWARF4 was significantly up-regulated (Additional file 1: Fig. S9A). The four signalling genes (OsBRI1-1, OsBRI1-2, OsBAK1-4, and OsBAK1-10) were dramatically decreased in RBSDV- and RSV-infected plants. The other four genes (OsI-BAK1, OsBAK1-2, OsBAK1-3, and OsBAK1-8) were severely downregulated, and two genes (OsBAK1-6 and OsBAK1-9) did not change (Additional file 1: Fig. S9B). The expression profiles of these genes were synchronously changed in the RBSDV- and RSV-infected plants. These results suggested that the BR pathway was suppressed in rice plants upon virus infection.

Seven phytohormone metabolism-related genes, including 154 genes, were analysed by qRT-PCR (Additional file 1: Fig. S10A). Of these genes, m6A modification was detected in 37 and 19 genes in RBSDV- and RSV-infected samples, respectively. m6A methylation was detected in nine genes in RBSDV- and RSV-infected plants (Additional file 1: Fig. S10B, Additional file 2: Table S14). m6A-IP-sequencing data mapped these m6A methylation sites to different regions of target genes. Most of the peaks were located in the coding sequence (CDS) region (Additional file 1: Fig. S10C). Furthermore, we analysed the relationship between m6A methylation sites on gene bodies and the relative gene expression levels. In most cases, the target gene was downregulated if the m6A site was located in the 5′-UTR, CDS, and CDS/3′-UTR, with the location in the 3′-UTR often being associated with up-regulation of the target gene (Additional file 1: Fig. S10D). Hence, a viral infection could regulate the expression of target genes, and the regulation mode was varied in different manners for the same genes under the two different viral infection conditions.

Widely integrated analyses of the relationship between m6A methylation regions and expression levels using RNA-seq and m6A-IP-seq

To investigate the relationship between the location sites of m6A methylation and the relative expression profiles, we analysed the common regulated peaks that appeared in both the RBSDV- and RSV-infected plants. We classified the reads that extended 5′-UTR to the start codon (5′-UTR/CDS), reads that extended the stop codon to the 3′-UTR (CDS/3′-UTR), and other reads that did not go beyond two obviously distinct regions to the corresponding 5′-UTR, CDS, or 3′-UTR (Additional file 2: Table S15).

In total, 3,734 and 3,336 peaks were analysed in RBSDV- and RSV-infected plants, respectively (Additional file 1: Fig. S10E and S10F). In RBSDV-infected plants, the relative expression level of 44.45% m6A modified genes was not changed, 29.73% were downregulated, and 25.82% were up-regulated. Most m6A sites were located in the CDS and 3′-UTR regions (Additional file 1: Fig. S10E). The number of genes that were upregulated, downregulated, and unchanged were increased gradually when the m6A methylation occurred in the 5′-UTR/CDS and CDS/3′-UTR region (Additional file 1: Fig. S10E). In RSV-infected plants, the majority of m6A modified genes (69.04%) were downregulated. Most m6A sites (49.34%) were located in the CDS region and 3′-UTR. CDS and 3′-UTR location of the m6A modified sites often indicated that the target genes were downregulated (Additional file 1: Fig. S10F), and the m6A methylation was mostly occurred in CDS and 3′-UTR region in both RSV- and RBSDV-infected samples. Taken together, these results indicated that the regulatory roles of post-transcriptional m6A modification differed upon RBSDV or RSV infection of rice plants and that certain m6A sites on specific genes may have specific functions. The common regular pattern between expression level and m6A methylation region, which are suitable for all m6A methylated genes, have not been found so far. This may be due to the fact that different viruses cause plant disease through hijacking different host’s signal pathways, and different genes own different sequences characteristics.

Discussion

The discovery of the m6A modification in different organisms has revealed a new and promising area of research in RNA epigenetics [4, 29, 31]. The understanding of the epigenetic roles of m6A modifications in eukaryotes has just begun. In our studies, to deepen the understanding of the fundamental roles of m6A in rice plants and the interactions between rice and viruses, we analysed the genome-wide m6A distributions of the transcriptome using the rice cultivar HD5 infected with RBSDV and RSV. We determined a series of genome-wide m6A distribution maps using m6A-IP data that corresponded to viruses and rice plants. These sites were identified by independent m6A-IP-qPCR and nucleic acid-based dot-blot assays using an m6A specific antibody and by quantitative LC-MS/MS. qRT-PCR comparison with transcriptome data permitted the comprehensive analysis of the relative expression levels of important antiviral pathway-related genes in rice. The m6A regulation mode was systematically analysed by correlating the positions of m6A modification on the gene body and the relative expression levels of target genes. These analyses have provided the first data on the genome-wide m6A modification distributions on rice, RSV, and RBSDV. Most of the consensus motifs were enriched in rice with different treatments. m6A modification was activated in virus-infected rice plants. m6A modification was strongly associated with genes in rice that were expressed at low levels upon virus infection. m6A modification plays essential roles in the regulation of the relative expression of the main antiviral pathway genes. The mode of regulation of gene expression by m6A modification was different for RBSDV and RSV infections.

High-throughput transcriptome-based genome-wide m6A distributions

m6A modification distribution maps of virus infection and the associations with host gene expression and virus infection have already been determined in mammalian cells. Differences in these maps upon infection with different viruses have been described [6, 11, 35, 50]. In plants, mapping m6A sites in the gene bodies of mRNAs has been investigated in Arabidopsis [17, 20, 32, 51,52,53]. However, little has been reported in maize [15] and rice [54]. Presently, genome-wide m6A distribution maps of rice leaves and in rice infected with RSV and RBSDV were obtained using high-throughput transcriptomic-based m6A-IP sequencing (Fig. 2). In total, we obtained > 20,000 unique m6A peaks in each treatment, which exceeds the numbers previously reported in callus and leaves (approximately 8000 and 14,000, respectively) [54]. Most of the peaks were common in mock-, RBSDV-, and RSV-infected rice, but some m6A peaks were randomly distributed and different in RBSDV- and RSV-infected rice, and the common m6A peaks were mainly distributed in rice chromosomes 1, 2, and 3 (Fig. 2). Besides, some common m6A peaks were simultaneously and specifically distributed in chromosomes 4, 5, and 10 in RBSDV- and RSV-infected rice (Fig. 2). In addition, the m6A peak density was high in telomeric regions in rice chromosomes 2, 3, 4, 9, and 10 (Fig. 2). These results imply that post-transcriptional m6A modification may play different roles in different viral infection conditions and chromatin conformation alternations. The chromatin state and transcription regulatory function of m6A modification in chromosome-associated regulation was revealed recently [4]. However, the detailed and real biological functions of m6A modification in plants remain largely unknown and require further study.

In animal viruses, the m6A modification of genomic RNA was determined several years ago. This modification plays significant roles in virus infection, replication, and gene expression processes [11, 33, 35, 50]. However, in plant viruses, the m6A modification distribution on viral genomic RNA is obscure. Only two plant viruses, AMV and TMV, are reportedly involved in m6A modification in Arabidopsis and tobacco, respectively, [33, 36]. In Arabidopsis, protein interactions between AMV CP and atALKBH98 (Eraser) resulted in successful infection, while the increased abundance of m6A in the AMV genomic RNA impaired systemic infection in conditions of AtALKBH98 suppression [23]. In tobacco, TMV infection reduced m6A levels by up-regulating the potential demethylase and down-regulation of the potential methylase [36]. The collective findings suggest that m6A has an important regulatory role in interactions between plants and viruses. However, the details remain unknown. We focused on the two most important and devastating viruses concerning rice production (RBSDV and RSV) to try to understand whether m6A was involved in the rice-virus interactions. The two genome-wide m6A distribution maps of RBSDV and RSV (Fig. 3) are novel. The precise and detailed positions of m6A sites on viral genomic RNA, especially RBSDV segment S5, S6, and S9 (Fig. 3C), or RSV RNA1, RNA2, and RNA3 (Fig. 3D). Although previous studies have shown that m6A modifications participate in plant virus infections, the precise m6A positions on the genomic RNA were not revealed. Our study is the first to show the genome-wide m6A distributions in the RBSDV and RSV. The candidate potential m6A sites may permit for a deeper understanding of the roles in virus infection or interactions between viruses and plants, especially in studies of relationships between RNA virus genomic structures and their infectivity. However, we located the m6A site only in a genomic region of approximately 200 bp, which was not sufficiently accurate and precise. Further studies are needed to identify the single nucleotide that was methylated. With the recent development of nanopore sequencing technology [32, 55], accurate and fine genome-wide m6A distribution mapping should be possible.

Most common consensus motifs assay

Earlier biochemical studies, high-throughput transcriptome-based m6A-IP sequencing, and recent nanopore sequencing of animal, viral, and plant mRNAs have shown that the consensus sequence RR[m6A]CH (R=A/G, H=A/C/U) is the most significantly enriched motif in m6A peaks in all eukaryotes to date [5, 15, 32, 56,57,58,59]. In agreement with these reports, we identified 114,546, 104,757, and 117,381 canonical RRACH motifs in mock-, RBSDV-, and RSV-infected samples, respectively. These accounted for 67.9%, 68.9%, and 67.8% of the respective total m6A peaks obtained from sequencing data (Additional file 2: Table S10). Additionally, we analysed the four other most common consensus motifs in rice (“URUAY”, “RAGRAG”, and “UGUAMM”, R=A/G, H=A/C/U, W=A/U, M=A/C) [14, 54] (Additional file 2: Table S10). These motifs also were significantly enriched in the three aforementioned treatments, implying that our sequencing data were confident and believable.

In the enriched m6A peaks (FC > 1.5, p < 0.05), we obtained the most common motifs using de novo prediction software for mock-, RBSDV-, and RSV-infected rice samples (Fig. 7). The “CGXCGXC” (X=A/U/C/G) motif was the most enriched motif in the significantly enriched m6A peaks in mock- and RBSDV-infected rice plants. The “CGXCGXCG” (X=A/U/C/G) motif was enriched in RSV-infected rice. These de novo predicted motifs may indicate potential consensus around the methylated adenine site that is common in the absence and presence of virus infection. Searches for other functional common consensus motifs for m6A modifications need to be performed.

Activation of m6A modification in rice under biotic stress

The potential functions of m6A epigenetic alterations in various eukaryotes and viruses, which are involved in RNA metabolism [29], plant development regulation [24], and RNA stability regulation [51], have been elucidated. However, the detailed and precise roles of m6A are obscure. Presently, the m6A modification process was activated upon virus infection in rice, as evidence by the viral infection associated increase of m6A peak number (Fig. 4A), the differential m6A peak number (Fig. 4C and Additional file 2: Table S6), and the m6A content (Fig. 8A, B). These results revealed the positive association of m6A with plant virus infection.

Viral infection in rice plants is a biotic stress. The m6A modification status during this stress was presently clarified for RSV and RBSDV. We investigated the relative expression of two candidate genes (OsAGO18 and OsSLRL1) and validated their m6A regions on the gene body using m6A-IP-qPCR with corresponding primers (Fig. 8C–F). The m6A modifications were strongly associated with the up- (Fig. 8C, E) or downregulation (Fig. 8D, F) of the target genes, which determined the direct correlation between m6A modification and relative expression. However, compared with the decreased m6A levels evident during TMV infection of tobacco [36], we obtained evidence using different examinations of the opposite result. The dichotomy may have resulted from the different virus infections in different plants. Viral characteristics are very distinct from each other, and their lifestyle in different hosts are often varied [60]. Therefore, different viruses have different effects on the host’s m6A modifications.

Regulation of the relative expression levels of target genes

Key discoveries of mRNA m6A modification functions have mostly concentrated on development [12, 24], such as embryo formation, growth, organ definition [61], apical dominance, trichome branches, gravitropic response, development of lateral roots and vasculature [18, 62], stem cell differentiation in the shoot apical meristem [17], maintenance of circadian and seasonal plant rhythms [32], and flowering [20], all these processes occurred due to the temporal and spatial expression of the key genes in development. Hence, how does m6A methylation regulate the expression of these key genes is crucial in its molecular functional dissections.

In most cases, the post-transcriptional gene regulations are considered to result from host small RNA-based regulation, which mediate the rapid change of relative expression of key pathway genes for better adaptation under stress [63]. Virus infection in rice plants is an important biological stress that alters (in an as yet unclear manner) the expression profiles of genes responsible for defence against the viral pathogen. However, the main group of genes that are m6A modified are unknown, as are the m6A modifications that participate in the response to the infection. The detailed mode of regulation of m6A modifications to the target mRNA are also unclear. Presently, we first focused on the statistical data of the relationship between the gene expression levels and m6A modification. The results showed that m6A mostly occurred in genes that were expressed at low levels during virus infection in plants (Fig. 5), and no other reports about this information in plants under biotic stress. Under abiotic stress conditions, the stability of salt-tolerant transcripts was enhanced by adding the m6A sites for tolerance improvement in salt-stressed Arabidopsis [51]. The WRITER (ALKBH10 family) and READER (ECT2 family) genes were increased, whereas the overall levels of m6A modifications were decreased in maize under drought stress [15]. The READER ECT2 was relocated to the stress granules in heat-stressed cells [24, 64]. These results indicate that the responsive m6A modifications are sensitive and complex under abiotic stress. Similarly, the present findings also indicate sensitivity and complexity of m6A in response to biotic stress, especially in genes expressed at low levels.

More specifically, several targeted pathways that were strongly associated with virus infections were further explored. We analysed the relative expression levels of the m6A modification machinery-related genes (Fig. 9A–D), the main antiviral RNA silencing pathway-related genes (Fig. 10A–C), potential resistance genes (Fig. 10D), and seven phytohormone metabolic pathway-related genes (Additional file 1: Fig. S3–S9) using qRT-PCR. Combined with m6A-IP-sequencing data, the m6A modifications were strongly associated with these pathways (Figs. 9E and 10E, Additional file 2: Table S14). In addition, we found that m6A modification was also involved in the top five KEGG pathway-related genes that originated from the RNA-seq data (Additional file 2: Table S9, S11, S12, S16, and S17). The collective results indicate that m6A actively participates in interactions between plants and viruses and has essential roles in the regulation of the relative expression of target genes.

With respect to the main components of the m6A modification machinery, we found that six component genes of the m6A machinery for RBSDV and RSV could be modified by m6A (Fig. 9E). Their relative expression levels also varied (Fig. 9A–D). These results suggest that the m6A machinery could regulate the expression of main component genes through m6A methylation, which influences the overall cellular m6A levels of the host under conditions of viral infection. In previous studies, the WRITER genes were downregulated, whereas the ERASER genes were down-regulated in TMV infection of tobacco [36]. The findings similarly showed the regulatory roles of virus infection in the main m6A machinery components. Suppression of AtALKBH9B expression, the systemic infection of Arabidopsis with AMV was impaired, while CMV infection did not influence [23]. Together with our findings, these results imply the potential of m6A modification in regulating plant antiviral defence processes.

As the main plant antiviral strategies, RNA silencing pathways and resistance genes play significant roles in countering virus infections [65, 66]. Presently, we first analysed the relative expression levels of nine OsDCLs genes, five OsRDRs genes, 17 OsAGOs genes, and nine reported resistance genes in rice infected with RBSDV or RSV (Fig. 10A–D). Based on the m6A-IP-sequencing date, we mapped the detailed m6A modified genes involved in the RNA silencing pathways (Fig. 10E). Surprisingly, many essential genes, including OsAGO18 [40] and OsAGO2 [67], could be m6A modified, and their relative expression levels were also affected. In addition, other RNA silencing pathway-related genes were also m6A modified during infection with RBSDV or RSV. Except for the small RNA-based gene expression regulatory strategies in plants, our results provide the first evidence that m6A could act as an important regulator of gene expression to participate in the regulation of the main antiviral RNA silencing pathway-related genes.

Phytohormones are the main chemical messengers in the regulation of the physiological processes from germination to senescence during the whole life cycle of plants. Phytohormones have significant roles in defending against biotic and abiotic stresses in plants, especially against viral infections [49, 68]. RBSDV and RSV are two important and devastating rice pathogens. These virus infections occur at epidemic proportions and severely affect the quality and yield of rice in the main production area of East Asia [69]. In recent years, that different types of phytohormones balancing plant growth and antiviral defence processes in plants have been demonstrated. They are JA [40, 70,71,72,73,74,75], SA [45, 71, 72], ABA [71], auxin [72], CTK [71], ET [70], and BR [40, 70, 74, 76], respectively. All these results imply that these seven phytohormones widely participate in the regulation of the physiological status of the host plant to coordinate response to virus infections. Similarly, we wanted to determine whether m6A modifications are also involved in regulating the metabolism of these seven phytohormones. First, the relative expression of these seven phytohormone metabolism pathway-related genes were determined by qRT-PCR (Additional file 1: Fig. S3–S9). Furthermore, m6A-IP sequencing identified the genes that were m6A modified (Additional file 2: Table S14). The data demonstrated that the m6A methylated genes were distributed among nearly all analysed specific phytohormone metabolism pathways during virus infection. These results imply the frequent and important roles of m6A modification in interactions between plants and viruses, especially in the use of phytohormones in defence. Numerous reports have described the relationships between specific phytohormones and virus infection. The present description of the m6A modification-mediated gene expression regulatory strategy in phytohormone metabolism is the first time.

Besides the analysed host m6A machinery, RNA silencing pathways, and seven phytohormone metabolism pathways, we further explored the top five KEGG pathways obtained by differentially expressed genes in RBSDV- and RSV-infected rice plants compared with mock treatment (Additional file 2: Table S11, S16, and S17). KEGG pathway analyses showed that phenylpropanoid biosynthesis (ko00940), plant hormone signal transduction (ko04075), amino sugar and nucleotide sugar metabolism (ko00520), starch and sucrose metabolism (ko00500), and cysteine and methionine metabolism (ko00270) pathways were differentially regulated under RBSDV and RSV infection of rice, with a total of 21 common genes methylated (Additional file 2: Table S17). Combined with the results of the pathways that correlated with virus infection, the collective findings of these analyses once again demonstrated that the m6A modification participates in almost every aspect of host biological process, especially the host’s responsive signal transduction of viral infection.

Taken together, the post-transcriptional mRNA m6A modifications, which were often unrecognised in previous gene expression analyses, play important roles in interactions between host and virus. The results imply that the general or possibly regulatory machinery of the up- or down-regulation of target genes in certain conditions could involve in m6A modifications, in addition to small RNA-based transcriptional regulation, RNA decay, and DNA transcriptional suppression.

Regulation-mode of the relative expression of target genes

Mediation of mRNA physiological status by m6A modification leads to fate determination and has been correlated with a wide range of mRNA metabolism, including changes in the secondary structure [26, 77,78,79], nuclear export [30, 80], pre-miRNA splicing and processing [26, 29, 78, 81,82,83,84], alternative polyadenylation site (APA) choice [85,86,87], RNA stability [17, 24, 88], and RNA maturation [14, 18, 61, 88]. The relationship between m6A modifications and mRNA translation has been intensively studied in recent years [2, 13]. Several groups demonstrated the stimulatory effects of m6A on translation [89,90,91,92,93], while others groups showed inhibitory effects [94,95,96]. The dichotomous effects of m6A on translation indicate the involvement of other, as yet unclear, factors. In terms of detail sites, m6A modification mainly affects the secondary structures of target RNA due to the addition of the methyl group to position 6 of adenylate at special gene body, and these chemical structures’ changes lead to other biological or activity alterations. Understanding the impact of the m6A modification on the relative expression levels of transcripts was the basis to uncover the mechanism of m6A translation.

To determine some regular pattern of m6A modifications to the transcripts’ relative expression levels, we first mapped the m6A modification sites on the gene body and then obtained data of the regulation pattern by using the m6A positions and their corresponding expression profiles from m6A-IP-seq and RNA-seq/qRT-PCR, respectively (Additional file 1: Fig. S10). Seven phytohormone metabolic pathway-related genes comprising a total of 154 genes, were used for m6A and relative expression profile determination analyses (Fig. 10A–D). m6A modification could lead to down-regulation of the target genes in most cases, especially when the m6A modification site was located in the 5′-UTR or CDS/3′-UTR region of the target gene. Previous studies, especially in animal systems, showed that higher levels of m6A modification in the last exon may affect the usage of APA [85], while Bartosovic et al. confirmed that the length of the 3′-UTR could be regulated by the regulation of APA in similar situations [86]. In plants, the rice mutant with disabled FIP37 (WRITER) resulted in chimeric mRNA formation of the two spatially adjacent genes (pair of AT4g30570/580 or AT1g71330/340) [17, 97]. These results confirmed the potential functions of m6A in APA alternatives, which directly determine the translational levels of the target genes. Our finding of the m6A location of the CDS/3′-UTR region may support the reported APA alternative mechanism from the perspective of transcription, and another finding of m6A location of the 5′-UTR regions, which supports that 5′-UTR m6A locations result in higher translational efficiency [13].

To ascertain the regulation pattern between positions of m6A modification sites and the target genes’ expression profile, we integrated the m6A-IP-seq data. We selected the common differentially regulated m6A peaks of both RBSDV- and RSV-infected rice samples that were compared with the mock treatments. Then, we annotated and mapped these differently regulated peaks to the corresponding target genes and obtained 2913 and 2624 genes in RBSDV- and RSV-infected rice samples, respectively (Additional file 2: Table S15). As shown in Fig. 10E, F, wherever the m6A modification sites were located, the target gene transcripts in most cases were downregulated (Fig. 10F) or unchanged (Fig. 10E), especially in the RBSDV-infected samples, consistent with our findings from m6A modifications in seven phytohormone metabolic pathway-related genes (Fig. 10A–D). The RSV-infected samples displayed smaller changes in gene expression under conditions of m6A modifications in any position of the target genes (Fig. 10E). All these results imply that m6A could directly regulate the secondary structures of mRNA by affecting APA usage and alternations. The results also indicate that infections by different viruses could coordinate and regulate different physiological states of plants for the best viral adaptation to host defence, to achieve a balance between plant defence and viral infection.

The present study is the first description of genome-wide m6A distribution maps of rice and of two significant viruses. We further analysed the common consensus motifs and validated the m6A modification in rice. m6A was activated in RBSDV- and RSV-infected rice, and the activation was predominantly in genes that were not actively expressed. Finally, m6A was strongly associated with and directly regulated the expression of the main antiviral pathway-related genes in rice plants. All these findings deepen our understanding of the significant roles of m6A in interactions between plants and viruses, especially the model plant (rice) and two significant model plant viruses (RBSDV and RSV). The findings highlight the need for more profound and detailed investigations of the functions of m6A on select target genes.

Method and materials

Plant materials

The rice japonica cultivar Huaidao 5 (HD-5) was planted in a bucket and cultured in a growth chamber under short-day (SD) conditions (11 h light/13 h dark) with a light intensity of 900 μmol/m2 s1 at 27 °C. HD-5 was highly susceptible to rice stripe virus or rice black-stripe dwarf virus in the field [98].

Artificial inoculation of RBSDV and RSV

Artificial inoculation of RBSDV and RSV was performed as previously described with minor modifications [99, 100]. For RSV inoculation, small brown planthoppers (SBPHs) that carring RSV, which was maintained at the Department of Plant Protection, Yangzhou University, was transferred to rice seedlings cultured in a glass beaker (1 L). Then, the seedlings were divided and transplanted to a 20 L plastic bucket after 3 days’ aspiration. For RBSDV inoculation, RBSDV-infected rice plants were collected from the field and exposed to non-virulent SBPHs for 4 days. As the SBPHs fed, they acquired the RBSDV. The viruliferous SBPH were transferred to a beaker that was cultured with rice seedlings. Three days later, the rice seedlings were transplanted to the same-size plastic bucket mentioned above.

For mock treatment, non-virulent SBPHs were transferred to a glass breaker and cultured with rice seedlings for three days. The seedlings were transplanted in a similar manner as described above for the viruliferous SBPH inoculated seedlings. The RSV- and RBSDV-infected rice plants were equipped with pest control nets and cultured in the growth chamber at Yangzhou University. Total protein was extracted from the planted rice seedlings and analysed. western blot was performed to detect the presence of RSV and RBSDV as described next.

Western blot detection of RSV and RBSDV

Western blot detection was performed as described previously, with minor modifications [101]. Briefly, total proteins of the mock- or virus-inoculated rice seedlings were extracted and dissolved in 2 × SDS protein loading buffer containing 5% β-mercaptoethanol. After 5 min of heating in boiling water and centrifugation at 13,000 g for 10 min, the samples were resolved by 12.5% SDS-PAGE. The total proteins were transferred to a piece of nitrocellulose membrane. Polyclonal antiserum against RSV NS3 and RBSDV p10 protein was used as the first antibody. This was followed by a secondary alkaline phosphatase-linked antibody (AP-A; Sigma-Aldrich, USA). The specific bands were visualised with a colour reaction by adding nitro blue tetrazolium/5-bromo-4-chloro-3-indolyl phosphate (NBT/BCIP) as the substrates (Sangon Biotech, China).

Total RNA exaction and reverse transcription-polymerase chain reaction (RT-PCR) detection

Trizol Reagent (TaKaRa Bio, Japan) was used for rice total RNA extraction as described previously [102]. RNase-free rDNase I (TaKaRa Bio, Japan) was added to remove contaminating genomic DNA. Reverse transcription reactions were performed as described previously with RSV- and RBSDV-specific primers, respectively (Additional file 2: Table S1) [101]. Synthetic cDNA was used as templates, and PCR reactions were performed with specific pairs of primers as described.

High-throughput RNA m6A-IP-seq and RNA-seq

Mock, RSV-, and RBSDV-infected rice plants were used for total RNA isolation. Ribosome RNA (rRNA) was captured and deleted using the RiboMinus™ Plant Kit for RNA-seq (Thermo Fisher Scientific, USA). RNA Fragmentation Reagents (Ambion, USA) were added and the RNA was randomly fragmented into approximately 200 nucleotide fragments. The fragmented RNA was divided into two identical portions. One was used in the m6A-IP experiment, the other was used as the input for the next RNA-seq.

Anti-m6A polyclonal antibody (#ab208577; Abcam, USA) was added to the obtained fragmented RNA and incubated at 4°C for 2 h in immunoprecipitation (IP) buffer (10 mM Tris-HCl, pH 7.4, 150 mM NaCl, 0.1% [v/v] octylphenoxypolyethoxyethanol [Igepal CA-630]) with the addition of RNase Plus RNase inhibitor (Promega, USA). The mixture was immunoprecipitated by incubation with protein-A conjugated beads at 4 °C for 2 h. After adequate washing and removal of non-specifically bound compounds, 0.8 mg/mL N6-methyladenosine in IP buffer was used to elute the bound RNA. TruSeq standard mRNA Sample Prep Kit (Illumina, USA) was used to construct the library of immunoprecipitated RNA and input RNA as described previously [15]. High-throughput sequencing was performed on the Hiseq X Ten platform (Illumina) with the libraries corresponding to immunoprecipitated RNA. Single-end reads of 50-bp-length were obtained. The libraries from the input RNA were also sequenced and 122-bp-length double-end reads were produced.

Analysis of the sequence data

Row reads from m6A-seq and RNA-seq were trimmed by adaptor sequence removal, followed by treatment with the Trimmomatic v0.36 tool [103] to remove the low-quality bases. The FastQC programme was used to examine the quality of the trimmed m6A-seq and RNA-seq clear reads (https://www.bioinformatics.babraham.ac.uk/projects/fastqc).

Distribution of m6A sites in rice and viral genomes

The obtained RNA-seq clear reads of mock-treated, RSV-, and RBSDV-infected rice were aligned to the rice (Oryza sativa. IRGSP-1.0), RBSDV, and RSV reference genome, respectively, using Tophat v2.1.1 [104]. The maximum intron length was set to 10 kb with default settings for other parameters. Cufflinks v2.2.1 [105] was used to analyse the unique mapping reads. The term FPKM mapped reads (FPKM = Counts of mapped fragments × 109/[length of transcript × total count of the mapped fragments]) was used for normalisation and estimation of the gene expression level. The Cuffdiff programme in Cufflinks v2.2.1 was used for differential analysis. The O. sativa reference genome sequences and annotation were download from Ensembl Plants (release 45; https://plants.ensembl.org)[106]. The reference genome sequences of RSV and RBSDV were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov/).

Peak calling analysis

The obtained m6A-seq clear reads of mock-treated, RSV-, and RBSDV-infected rice samples were aligned to their respective reference genome as described above. The m6A enriched peaks in each m6A-IP sample were determined by MeTDiff peak calling software [107] with the corresponding input sample serving as control. The running options of the MeTDiff software were set as: FRAGMENT_LENGTH = 200, PEAK_CUTOFF_PVALUE = 0.01, PEAKS_CUTOFF_FDR = 0.05. After the running, the m6A peak was detected, and the called peaks were annotated by intersection with gene architecture by ChIPseeker [108].

Identification of enriched motifs within m6A peaks

The MEME (http://meme-suite.org/tools/dreme) [109] and Discriminative Regular Expression Motif Elicitation (DREME) [110] tools suites located in the MEME suite were used to perform motif enhancement analyses. The DREME suite was used to discover short (≤ 8 bp) and ungapped motifs enriched within series target m6A peak sequences relative to control sequences (shuffled m6A peak sequences). The series target m6A peak sequences were extracted from the rice reference genome (O. sativa. IRGSP-1.0) using the fastaFromBed function in BEDTools software v2.28 [111]. The control sequences were randomly produced by shuffling each m6A peak sequence while retaining the nucleotide frequencies. The performance was executed by “fasta-shuffle-letters” utility (k = 1) in the MEME suite. These two aforementioned types of sequences were loaded into DREME for discovery motifs. The parameter settings were as follows: meme-minw, 8; meme-maxw 30; meme-nmotifs, 10; E dreme e-value threshold, 0.05; Fimo-skip-spamo-skip-centrimo-score, 5.0; centrimo-ethresh, 10.0; The AME tools (http://meme-suite.org/tools/ame) package in the MEME suite was used to calculate the significance level of relative enrichment of one special motif within target sequences relative to the control sequences.

Gene expression quantification

Total RNA was extracted as described previously [102]. After treatment with RNase-free rDNase I (TaKaRa Bio, Japan), the first-strand cDNA was produced using PrimeScript II First Strand cDNA Synthesis Kit (TaKaRa Bio, Japan) according to the manufacturer’s instructions. Quantitative PCR (qPCR) was performed using the CFX96 Real-Time PCR Detection System with SsoFast EvaGreen Supermix (Bio-Rad, USA). The quantitative experiment was performed with three biological replicates and three technical replicates. Rice eEF1α was chosen as the internal reference gene [112]. The gene expression levels were analysed using the Bio-Rad CFX Manager software. The pairs of primers used for gene expression quantification are listed in Additional file 2: Tables S18 and S19.

Total mRNA isolation and Dot-blot analyses

The extracted total RNAs of the rice plant by Trizol Reagent (TaKaRa Bio, Japan). Then, the Ambion® Poly(A)Purist™ MAG kit (InvitrogenTM, USA) was taken for total mRNA isolation according to the manufacturer's instructions. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Dot blotting was performed with minor modification as described previously [15]. Briefly, the extracted total mRNA was treated with DNase and quantified with a NanoDrop™ 1000 Spectrophotometer (Thermo Fisher Scientific, USA). The loading volume to each dot was calculated according to the concentration of the purified total RNA. The RNAs were denatured by heating at 75 °C for 5 min, followed by chilling on ice for at least 10 min before loading. The samples were spotted on a Hybond-N+ membrane (Amersham, USA), followed by drying at 25 °C for 20 min. After ultraviolet crosslinking, the membrane was blocked with 3% skimmed milk in 1 × PBST buffer (PBS containing 0.05% Tween-20) for 3 h. Then, the membrane was washed three times for 10 min interval with 1× PBST buffer. The membrane was incubated with anti-m6A antibody (#ab208577; Abcam) in 1× PBST buffer for 1 h at 4 °C. After washing three times, the membrane was incubated with alkaline phosphatase (AP)-conjugated anti-mouse immunoglobulin G secondary antibody (#ab98712; Abcam). The colour reaction was performed as in the western blot.

Liquid chromatography-tandem mass spectrometry (LC-MS/MS) analysis of m6A contents

RNA samples isolated by the kit above were treated with a digestion cocktail containing nuclease P1 (Sigma-Aldrich, USA) and phosphodiesterase I (Sigma-Aldrich). The digestion was performed in buffer (10 mM Tris-HCl pH 6.8, 30 mM sodium acetate, 2 mM ZnCl2) with 2 U and 0.01 U of each corresponding enzyme at 37 °C for 12 h, followed by further digestion with 2 U AP for 1 h at 37 °C. The digestion enzymes were removed using a 10 kDa filtration tube (Amicon Ultra, USA). RNA solutions were diluted for 10 times, and 10 μL of the dilution was used for LC-MS/MS. Reverse-phase high-performance liquid chromatography (Agilent, USA) was used to separate the nucleosides using an Agilent C 18 column (5 μm size, 150 mm × 2.1 mm), which was combined with MS detection using a QTRAP 5500 device (AB Sciex, USA). The detailed parameters were set as previously described [113].

Validation of m6A site using m6A-IP-qPCR

Total RNA was extracted from mock-, RSV-, and RBSDV-infected samples at different times (0, 1, 2, 4, 8, and 16 days post-inoculation). After rRNA removal by the Ambion® Poly(A)Purist™ MAG kit, RNA fragmentation, and m6A-IP experiments were performed as described above [16]. The precipitated RNAs from N6-methyladenosine eluted solution were used as templates. The fold enrichment of each fragment was determined by qRT-PCR. The designed specific primers are shown in Additional file 2: Table S18.

Availability of data and materials

All the datasets supporting the conclusions of this article are available in the Additional file 2: Tables included in this manuscript. The raw sequencing data were uploaded and deposited in NCBI (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA732795) [114]. The accession number was from SRR14675282 to SRR14675293. The details were showed as below: RNA-sequencing data were SRR14675288 (Mock_Rice_1), SRR14675292 (Mock_Rice_2), SRR14675284 (RBSDV_Rice_1), SRR14675286 (RBSDV_Rice_2), SRR14675282 (RSV_Rice_1), and SRR14675290 (RSV_Rice_2), respectively. m6A-IP-sequencing data were SRR14675289 (Mock_Rice_1), SRR14675293 (Mock_Rice_2), SRR14675285 (RBSDV_Rice_1), SRR14675287 (RBSDV_Rice_2), SRR14675283 (RSV_Rice_1), and SRR14675291 (RSV_Rice_2), respectively. The uncropped western blot images were shown in Additional file 1: Fig. S12.

References

  1. 1.

    Vanyushin BF, Tkacheva SG, Belozersky AN. Rare bases in animal DNA. Nature. 1970;225(5236):948–9. https://doi.org/10.1038/225948a0.

    CAS  Article  PubMed  Google Scholar 

  2. 2.

    Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3' UTRs and near stop codons. Cell. 2012;149(7):1635–46. https://doi.org/10.1016/j.cell.2012.05.003.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  3. 3.

    Zheng G, Dahl JA, Niu Y, Fedorcsak P, Huang CM, Li CJ, et al. ALKBH5 is a mammalian RNA demethylase that impacts RNA metabolism and mouse fertility. Mol Cell. 2012;49:18–29.

  4. 4.

    Liu J, Dou X, Chen C, Chen C, Liu C, Xu MM, et al. N6-methyladenosine of chromosome-associated regulatory RNA regulates chromatin state and transcription. Science. 2020;367(6477):580–6. https://doi.org/10.1126/science.aay6018.

  5. 5.

    Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. 2012;485(7397):201–6. https://doi.org/10.1038/nature11112.

  6. 6.

    Lichinchi G, Gao S, Saletore Y, Gonzalez GM, Bansal V, Wang Y, et al. Dynamics of the human and viral m6A RNA methylomes during HIV-1 infection of T cells. Nature Microbiol. 2016;1(4):16011. https://doi.org/10.1038/nmicrobiol.2016.11.

  7. 7.

    Perry RP, Kelley DE, Friderici K, Rottman F. The methylated constituents of L cell messenger RNA: Evidence for an unusual cluster at the 5′ terminus. Cell. 1975;4(4):387–94. https://doi.org/10.1016/0092-8674(75)90159-2.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Desrosiers RC, Friderici KH, Rottman FM. Characterization of Novikoff hepatoma mRNA methylation and heterogeneity in the methylated 5' terminus. Biochemistry. 1975;14(20):4367–74. https://doi.org/10.1021/bi00691a004.

    CAS  Article  PubMed  Google Scholar 

  9. 9.

    Mondo SJ, Dannebaum RO, Kuo RC, Louie KB, Bewick AJ, LaButti K, et al. Widespread adenine N6-methylation of active genes in fungi. Nat Genet. 2017;49(6):964–8. https://doi.org/10.1038/ng.3859.

  10. 10.

    Chao Z, Changshi W, Hongbo L, Qiangwei Z, Qian L, Yan G, et al. Identification and analysis of adenine N6-methylation sites in the rice genome. Nature Plants. 2018;4:554–63.

  11. 11.

    Gianluigi, Lichinchi, Boxuan Simen, Zhao, Yinga, Wu, Zhike, Lu, Yue, Qin: Dynamics of human and viral RNA methylation during zika virus infection. Cell Host Microbe 2016, 20:666-673.

  12. 12.

    Arribas-Hernández L, Bressendorff S, Hansen MH, Poulsen C, Erdmann S, Brodersen P. An m6A-YTH module controls developmental timing and morphogenesis in Arabidopsis. Plant Cell. 2018;30(5):952–67. https://doi.org/10.1105/tpc.17.00833.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Luo JH, Wang Y, Wang M, Zhang LY, Peng HR, Zhou YY, et al. Natural variation in RNA m6A methylation and its relationship with translational status. Plant Physiol. 2019;182:332–44.

  14. 14.

    Zhang F, Zhang YC, Liao JY, Yu Y, Chen YQ. The subunit of RNA N6-methyladenosine methyltransferase OsFIP regulates early degeneration of microspores in rice. Plos Genet. 2019;15(5):e1008120. https://doi.org/10.1371/journal.pgen.1008120.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Miao Z, Zhang T, Qi Y, Song J, Han Z, Ma C. Evolution of the RNA N6-methyladenosine methylome mediated by genomic duplication. Plant Physiol. 2020;182(1):345–60. https://doi.org/10.1104/pp.19.00323.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Zhang Q, Liang Z, Cui X, Ji C, Li Y, Zhang P, et al. N6-Methyladenine DNA methylation in japonica and indica rice genomes and its association with gene expression, plant development, and stress responses. Mol Plant. 2018;11(12):1492–508. https://doi.org/10.1016/j.molp.2018.11.005.

  17. 17.

    Shen L, Liang Z, Gu X, Chen Y, Teo ZWN, Hou X, Cai WM, Dedon PC, Liu L, Yu H. N6-Methyladenosine RNA modification regulates shoot stem cell fate in Arabidopsis. Dev Cell. 2016;38(2):186–200. https://doi.org/10.1016/j.devcel.2016.06.008.

  18. 18.

    Ruzicka K, Zhang M, Campilho A, Bodi Z, Kashif M, Saleh M, et al. Identification of factors required for m6A mRNA methylation in Arabidopsis reveals a role for the conserved E3 ubiquitin ligase HAKAI. New Phytol. 2017;215(1):157–72. https://doi.org/10.1111/nph.14586.

  19. 19.

    Liang Z, Geng Y, Gu X. Adenine methylation: new epigenetic marker of DNA and mRNA. Molecular Plant. 2018;11(10):1219–21. https://doi.org/10.1016/j.molp.2018.08.001.

    CAS  Article  PubMed  Google Scholar 

  20. 20.

    Duan HC, Wei LH, Zhang C, Wang Y, Chen L, Lu Z, et al. ALKBH10B is an RNA N6-methyladenosine demethylase affecting Arabidopsis floral transition. Plant Cell. 2017;29(12):2995–3011. https://doi.org/10.1105/tpc.16.00912.

  21. 21.

    Cantoni GL. Biological methylation: selected aspects. Annu Rev Biochem. 1975;44(1):435–51. https://doi.org/10.1146/annurev.bi.44.070175.002251.

    CAS  Article  PubMed  Google Scholar 

  22. 22.

    Sekula B, Ruszkowski M, Dauter Z. S-adenosylmethionine synthases in plants: Structural characterization of type I and II isoenzymes from Arabidopsis thaliana and Medicago truncatula. Int J Biol Macromol. 2020;151:554–65. https://doi.org/10.1016/j.ijbiomac.2020.02.100.

    CAS  Article  PubMed  Google Scholar 

  23. 23.

    M M-P, Aparicio F, MP L-G, JM B, JA S-N, V P: Arabidopsis m6A demethylase activity modulates viral infection of a plant virus and the m6A abundance in its genomic RNAs. Proc Natl Acad Sci U S A 2017, 114:10755.

  24. 24.

    Wei LH, Song P, Wang Y, Lu Z, Tang Q, Yu Q, et al. The m6A reader ECT2 controls trichome morphology by affecting mRNA stability in Arabidopsis. Plant Cell. 2018;30:tpc.00934.02017.

  25. 25.

    Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature. 2014;505(7481):117–20. https://doi.org/10.1038/nature12730.

  26. 26.

    Liu N, Dai Q, Zheng G, He C, Parisien M, Pan T. N6-methyladenosine-dependent RNA structural switches regulate RNA-protein interactions. Nature. 2015;518(7540):560–4. https://doi.org/10.1038/nature14234.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Zhou J, Wan J, Gao X, Zhang X, Jaffrey SR, Qian SB. Dynamic m6A mRNA methylation directs translational control of heat shock response. Nature. 2015;526(7574):591–4. https://doi.org/10.1038/nature15377.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Fustin JM, Doi M, Yamaguchi Y, Hida H, Nishimura S, Yoshida M, et al. RNA-methylation-dependent RNA processing controls the speed of the circadian clock. Cell. 2013;155(4):793–806. https://doi.org/10.1016/j.cell.2013.10.026.

  29. 29.

    Bhat SS, Bielewicz D, Gulanicz T, Bodi Z, Yu X, Anderson SJ, et al. mRNA adenosine methylase (MTA) deposits m6A on pri-miRNAs to modulate miRNA biogenesis in Arabidopsis thaliana. Proc Natl Acad Sci U S A. 2020;117(35):21785–95. https://doi.org/10.1073/pnas.2003733117.

  30. 30.

    Edens BM, Vissers C, Su J, Arumugam S, Xu Z, Shi H, et al. FMRP modulates neural differentiation through m6A-dependent mRNA nuclear export. Cell Rep. 2019;28:845–854.e845.

  31. 31.

    Arribas-Hernández L, Brodersen P. Occurrence and functions of m6A and other covalent modifications in plant mRNA. Plant Physiol. 2020;182(1):79–96. https://doi.org/10.1104/pp.19.01156.

    CAS  Article  PubMed  Google Scholar 

  32. 32.

    Parker MT, Knop K, Sherwood AV, Schurch NJ, Mackinnon K, Gould PD, et al. Nanopore direct RNA sequencing maps the complexity of Arabidopsis mRNA processing and m6A modification. eLife. 2020;9:e49658. https://doi.org/10.7554/eLife.49658.

  33. 33.

    Edward M, Kennedy, Hal P, Bogerd, Anand VR, Kornepati, Dong, Kang, Delta, Ghoshal: Posttranscriptional m6A editing of HIV-1 mRNAs enhances viral gene expression. Cell Host Microbe 2016, 19:675-685.

  34. 34.

    Tirumuru N, Zhao BS, Lu W, Lu Z, Wu L. Correction: N6-methyladenosine of HIV-1 RNA regulates viral infection and HIV-1 Gag protein expression. eLife. 2017;6. https://doi.org/10.7554/eLife.31482.

  35. 35.

    Gokhale N, Mcintyre AR, Mcfadden M, Roder A, Horner S. N6-Methyladenosine in flaviviridae viral RNA genomes regulates infection. Cell Host Microbe. 2016;20(5):654–65. https://doi.org/10.1016/j.chom.2016.09.015.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Li Z, Shi J, Yu L, Zhao X, Ran L, Hu D, et al. N6-methyl-adenosine level in Nicotiana tabacum is associated with tobacco mosaic virus. Virol J. 2018;15(1):87. https://doi.org/10.1186/s12985-018-0997-4.

  37. 37.

    Bratlie MS, Drablos F. Bioinformatic mapping of AlkB homology domains in viruses. BMC Genomics. 2005;6(1):1. https://doi.org/10.1186/1471-2164-6-1.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  38. 38.

    van den Born E, Omelchenko MV, Bekkelund A, Leihne V, Koonin EV, Dolja VV, et al. Viral AlkB proteins repair RNA damage by oxidative demethylation. Nucleic Acids Res. 2008;36(17):5451–61. https://doi.org/10.1093/nar/gkn519.

  39. 39.

    Wu J, Yang Z, Wang Y, Zheng L, Ye R, Ji Y, et al. Viral-inducible Argonaute18 confers broad-spectrum virus resistance in rice by sequestering a host microRNA. eLife. 2015;4:e05733. https://doi.org/10.7554/eLife.05733.

  40. 40.

    Yang Z, Huang Y, Yang J, Yao S, Zhao K, Wang D, et al. Jasmonate signaling enhances RNA silencing and antiviral defense in rice. Cell Host Microbe. 2020;28:89–103.e108.

  41. 41.

    Um TY, Lee HY, Lee S, Chang SH, Chung PJ, Oh KB, et al. Jasmonate Zim-domain protein 9 interacts with slender rice 1 to mediate the antagonistic interaction between jasmonic and gibberellic acid signals in rice. Front Plant Sci. 2018;9:1866. https://doi.org/10.3389/fpls.2018.01866.

  42. 42.

    De Vleesschauwer D, Seifi HS, Filipe O, Haeck A, Huu SN, Demeestere K, et al. The DELLA protein SLR1 integrates andamplifies salicylic acid- and jasmonic acid-dependent innate immunity in rice. Plant Physiol. 2016;170(3):1831–47. https://doi.org/10.1104/pp.15.01515.

  43. 43.

    Mo W, Tang W, Du Y, Jing Y, Bu Q, Lin R. PHYTOCHROME-INTERACTING FACTOR-LIKE14 and SLENDER RICE1 interaction controls seedling growth under salt stress. Plant Physiol. 2020;184(1):506–17. https://doi.org/10.1104/pp.20.00024.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Liu Q, Deng S, Liu B, Tao Y, Ai H, Liu J, et al. A helitron-induced RabGDIα variant causes quantitative recessive resistance to maize rough dwarf disease. Nat Commun. 2020;11(1):495. https://doi.org/10.1038/s41467-020-14372-3.

  45. 45.

    Wang Q, Liu Y, He J, Zheng X, Hu J, Liu Y, et al. STV11 encodes a sulphotransferase and confers durable resistance to rice stripe virus. Nature Commun. 2014;5(1):4768. https://doi.org/10.1038/ncomms5768.

  46. 46.

    Hayano-Saito Y, Hayashi K. Stvb-i, a rice gene conferring durable resistance to rice stripe virus, protects plant growth from heat stress. Front Plant Sci. 2020;11:519. https://doi.org/10.3389/fpls.2020.00519.

    Article  PubMed  PubMed Central  Google Scholar 

  47. 47.

    Santner A, Estelle M. Recent advances and emerging trends in plant hormone signalling. Nature. 2009;459(7250):1071–8. https://doi.org/10.1038/nature08122.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Durbak A, Yao H, McSteen P. Hormone signaling in plant development. Curr Opin Plant Biol. 2012;15(1):92–6. https://doi.org/10.1016/j.pbi.2011.12.004.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Collum TD, Culver JN. The impact of phytohormones on virus infection and disease. Curr Opin Virol. 2016;17:25–31. https://doi.org/10.1016/j.coviro.2015.11.003.

    CAS  Article  PubMed  Google Scholar 

  50. 50.

    Kane SE, Beemon K. Precise localization of m6A in rous sarcoma virus RNA reveals clustering of methylation sites: implications for RNA processing. Mol Cell Biol. 1985;5(9):2298–306. https://doi.org/10.1128/mcb.5.9.2298-2306.1985.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Anderson SJ, Kramer MC, Gosai SJ, Yu X, Vandivier LE, Nelson ADL, et al. N6-methyladenosine inhibits local ribonucleolytic cleavage to stabilize mRNAs in Arabidopsis. Cell Rep. 2018;25:1146–1157.e1143.

  52. 52.

    Luo G-Z, MacQueen A, Zheng G, Duan H, Dore LC, Lu Z, et al. Unique features of the m6A methylome in Arabidopsis thaliana. Nat Commun. 2014;5(1):5630. https://doi.org/10.1038/ncomms6630.

  53. 53.

    Wan Y, Tang K, Zhang D, Xie S, Zhu X, Wang Z, et al. Transcriptome-wide high-throughput deep m6A-seq reveals unique differential m6A methylation patterns between three organs in Arabidopsis thaliana. Genome Biology. 2015;16(1):272. https://doi.org/10.1186/s13059-015-0839-2.

  54. 54.

    Li Y, Wang X, Li C, Hu S, Yu J, Song S. Transcriptome-wide N6-methyladenosine profiling of rice callus and leaf reveals the presence of tissue-specific competitors involved in selective mRNA modification. RNA Biology. 2014;11(9):1180–8. https://doi.org/10.4161/rna.36281.

    Article  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Garalde DR, Snell EA, Jachimowicz D, Sipos B, Lloyd JH, Bruce M, et al. Highly parallel direct RNA sequencing on an array of nanopores. Nat Methods. 2018;15(3):201–6. https://doi.org/10.1038/nmeth.4577.

  56. 56.

    Schibler U, Kelley DE, Perry RP. Comparison of methylated sequences in messenger RNA and heterogeneous nuclear RNA from mouse L cells. J Mol Biol. 1977;115(4):695–714. https://doi.org/10.1016/0022-2836(77)90110-3.

    CAS  Article  PubMed  Google Scholar 

  57. 57.

    Canaani D, Kahana C, Lavi S, Groner Y. Identification and mapping of N6-methyladenosine containing sequences in simian virus 40 RNA. Nucleic Acids Res. 1979;6(8):2879–99. https://doi.org/10.1093/nar/6.8.2879.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  58. 58.

    Li B, Wang X, Li Z, Lu C, Zhang Q, Chang L, et al. Transcriptome-wide analysis of N6-methyladenosine uncovers its regulatory role in gene expression in the lepidopteran Bombyx mori. Insect Mol Biol. 2019;28(5):703–15. https://doi.org/10.1111/imb.12584.

  59. 59.

    Zhao BS, Wang X, Beadell AV, Lu Z, Shi H, Kuuspalu A, et al. m6A-dependent maternal mRNA clearance facilitates zebrafish maternal-to-zygotic transition. Nature. 2017;542(7642):475–8. https://doi.org/10.1038/nature21355.

  60. 60.

    King AMQ, Adams MJ, Eric B, Carstens LEJ. Virus taxonomy: Ninth report of the international committee on taxonomy of viruses. San Diego: Elsevier; 2011.

    Google Scholar 

  61. 61.

    Zhong S, Li H, Bodi Z, Button J, Vespa L, Herzog M, et al. MTA is an Arabidopsis messenger RNA adenosine methylase and interacts with a homolog of a sex-specific splicing factor. Plant Cell. 2008;20(5):1278–88. https://doi.org/10.1105/tpc.108.058883.

  62. 62.

    Bodi Z, Zhong S, Mehra S, Song J, Graham N, Li H, et al. Adenosine methylation in Arabidopsis mRNA is associated with the 3′ end and reduced levels cause developmental defects. Front Plant Sci. 2012;3:48.

  63. 63.

    Sunkar R, Chinnusamy V, Zhu J, Zhu JK. Small RNAs as big players in plant abiotic stress responses and nutrient deprivation. Trends Plant Sci. 2007;12(7):301–9. https://doi.org/10.1016/j.tplants.2007.05.001.

    CAS  Article  PubMed  Google Scholar 

  64. 64.

    Scutenaire J, Deragon JM, Jean V, Benhamed M, Raynaud C, Favory JJ, et al. The YTH domain protein ECT2 is an m6A reader required for normal trichome branching in Arabidopsis. Plant Cell. 2018;30:tpc.00854.02017.

  65. 65.

    Jin Y, Zhao J-H, Guo H-S. Recent advances in understanding plant antiviral RNAi and viral suppressors of RNAi. Curr Opin Virol. 2021;46:65–72. https://doi.org/10.1016/j.coviro.2020.12.001.

    CAS  Article  PubMed  Google Scholar 

  66. 66.

    Guo Z, Li Y, Ding SW. Small RNA-based antimicrobial immunity. Nature Rev Immunol. 2019;19(1):31–44. https://doi.org/10.1038/s41577-018-0071-x.

    CAS  Article  Google Scholar 

  67. 67.

    Zheng S, Li J, Ma L, Wang H, Zhou H, Ni E, et al. OsAGO2 controls ROS production and the initiation of tapetal PCD by epigenetically regulating OsHXK1 expression in rice anthers. Proc Natl Acad Sci U S A. 2019;116(15):7549–58. https://doi.org/10.1073/pnas.1817675116.

  68. 68.

    Qin J, Wang C, Wang L, Zhao S, Wu J. Defense and counter-defense in rice–virus interactions. Phytopathol Res. 2019;1(1):34. https://doi.org/10.1186/s42483-019-0041-7.

    Article  Google Scholar 

  69. 69.

    Hibino H. Biology and epidemiology of rice viruses. Annu Rev Phytopathol. 1996;34(1):249–74. https://doi.org/10.1146/annurev.phyto.34.1.249.

    CAS  Article  PubMed  Google Scholar 

  70. 70.

    He Y, Zhang H, Sun Z, Li J, Hong G, Zhu Q, et al. Jasmonic acid-mediated defense suppresses brassinosteroid-mediated susceptibility to rice black streaked dwarf virus infection in rice. New Phytol. 2017;214(1):388–99. https://doi.org/10.1111/nph.14376.

  71. 71.

    Renyan H, Yueyue L, Guilin T, Shugang H, Zeyu Y, Junwei Z, et al. Dynamic phytohormone profiling of rice upon rice black-streaked dwarf virus invasion. J Plant Physiol. 2018;228:92–100.

  72. 72.

    Zhang H, Tan X, Li L, He Y, Hong G, Li J, et al. Suppression of auxin signalling promotes rice susceptibility to rice black streaked dwarf virus infection. Mol Plant Pathol. 2019;20(8):1093–104. https://doi.org/10.1111/mpp.12814.

  73. 73.

    Han K, Huang H, Zheng H, Ji M, Yuan Q, Cui W, et al. Rice stripe virus coat protein induces the accumulation of jasmonic acid, activating plant defence against the virus while also attracting its vector to feed. Mol Plant Pathol. 2020;21(12):1647–53. https://doi.org/10.1111/mpp.12995.

  74. 74.

    Hu J, Huang J, Xu H, Wang Y, Li C, Wen P, et al. Rice stripe virus suppresses jasmonic acid-mediated resistance by hijacking brassinosteroid signaling pathway in rice. Plos Pathog. 2020;16(8):e1008801. https://doi.org/10.1371/journal.ppat.1008801.

  75. 75.

    He L, Chen X, Yang J, Zhang T, Li J, Zhang S, et al. Rice black-streaked dwarf virus-encoded P5-1 regulates the ubiquitination activity of SCF E3 ligases and inhibits jasmonate signaling to benefit its infection in rice. New Phytol. 2020;225(2):896–912. https://doi.org/10.1111/nph.16066.

  76. 76.

    He Y, Hong G, Zhang H, Tan X, Li L, Kong Y, et al. The OsGSK2 kinase integrates brassinosteroid and jasmonic acid signaling by interacting with OsJAZ4. Plant Cell. 2020;32(9):2806–22. https://doi.org/10.1105/tpc.19.00499.

  77. 77.

    Roost C, Lynch SR, Batista PJ, Qu K, Chang HY, Kool ET. Structure and thermodynamics of N6-methyladenosine in RNA: a spring-loaded base modification. J Am Chem Soc. 2015;137(5):2107–15. https://doi.org/10.1021/ja513080v.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Liu N, Zhou KI, Parisien M, Dai Q, Diatchenko L, Pan T. N6-methyladenosine alters RNA structure to regulate binding of a low-complexity protein. Nucleic Acids Res. 2017;45(10):6051–63. https://doi.org/10.1093/nar/gkx141.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  79. 79.

    Spitale RC, Flynn RA, Zhang QC, Crisalli P, Lee B, Jung JW, et al. Structural imprints in vivo decode RNA regulatory mechanisms. Nature. 2015;519(7544):486–90. https://doi.org/10.1038/nature14263.

  80. 80.

    Roundtree IA, Luo GZ, Zhang Z, Wang X, Zhou T, Cui Y, et al. YTHDC1 mediates nuclear export of N6-methyladenosine methylated mRNAs. eLife. 2017;6. https://doi.org/10.7554/eLife.31311.

  81. 81.

    Haussmann IU, Bodi Z, Sanchez-Moran E, Mongan NP, Archer N, Fray RG, et al. m6A potentiates Sxl alternative pre-mRNA splicing for robust Drosophila sex determination. Nature. 2016;540(7632):301–4. https://doi.org/10.1038/nature20577.

  82. 82.

    Lence T, Akhtar J, Bayer M, Schmid K, Spindler L, Ho CH, et al. m6A modulates neuronal functions and sex determination in Drosophila. Nature. 2016;540(7632):242–7. https://doi.org/10.1038/nature20568.

  83. 83.

    Xiao W, Adhikari S, Dahal U, Chen YS, Hao YJ, Sun BF, et al. Nuclear m6A reader YTHDC1 regulates mRNA splicing. Mol Cell. 2016;61(4):507–19. https://doi.org/10.1016/j.molcel.2016.01.012.

  84. 84.

    Pendleton KE, Chen B, Liu K, Hunter OV, Xie Y, Tu BP, et al. The U6 snRNA m6A methyltransferase METTL16 regulates SAM synthetase intron retention. Cell. 2017;169:824–835.e814.

  85. 85.

    Ke S, Alemu EA, Mertens C, Gantman EC, Fak JJ, Mele A, et al. A majority of m6A residues are in the last exons, allowing the potential for 3′ UTR regulation. Genes Dev. 2015;29(19):2037–53. https://doi.org/10.1101/gad.269415.115.

  86. 86.

    Marek, Bartosovic, Helena, Covelo, Molares, Pavlina, Gregorova, Dominika, Hrossova, Grzegorz: N6-methyladenosine demethylase FTO targets pre-mRNAs and regulates alternative splicing and 3'-end processing. Nucleic Acids Res 2017, 45:11356-11370.

  87. 87.

    Yue Y, Liu J, Cui X, Cao J, Luo G, Zhang Z, et al. VIRMA mediates preferential m6A mRNA methylation in 3'UTR and near stop codon and associates with alternative polyadenylation. Cell Discov. 2018;4(1):10. https://doi.org/10.1038/s41421-018-0019-0.

  88. 88.

    Hofmann NR. Epitranscriptomics and flowering: mRNA methylation/demethylation regulates flowering time. Plant Cell. 2017;29(12):2949–50. https://doi.org/10.1105/tpc.17.00929.

    CAS  Article  PubMed  Google Scholar 

  89. 89.

    Meyer KD, Patil DP, Zhou J, Zinoviev A, Skabkin MA, Elemento O, et al. 5' UTR m6A promotes cap-independent translation. Cell. 2015;163(4):999–1010. https://doi.org/10.1016/j.cell.2015.10.012.

  90. 90.

    Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, et al. N6-methyladenosine modulates messenger RNA translation efficiency. Cell. 2015;161(6):1388–99. https://doi.org/10.1016/j.cell.2015.05.014.

  91. 91.

    Coots RA, Liu XM, Mao Y, Dong L, Zhou J, Wan J, et al. m6A facilitates eIF4F-independent mRNA translation. Mol Cell. 2017;68:504–514.e507.

  92. 92.

    Li A, Chen YS, Ping XL, Yang X, Xiao W, Yang Y, et al. Cytoplasmic m6A reader YTHDF3 promotes mRNA translation. Cell Res. 2017;27(3):444–7. https://doi.org/10.1038/cr.2017.10.

  93. 93.

    Shi H, Wang X, Lu Z, Zhao BS, Ma H, Hsu PJ, et al. YTHDF3 facilitates translation and decay of N6-methyladenosine-modified RNA. Cell Res. 2017;27(3):315–28. https://doi.org/10.1038/cr.2017.15.

  94. 94.

    Choi J, Ieong KW, Demirci H, Chen J, Petrov A, Prabhakar A, et al. N6-methyladenosine in mRNA disrupts tRNA selection and translation-elongation dynamics. Nat Struct Mol Biol. 2016;23(2):110–5. https://doi.org/10.1038/nsmb.3148.

  95. 95.

    Qi ST, Ma JY, Wang ZB, Guo L, Hou Y, Sun QY. N6-methyladenosine sequencing highlights the involvement of mRNA methylation in oocyte meiotic maturation and embryo development by regulating translation in Xenopus laevis. J Biol Chem. 2016;291(44):23020–6. https://doi.org/10.1074/jbc.M116.748889.

    CAS  Article  PubMed  Google Scholar 

  96. 96.

    Slobodin B, Han R, Calderone V, Vrielink J, Loayza-Puch F, Elkon R, et al. Transcription impacts the efficiency of mRNA translation via co-transcriptional N6-adenosine methylation. Cell. 2017;169:326–337.e312.

  97. 97.

    Pontier D, Picart C, El Baidouri M, Roudier F, Xu T, Lahmy S, et al. The m6A pathway protects the transcriptome integrity by restricting RNA chimera formation in plants. Life Sci Alliance. 2019;2(3):e201900393. https://doi.org/10.26508/lsa.201900393.

  98. 98.

    Feng Z, Kang H, Li M, Zou L, Wang X, Zhao J, et al. Identification of new rice cultivars and resistance loci against rice black-streaked dwarf virus disease through genome-wide association study. Rice. 2019;12(1):49. https://doi.org/10.1186/s12284-019-0310-1.

  99. 99.

    Zhou T, Du L, Wang L, Wang Y, Gao C, Lan Y, et al. Genetic analysis and molecular mapping of QTLs for resistance to rice black-streaked dwarf disease in rice. Sci Rep. 2015;5(1):10509. https://doi.org/10.1038/srep10509.

  100. 100.

    Yan F, Zhang H, Adams MJ, Yang J, Peng J, Antoniw JF, et al. Characterization of siRNAs derived from rice stripe virus in infected rice plants by deep sequencing. Arch Virol. 2010;155(6):935–40. https://doi.org/10.1007/s00705-010-0670-8.

  101. 101.

    Zhang K, Zhang Y, Yang M, Liu S, Li Z, Wang X, et al. The barley stripe mosaic virus γb protein promotes chloroplast-targeted replication by enhancing unwinding of RNA duplexes. Plos Pathog. 2017;13(4):e1006319. https://doi.org/10.1371/journal.ppat.1006319.

  102. 102.

    Zhang K, Niu S, Di D, Shi L, Liu D, Cao X, et al. Selection of reference genes for gene expression studies in virus-infected monocots using quantitative real-time PCR. J Biotechnol. 2013;168(1):7–14. https://doi.org/10.1016/j.jbiotec.2013.08.008.

  103. 103.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  104. 104.

    Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology. 2013;14(4):R36. https://doi.org/10.1186/gb-2013-14-4-r36.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  105. 105.

    Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol. 2013;31(1):46–53. https://doi.org/10.1038/nbt.2450.

    CAS  Article  PubMed  Google Scholar 

  106. 106.

    Kersey PJ, Allen JE, Allot A, Barba M, Boddu S, Bolt BJ, et al. Ensembl Genomes 2018: an integrated omics infrastructure for non-vertebrate species. Nucleic Acids Res. 2017;46:D802–8.

  107. 107.

    Cui X, Zhang L, Meng J, Rao MK, Chen Y, Huang Y. MeTDiff: A novel differential RNA methylation analysis for MeRIP-seq Data. IEEE/ACM Transact Comput Biol Bioinformatics. 2018;15(2):526–34. https://doi.org/10.1109/TCBB.2015.2403355.

    CAS  Article  Google Scholar 

  108. 108.

    Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 2015;31(14):2382–3. https://doi.org/10.1093/bioinformatics/btv145.

    CAS  Article  PubMed  Google Scholar 

  109. 109.

    Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, et al. MEME Suite: tools for motif discovery and searching. Nucleic Acids Res. 2009;37(Web Server):W202–8. https://doi.org/10.1093/nar/gkp335.

  110. 110.

    Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38(4):576–89. https://doi.org/10.1016/j.molcel.2010.05.004.

  111. 111.

    Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2. https://doi.org/10.1093/bioinformatics/btq033.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  112. 112.

    Jain M. Genome-wide identification of novel internal control genes for normalization of gene expression during various stages of development in rice. Plant Sci. 2009;176(5):702–6. https://doi.org/10.1016/j.plantsci.2009.02.001.

    CAS  Article  Google Scholar 

  113. 113.

    Zhu Y, Zhou G, Yu X, Xu Q, Wang K, Xie D, et al. LC-MS-MS quantitative analysis reveals the association between FTO and DNA methylation. Plos One. 2017;12(4):e0175849. https://doi.org/10.1371/journal.pone.0175849.

  114. 114.

    Zhang K, Zhuang XJ, Dong, ZZ, Xu K, Chen XJ, Liu F, He Z: The dynamics of N6-methyladenine RNA modification in interactions between rice and plant viruses. Datasets. NCBI Bioproject. 2021. https://www.ncbi.nlm.nih.gov/bioproject/PRJNA732795. Accessed 31 May 2022.

Download references

Acknowledgements

We would very grateful to professor Chenggui Han (China Agricultural University, Beijing, China) for generously providing the RBSDV p10 antiserum in serological detection of RBSDV and the RBSDV-infected rice plants. Besides, we also very grateful to Song Shi and Ling Liu’s (OEbiotech, Shanghai, China) assistances in the sequencing data handling and arrangement.

Review history

The review history is available as Additional file 3.

Peer review information

Wenjing She was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Funding

This work was supported by the National Natural Science Foundation of China (31801699, 31601604, and 81903764), the National Science Foundation of Jiangsu Province (BK20180904), the Project of Science and Technology Development Plan for Traditional Chinese Medicine of Jiangsu Province (YB201993), the Qing Lan Project of Yangzhou University, the programme “Promoting Talen’s” launched by Jiangsu Association for Science and Technology.

Author information

Affiliations

Authors

Contributions

Zhen He and Kun Zhang conceived the project, designed and performed the experiments, analysed the data, and wrote the manuscript. Zhen He, Kun Zhang, and Fang Liu supervised the project. All other co-authors helped in performing the experiments and analysing the data. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Fang Liu or Zhen He.

Ethics declarations

Ethics approval and consent to participate

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: Fig. S1

Distribution pattern of sequenced m6A-IP-seq reads along the transcripts. Fig. S2 Intersection among m6A peaks identified in two biological replicates of three treatments. Fig. S3 Relative expression of jasmine acid (JA) biosynthesis and response genes in rice infected with viruses. Fig. S4 Relative expression of salicylic acid (SA) biosynthesis and response genes in rice infected with viruses. Fig. S5 Relative expression levels of abscisic acid (ABA) biosynthesis and deactivation genes in rice infected with viruses. Fig. S6 Relative expression levels of auxin biosynthesis, transportation, and signaling genes in rice infected with viruses. Fig. S7. Relative expression levels of cytokinin (CTK) biosynthesis, oxidation, and response genes in rice infected with viruses. Fig. S8 Relative expression levels of ethylene (ET) biosynthesis genes in rice infected with viruses. Fig. S9 Relative expression levels of brassinosteroids (BR) biosynthesis and signaling genes in rice infected with viruses. Fig. S10 Integrated analyses of the seven main phytohormone related genes with m6A methylation and relative expression levels in rice infected with viruses. Fig. S11 Loading control that corresponding to the dot-blot analyses. Fig. S12 The uncropped western blot membranes that showed in Fig. 1B and C.

Additional file 2: Table S1.

Primers used for RSV and RBSDV detections in this study. Table S2. Sequenced and rice genome mapped reads in m6A-IP-seq, input RNA-seq rice samples. Table S3. Nucleotide localization and enrichment of the top 10 m6A peaks identified in Mock-, RSV-, and RBSDV-infected rice transcripts by m6A-IP-Seq. Table S4. The category of the differential m6A peaks upon two viruses' infection in rice. Table S5. Nucleotide localization and enrichment of the m6A peaks identified in RSV and RBSDV genomics by m6A-IP-Seq. Table S6. Gene ID and their fpkm analyses. Table S7. Different peaks statistic. Table S8. The m6A peaks appeared in different treatments. Table S9. The most abundant consensus motif in Mock, RSV-, and RBSDV-infected rice plant using suite of DREME and MEME. Table S10. Analyses of the m6A peaks that containing most four common consensus appeared in other species. Table S11. Detail information of ListHits Gene and m6A methylated genes under rice viruses’ infection derived from Additional file 2: Table S16. Table S12. Integrated analyses of the m6A methylation related genes with m6A modifications and expression profiles. Table S13. Integrated analyses of the anti-viral RNA silencing pathway related genes with m6A modification and expression profile. Table S14. Integrated analyses of the plant hormone metabolic genes with m6A modifications and expression profiles. Table S15. Integrated analyses of the relationship betwwen relative expression and m6A positions. Table S16. Summary of the m6A RNA methylation level in enriched top 5 KEGG pathways related genes from RNA-seq under rice viruses’ infection. Table S17. Detail information of the common m6A methylated genes appeared in enriched top 5 KEGG pathways under rice viruses' infection. Table S18. Primers used for qRT-PCR validation of the methylation of OsAGO18 and OsSLRL1 genes. Table S19. Primers used in qRT-PCR qualification of the candidate genes.

Additional file 3:.

Review history.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Zhang, K., Zhuang, X., Dong, Z. et al. The dynamics of N6-methyladenine RNA modification in interactions between rice and plant viruses. Genome Biol 22, 189 (2021). https://doi.org/10.1186/s13059-021-02410-2

Download citation

Keywords

  • N 6-methyladenosine
  • Rice
  • Plant viruses
  • Interactions