Transcriptome-wide N6-methyladenosine (m6A) methylation in watermelon under CGMMV infection
BMC Plant Biology volume 21, Article number: 516 (2021)
Cucumber green mottle mosaic virus (CGMMV) causes substantial global losses in cucurbit crops, especially watermelon. N6-methyladenosine (m6A) methylation in RNA is one of the most important post-transcriptional modification mechanisms in eukaryotes. It has been shown to have important regulatory functions in some model plants, but there has been no research regarding m6A modifications in watermelon.
We measured the global m6A level in resistant watermelon after CGMMV infection using a colorimetric method. And the results found that the global m6A level significantly decreased in resistant watermelon after CGMMV infection. Specifically, m6A libraries were constructed for the resistant watermelon leaves collected 48 h after CGMMV infection and the whole-genome m6A-seq were carried out. Numerous m6A modified peaks were identified from CGMMV-infected and control (uninfected) samples. The modification distributions and motifs of these m6A peaks were highly conserved in watermelon transcripts but the modification was more abundant than in other reported crop plants. In early response to CGMMV infection, 422 differentially methylated genes (DMGs) were identified, most of which were hypomethylated, and probably associated with the increased expression of watermelon m6A demethylase gene ClALKBH4B. Gene Ontology (GO) analysis indicated quite a few DMGs were involved in RNA biology and stress responsive pathways. Combined with RNA-seq analysis, there was generally a negative correlation between m6A RNA methylation and transcript level in the watermelon transcriptome. Both the m6A methylation and transcript levels of 59 modified genes significantly changed in response to CGMMV infection and some were involved in plant immunity.
Our study represents the first comprehensive characterization of m6A patterns in the watermelon transcriptome and helps to clarify the roles and regulatory mechanisms of m6A modification in watermelon in early responses to CGMMV.
N6-methyladenosine (m6A) methylation is the most prevalent internal modifier of eukaryote RNA and is widely present. It plays multiple roles in RNA metabolism, including stability and degradation , alternative splicing , translation efficiency , and nuclear export . In plants, the dynamic patterning of m6A modifications in mRNA is mediated by three clades of elements, including m6A writer, eraser, and reader proteins . The m6A methylation is formed by the m6A methyltransferase writer, which is complex composed of N6-adenosine-methyltransferase MT-A70-like (MT) proteins MTA, MTB, and other factors. In Arabidopsis, removal of m6A is catalyzed by m6A demethylases of alkylated DNA repair proteins AlkB homologs (ALKBHs) such as ALKBH9B and ALKBH10B. The discovery of methyltransferases and demethylases implies the m6A modification is a reversible and dynamic RNA methylation. A clade of proteins containing YTH-domain serves as m6A readers such as Arabidopsis EVOLUTIONARILY CONSERVED C-TERMINAL REGION2/3/4 (ECT2/3/4) and Cleavage and Polyadenylation Specificity Factor 30 (CPSF30) proteins, which could bind specifically to m6A sites (such as RRACH) and mediate specific functions.
Evidence indicates that m6A RNA methylation is highly conserved in plants and involved in multiple developmental and biological processes, such as embryonic development [6, 7], trichome branch development [8,9,10], leaf initiation , flowering and root development . Moreover, a few reports suggest that m6A methylation plays important roles in plant antiviral immunity. For example, Arabidopsis m6A demethylase ALKBH9B can remove m6A modifications from viral RNA and repress Alfalfa mosaic virus (AMV) infections . In tobacco, the global m6A level was reduced under Tobacco mosaic virus (TMV) infection, probably associated with decreased m6A methyltransferase and increased demethylase expression . Furthermore, m6A levels in susceptible rice increased in response to Rice black-streaked dwarf virus (RBSDV) and Rice stripe virus (RSV) infection, as m6A methylation mainly depends on genes which are not actively expressed in virus-infected rice plants via m6A-seq analyses . In wheat, the transcriptome-wide m6A profile of two varieties with different resistance to Wheat yellow mosaic virus (WYMV) revealed that many genes related to plant defense responses and plant-pathogen interaction significantly changed in both m6A and mRNA levels . Overall, these studies suggested that m6A modification plays important roles in plant response to viral infection, but the relevant regulatory mechanism is complex and still needs to be elucidated.
Cucumber green mottle mosaic virus (CGMMV), a member of the genus Tobamovirus, has spread rapidly and caused huge economic losses in cucurbits especially watermelon. Consequently, it is important to determine host response mechanisms of watermelon to CGMMV. Previously, the transcriptome, and proteomic, microRNA (miRNA), and virus derived small interfering RNAs (vsiRNAs) profiles via high-throughput sequencing have been performed, and some CGMMV-responsive genes, miRNAs, and vsiRNAs were identified in watermelon [17,18,19]. However, as the most important post-transcriptional mechanism, m6A modification profiles in watermelon following CGMMV infection have not been explored.
In this study, we found that the global m6A level in resistant watermelon clearly decreased after CGMMV infection by colorimetry. Therefore, we analyzed watermelon m6A methylation profiles in response to CGMMV infection via whole-genome m6A-seq and bioinformatics approaches. Meanwhile, combined with RNA-sequencing, the response patterns and putative functions of differentially expressed m6A-modified genes (DMGs) derived from the transcriptomes of CGMMV-infected watermelon leaves were analyzed. Furthermore, we verified some target genes were verified by m6A-immunoprecipitation (IP)-qPCR (m6A-IP-qPCR). Our study represents the first comprehensive characterization of m6A modification in the watermelon transcriptome. Furthermore, we attempted to elucidate the potential m6A-mediated roles and regulatory mechanisms in watermelon responses to CGMMV.
CGMMV infection decreased the m6A global level in resistant watermelon plants
To explore whether m6A modification responded to CGMMV infection, we assessed the mRNA m6A methylation levels in the control and CGMMV-infected leaves from resistant and susceptible watermelon (R1288 and S1252) by colorimetry (Fig. 1). The m6A level weakly decreased to 0.923-fold in resistant watermelon plants but significantly increased to 1.397-fold in susceptible watermelon plants 24 h after CGMMV infection. However, at 48 h, the m6A level in resistant watermelons significantly decreased to 0.757-fold, and the m6A level in susceptible watermelons weakly decreased compared to control samples. Therefore, CGMMV infection reduced leaf m6A levels of resistant watermelon.
Transcriptome-wide mapping of m6A in watermelon
To investigate the transcriptome-wide map of m6A methylation in watermelon, a series of m6A- immunoprecipitation (IP) and matched input (non-IP control) libraries were constructed and sequenced (Supplementary file 5). This series included control and CGMMV-infected watermelon leaves, with two biological replicates each. High Pearson correlation coefficient (PCC) were determined for the abundance of confident peaks between biological replicates, representing highly reproducible results (Supplementary file Figure. S1). A total of 65.59-86.12 million reads were generated for each library, and after filtering out low-quality data, 65.11-85.41 million high-quality reads were mapped to watermelon reference genome (Charleston Gray). Furthermore, 63.03-82.59 million distinct reads were uniquely aligned to the genome and 22.85-30.35 million were mapped to splice reads (Table 1). Only m6A peaks consistently detected in both sample biological replicates were used for subsequent analyses.
m6A methylation is a common feature of watermelon mRNAs
At the genome level, 13,004 and 13,123 m 6A peaks with high confidence, distributed in 8,899 and 9,200 genes, were identified from the control and CGMMV-infected leaves, respectively, with an average of 1.45 m6A peaks within transcription units for each modified gene (Fig. 2A, Supplementary file 6). The enrichment multiple distribution analyses showed that more reads were enriched in the peaks of the CGMMV-infected watermelon than those of the control samples (Fig. 2B). The average length of the peaks was 671.38 and 689.34 bp in the control and CGMMV-infected samples, respectively (Fig. 2C).
The frequency of m6A modifications was unevenly distributed within RNA, being particularly highly enriched in mature mRNAs [20, 21]. The watermelon m6A peaks were highly abundant around the coding sequences (CDS) (58.89%-64.25%) and 3’-untranslated regions (3’-UTRs) (29.98%-36.22%) in all experimental samples (Fig. 3A and B). To identify enriched sequence motifs within watermelon m6A peaks, hypergeometric optimization of motifs enrichment were applied using the HOMER suite. The results indicated the most significant motif identified in watermelon m6A peaks was an RRACH motif (where R represents A/G, A is m6A, H represents A/C/U, and Y represents A, G, U, or C) with an E-value of 1e-44 and 1e-93 in control and CGMMV infected samples, respectively. The motif was found in 75.26% of all watermelon m6A peaks. Meanwhile, a URUAY motif was significantly enriched in 72.38% of watermelon m6A peaks (where Y represents C/U), with an E-value of 1e-59 and 1e-17 in control and CGMMV infected samples, respectively (Fig. 3C). These results suggest that m6A distribution and motifs in watermelon leaves are conserved, as reported in other species [11, 22, 23].
m6A-modified genes in watermelon are involved in multiple signaling pathways and cellular processes
Gene ontology (GO) analysis indicated that watermelon m6A-modified transcripts were grouped into three categories: biological process, cellular component, and molecular_function (Supplementary file 8; Figure S2). Furthermore, these m6A-modified transcripts participated in multiple signaling pathways and cellular processes. The significantly enriched GO terms in biological process were involved in mRNA metabolism including mRNA processing (GO:0006397), rRNA processing (GO:0006364), mRNA splicing (GO:0000398, GO:0008380, and GO:0045292). For the cellular component, most of these modified genes were significantly enriched in the cytosol (GO:0005829), nucleolus (GO:0005730), and nucleus (GO:0005634) (Supplementary file 7; Figure S2). A The most enriched molecular functions were ATP binding (GO:0005524), RNA binding (GO:0003723), mRNA binding (GO:0003729), protein serine/threonine kinase activity (GO:0004674) and so on.
Most watermelon differentially methylated peaks (DMPs) were hypo-methylated in response to CGMMV infection
MeTDiff was used to compare the abundance of m6A peaks between control and CGMMV-infected samples, and DMPs were defined based on a cutoff criterion of FPKM fold change ≥ 1.5 and P value < 0.05. Accordingly, 430 DMPs were identified, they were distributed among 422 differentially methylated genes (DMGs) (Supplementary file 8), and of these DMPs, 255 were hypo-methylated m6A peaks, whereas 175 were hyper-methylated m6A peaks. Moreover, 295, 115, and 21 DMPs were located in CDS, 3’-UTRs, and 5’-UTRs, respectively (Supplementary file 8; Figure S3).
We performed GO analysis of the DMGs to explore their potential functional roles in response to CGMMV infection, and found that DMGs were involved in some important biological pathways such as cellular process, metabolic process, response to stimulus, cellular component organization or biogenesis, biological regulation, and so on (Fig. 4). Notably, 11 DMGs were involved in immune responses (Supplementary file 9-1), and were annotated as follows: zinc finger protein (ClCG01G014940), eukaryotic aspartyl protease family protein (ClCG05G026860), ankyrin repeat-containing protein (ClCG02G016860), 30 kDa cleavage and polyadenylation specificity factor subunit 4-like protein (ClCG03G006290), RNA-binding protein (ClCG09G021330), 26S protease regulatory subunit (ClCG00G004170), MLO-like protein (ClCG00G002750), CPR5 protein (ClCG01G019270), tetratricopeptide repeat protein (ClCG01G017250), ABC transporter family protein (ClCG07G013290), and late embryogenesis abundant (LEA) protein (ClCG10G018420). All of their m6A methylation levels were significantly downregulated from 0.042- to 0.666-fold except the genes encoding MLO-like, CPR5, and LEA genes (Supplementary file 9-1). Moreover, three DMGs (ClCG03G016200, ClCG05G016300, and ClCG10G003630), were annotated as DNA-directed RNA polymerase, RNA-dependent RNA polymerase, and DNA-directed RNA polymerase IV subunit 1 protein, respectively (Supplementary file 9-2). They all participated in the production of siRNA involved in RNA interference (GO:0030422) and m6A methylation levels increased from 1.528- to 2.042-fold after CGMMV infection (Supplementary file 9-2).
Most watermelon differentially expressed genes (DEGs) were induced by CGMMV infection
We performed RNA-seq analysis to quantify the expression levels of m6A elements in resistant watermelon leaves in response to CGMMV infection. Genes with fold change ≥ 1.5 and P value < 0.05 between libraries were considered to be significantly differentially expressed genes (DEGs). A total of 1,855 DEGs were identified, of which 1,230 were upregulated and 625 were downregulated 48 h after CGMMV infection (Supplementary file 10). Among these DEGs, 45.88% contained m6A peaks, which was significantly higher than the rate for all of watermelon transcripts (39.43% - 40.76%). GO enrichment analysis (top 10 enrichment scores) indicated that these DEGs were significantly enriched in various biological process pathways including defense response (GO:0006952), response to oxidative stress (GO:0006979), defense response to bacterium (GO:1900426 and GO:0042742), response to chitin (GO:0010200), DNA replication initiation (GO:0006270), and so on. In molecular_function, the DEGs were significantly enriched in pathways related to heme binding (GO:0020037), iron ion binding (GO:0005506), calmodulin binding (GO:0005516), peroxidase activity (GO:0004601), and chitinase activity (GO:0004568). Besides, they highly located in plasma membrane (GO:0005886) (Fig. 5).
m6A machinery components involved in regulation of m6A modification in virus-infected watermelon plants
All watermelon m6A machinery components including three m6A methylases (ClMTA, ClMTB, and ClMTC), twelve demethylases (ClALKBH1-ClALKBH12), and five m6A readers (ClECT1- ClECT5), were identified from watermelon genome (Supplementary file 11). None of these m6A gene expressions showed statistically significant differences in response to CGMMV infection via RNA-seq analysis, except two m6A demethylases, ClALKBH4B (ClCG04G006530) and ClALKBH11B (ClCG11G013640). After infection, the expression level of ClALKBH11B dropped to 0.493-fold, but that of ClALKBH4B expression increased to 1.519-fold. Meanwhile, the expression patterns of these two genes were also verified by qRT-PCR in resistant and susceptible watermelon varieties under CGMMV infection. Similarly, ClALKBH4B expression increased 7.086-fold in resistant varieties, but did not change (0.931-fold) in susceptible varieties 48 h after CGMMV infection. ClALKBH11B expression decreased to 0.142-fold in resistant watermelon, but increased to 1.828-fold in susceptible watermelon (Fig. 6). Moreover, a m6A reader, ClECT1 (ClCG03G006290), exhibited significantly decreased m6A (0.666-fold) under CGMMV infection in m6A-seq (Supplementary file 9).
Conjoint analysis of DMGs and DEGs responsive to CGMMV infection
To explore a potential correlation between transcript levels and m6A mRNA methylation in response to CGMMV infection, a conjoint analysis of DMGs and DEGs was conducted and 59 differentially methylated and expressed genes (DMEGs) were identified. The DMEGs were mainly divided into four groups, including 22 hypomethylated but upregulated genes (hypo-up) (37.29%), 16 hypermethylated as well as upregulated genes (hyper-up) (27.12%), 17 hypomethylated as well as downregulated genes (hypo-down) (28.81%), and 6 hypermethylated but downregulated genes (hyper-down) (10.17%) (Supplementary file 12). Functional annotation analyses indicated that these genes were involved in multiple roles and signaling pathways such as resistance response, secondary biosynthesis and metabolism, and RNA processes. Among them, six transcript factors (TFs) including WRKY DNA-binding protein 29 (WRKY29) (ClCG07G007870), ethylene-responsive transcription factor (ERF) (ClCG10G017140), MYB domain protein 68 (MYB68) (ClCG08G016830), heat stress transcription factor (HSF) (ClCG01G016810), GATA transcription factor (ClCG07G008920), and zinc finger (bZIP) (ClCG08G016400) were identified. Except for WRKY29, all of their m6A methylation levels decreased from 0.55- to 0.65-fold, but transcript patterns displayed diverse responses. The transcript levels of WRKY29, HSF, and bZIP were upregulated 1.85- to 3.39-fold, but MYB68, ERF, and GATA were downregulated from 0.14- to 0.57-fold. Seven genes involved in sugar metabolism and signaling encoding were identified , including alpha-l-fucosidase (ClCG01G001290), glycosyltransferase family protein (ClCG01G025940), C2 Calcium/lipid-binding plant phosphoribosyl transferase family protein (ClCG04G005410), O-fucosyltransferase family protein (ClCG10G021720), glutamine-fructose-6-phosphate amino transferase (ClCG11G009710), UDP-D-GLUCURONATE 4-EPIMERASE 3 family protein (ClCG10G003090), and endo-1,4-beta-glucanase (ClCG11G011070). Among these seven genes, the m6A modification levels of four genes including ClCG01G001290, ClCG04G005410, ClCG01G025940, and ClCG10G003090 were significantly downregulated from 0.40- to 0.60-fold, while the other three genes were upregulated from 1.55- to 2.32-fold. Five plant resistance genes were identified, including encoded receptor lectin kinase-like protein (ClCG10G022400), late embryogenesis abundant (LEA) hydroxyproline-rich glycoprotein protein (ClCG10G018420), ankyrin repeat-containing proteins (ClCG02G016860 and ClCG01G009830), and ABC transporter family pleiotropic drug resistance protein (ClCG07G013290). Among them, the m6A methylation levels of ClCG10G018420 and ClCG02G003390 significantly increased but the other three genes were decreased. Meanwhile, except for ClCG01G009830, all the four gene expressions were induced. Five genes participated in genetic information processing, among them, ClCG09G005490 and ClCG08G009480 were involved in transcription, with increased m6A modification and decreased transcript levels, And ClCG01G025960 involved in replication and repair, showed increased m6A modification and transcript levels (Supplementary file 12).
Validation of m6A methylation and transcript levels of DEGs
We confirmed the m6A methylation and transcript levels of WRKY29, MYB68, and bZIP in watermelon by m6A-IP-qPCR and qRT-PCR. As expected, methylation levels of MYB68 and bZIP genes were downregulated 0.19- and 0.25-fold, respectively, whereas the methylation level of WRKY29 was upregulated 2.22-fold. Meanwhile, the expression levels of WRKY29 and bZIP were upregulated 5.62- and 3.12-fold (Fig. 7), but MYB68 expression was repressed (0.31-fold). These results confirm that our m6A-seq and RNA-seq data were accurate and robust.
Plants have developed complex molecular mechanisms to resist viruses. Some RNA epigenetic modifications involved in the post-transcriptional regulation of virus-induced genes play significant roles in plant antiviral processes. As one of the most important RNA methylations in eukaryotes, m6A methylation is important for distinct RNA processing steps  and was recently identified regulating virus infection in plants [13, 14]. However, presently there is no transcriptome-wide study on the m6A modification in early response to virus infection in plants. CGMMV, a species of the Tobamovirus genus, has caused considerable yield losses in commercial watermelon production worldwide. To determine the response patterns of RNA m6A modification to CGMMV infection in watermelon, libraries were constructed from watermelon leaves collected 48 h after CGMMV to identify the responsive genes at the m6A methylation level using high-throughput sequencing technology (m6A-seq).
The m6A modification features in watermelon
Our study is the first to generate and globally characterize the transcriptome-wide RNA m6A modification profiles in watermelon. For the watermelon genome, out of all transcripts, 8,899 to 9,200 genes (39.43-40.76%) were identified as m6A modified genes (Supplementary file 6). The proportion of modified genes in watermelon was higher than that in Arabidopsis (19.58-22.74%) and tomato (26.48-27.71%) [22, 23]. Additionally, over 30% of m6A-modified watermelon genes had more than one peak, which was obviously higher than the rate in tomato (8.27%) and maize (5.88%). To date, two conserved motifs, RRACH and URUAY, have been identified in plants. Specifically, URUAY is a plant-specific consensus motif enriched within m6A peaks generated for Arabidopsis, tomato, rice, and maize [11, 22, 23, 25]. In our study, m6A peaks in watermelon also contained the motifs of RRACH and URUAY, with frequency of 75.26 and 72.38%, respectively. Further studies indicate that the enrichment E-values of the RRACH motif were clearly lower than those of the URUAY motif in watermelon and Arabidopsis  (Fig. 3C), but the opposite result was observed in maize. The different m6A site biases between monocot (maize) and dicot (Arabidopsis and watermelon) plants is probably related to their unique differences in structure and function of m6A methylation. Almost all watermelon m6A peaks (94.22-95.11%) were located in 3’-UTRs and CDS (Fig. 3A and B), which is similar to the results in Arabidopsis, tomato, rice, and maize [11, 22, 23]. GO analysis indicated that m6A modified-genes in watermelon were significantly enriched in multiple signaling pathways and cellular processes, implying their multiple roles (Fig. 4, Supplementary file 7). Collectively, m6A is a common feature in a substantial fraction of watermelon mRNA. m6A modification distributions and characteristics were highly conserved in watermelon, but more abundant than in other reported plants.
Potential roles of m6A modification in early response to CGMMV in watermelon
Studies have found that the m6A methylation level of plant hosts is reprogrammed in response to viral infection [13, 14]. We found the global m6A level was regulated by CGMMV infection (Fig. 1), it was more significantly downregulated in resistant, than susceptible, watermelons. This implies that a global m6A-hypomethylation probably activates watermelon defense responses at the early stage of CGMMV infection. Furthermore, 422 DMGs responsive to CGMMV infection were screened out by m6A-seq, of which 60.42% were the hypomethylated after the CGMMV infection. These results suggest that m6A modifications might negatively affect watermelon responses to CGMMV even further. Meanwhile, 1,855 DEGs responsive to CGMMV infection were identified in watermelon via RNA-seq, of which 66.31% were induced. This is consistent with previous findings that m6A methylation is generally negatively correlated with transcript abundance .
The m6A methylation is directly controlled by m6A methylases, demethylases, and binding (reader) proteins. In watermelon, we found a m6A reader, homologous with Arabidopsis ECT1/CPSF30, was hypomethylated by m6A modification after CGMMV infection (Supplementary file 9). In Arabidopsis, ECT1/CPSF30 played important roles in immune responses by regulating the splicing of the 3’ end of mRNA . In tobacco, the m6A global level was decreased because the expression of some ALKBHs were increased after TMV infection . In watermelon, the transcript levels of two m6A demethylases, ClALKBH4B and ClALKBH11B, were significantly changed in RNA-seq and qRT-PCR data under CGMMV infection (Fig. 6, Supplementary file 11). ClALKBH4B expression was significantly induced by CGMMV infection in resistant watermelon plants, but seemed to be unaffected in susceptible plants. Therefore, we suspected that the decreased expression of ClALKBH4B resulted in the hypomethylation of most DMGs in m6A-seq.
Virus-induced gene silencing, as one type of posttranscriptional gene silencing (PTGS), is a common antiviral strategy used by plants. Virus infection is initiated through the transcription catalyzed by an RNA-dependent RNA polymerase to produce double-stranded RNA (dsRNA) [27, 28]. The dsRNA is cleaved to generate small interfering RNA (siRNA), which is recruited by host RNA-induced silencing complexes to repress viral gene transcription and translation. Arabidopsis AtALKB9B was found to colocalize with siRNA-body and is implicated in PTGS through the synthesis of dsRNAs to generate siRNAs (Martínez-Pérez et al. 2017). We found three DMGs including ClCG03G016200, ClCG05G016300, and ClCG10G003630, participated in the production of siRNA and all of them were m6A-hypermethylated after CGMMV infection (Supplementary file 9). This suggests that m6A methylation in watermelon was probably involved in virus-induced gene silencing to resist CGMMV infection.
In terms of plant-virus interactions, plants have evolved immune systems that recognize the presence of viruses and initiate the expression of effective defense response factors. Recently, two studies have reported on transcriptome and proteomic profiles of CGMMV-infected and mock-inoculated watermelon leaves or fruit [19, 29]. We compared them with our RNA-seq data and found that the number of DEGs identified in previous RNA-seq data (1641 DEGs) was similar to our study (1855 DEGs). Most of the DEGs or differentially accumulated proteins (DAPs) in the three studies were induced by CGMMV infection and functional analyses indicated they enriched in the pathways involved in stress response, photosynthesis, secondary metabolism, and so on [19, 29]. These results further validated the common pathways in watermelon involved in CGMMV response. To further explore the relationship between the m6A methylation extent and the transcript level in response to CGMMV infection, we identified 59 genes as DMEGs and they were important candidate CGMMV-responsive genes at transcriptional and post-transcriptional (m6A modification) levels. Most of their homologous genes have been found to participated in plant innate immunity. For example, six TFs were identified and all of their m6A methylation levels decreased except WRKY29, and WRKY29, HSF, and bZIP transcript levels increased in response to CGMMV infection (Supplementary file 12). These TFs have been proven to participate in plant-pathogen interaction pathways. WRKY TFs were involved in the defense signaling pathway associated with Tomato yellow leaf curly virus infection . Arabidopsis HSF was involved in the regulation of Pdf1.2 expression and pathogen resistance . And bZIP TFs were important for plant innate immunity as they regulated the expression of genes associated with PAMP- and, effector-triggered immunity, and hormone signaling networks [32, 33].
Plant antiviral mechanisms result in the regulation of plant carbohydrate allocation and signaling. Soluble sugars and starch accumulate in the leaves to initiate viral replication, resulting in decreased photosynthesis and increased respiration [34, 35]. We identified seven DMEGs involved in sugar metabolism and signaling (Supplementary file 12). Most of them were m6A-hypomethylated but the transcript levels were upregulated after CGMMV infection. Most of their homologous genes were found to have participated in plant innate immunity. Notably, we compared these DMEG with the previous RNA-seq results and found two DMEGs (ClCG11G011070 and ClCG10G021720), both involved in suger metabolism and signalling, were also significantly induced by CGMMV infection in the previous RNA-seq data . These m6A-regulated genes probably affect soluble sugar and starch accumulation to participate in CGMMV response.
It's worth mentioning that, only two genes (ClCG02G022850 and ClCG09G016700), both encoded ABC transporters, out of all DMEGs contained two DMPs, and their transcript levels were induced under CGMMV infection, while the two m6A peaks in ClCG02G022850 were significantly hypermethylated, but only one peak in ClCG09G016700 were hypermethylated in response to CGMMV infection (Supplementary file 12). ABC transporter proteins have been found to play important roles in resisting cell wall penetration and haustorium formation, and improving the resistance to some fungal and oomycete pathogens . These studies suggested that the m6A modified ABC transporter in watermelon probably play important roles in CGMMV responses via a complex regulatory mechanism.
m6A modification is essential for plant growth and development and stress responses, including plant antiviral immunity. In this study, we performed m6A-seq to analyze transcriptome-wide m6A methylation profile of watermelon leaves under CGMMV infection. And the results indicated the distributions and motifs of m6A methylation were highly conserved while the modification rates was obviously higher in watermelon than that in other reported plants. In further, we proposed a preliminary hypothesis to explain the mechanism of watermelon resistance to viral infection via the m6A modification; m6A demethylase gene ClALKBH4B could be significantly induced in early response to CGMMV in resistant watermelon, and the decreased expression of ClALKBH4B resulted in the methylation of numerous target modified-genes were downregulated. The m6A methylation of the transcripts was generally negatively correlated with transcriptional level as reported previously , so some downstream defense response factors involved in virus-induced gene silencing, TFs, plant carbohydrate allocation and signaling, and so on, were induced at transcript level via RNA-seq analysis and a series of plant immune responses were activated at the early stage of CGMMV infection. This study firstly provide insights into the m6A modification features in watermelon leaf and reveal the roles and regulatory mechanisms of m6A modification under CGMMV infection in watermelon.
Plant materials and CGMMV inoculation
We used two watermelon germplasms, namely, R1288 (a highly CGMMV-resistant watermelon germplasm, R) and S1252 (a highly CGMMV-susceptible watermelon germplasm, S), that were identified through resistance screening (Supplementary file Fig. S4). Watermelon seedlings were grown under standard culture conditions (25 °C and a 16-h light/8-h dark cycle) in a greenhouse. For CGMMV inoculation, seedlings at the two true-leaf stage of the two watermelon germplasms were agroinfiltrated with Agrobacterium carrying CGMMV infectious clone, as described previously . Meanwhile, the seedlings of the two watermelon germplasms were inoculated with sodium phosphate buffer (pH 7.2) without virus and were used as control (Mock). At 24 h post-inoculation, at least three leaves each from un-infected control and infected plants were collected and mixed to a biological replicate and three biological replicates were obtained for watermelon germplasm R1288 and S1252, respectively. Similarly, samples were collected at 48 h post-inoculation. Samples collected from the un-infected and infected plants of R1288 and S1252 at 24 and 48 h post-inoculation, respectively, were used for global m6A methylation quantification and quantitative RT-PCR (qRT-PCR) analyses. Samples collected from from the un-infected and infected R1288 plants at 48 h post-inoculation were used for high-throughput m6A-seq and RNA-seq.
Global mRNA m6A methylation assay
Global m6A methylation assay was performed as reported . Firstly, total RNA was extracted from harvested leaves by using RNAiso Plus (TaKaRa). And PolyA+mRNA selection was isolated using a Dyna beads mRNA Purification Kit (Life Technologies). Then we measured the global m6A levels in transcripts with an EpiQuik m6A RNA Methylation Quantification Kit (Colorimetric) (Epigentek) according to the manufacturer’s protocol.
Identification of watermelon m6A elements
The hidden Markov model (HMM) profile for the MT70 superfamily (PF05063), ALKBH superfamily (PF13532, clavaminate synthase-like domain), and YTH domain-containing protein (PF04146) sequences were downloaded from the PFAM website (http://pfam.xfam.org/), and the HMMER tool was applied to identify orthologous proteins in the watermelon reference genome (Charleston Gray) (http://cucurbitgenomics.org/). These m6A genes were named according to their counterparts in Arabidopsis.
High-throughput m6A-seq and RNA-seq
Resistant watermelon leaves (R1288) were collected at 48 h post-CGMMV inoculation and mock and used for the m6A-seq, which was carried out as reported . Specifically, total RNA was extracted and fragmented into approximately 200 nucleotide-long fragments. These fragments were incubated and purified as previous reported . m6A-seq was performed on the Illumina HiSeq platform at Shanghai OE Biotech Co., Ltd.
Sequencing data analysis
The raw reads were trimmed to remove adapters and low-quality bases with the Trimmomatic 0.36 , and their quality was assessed with the FastQC program (https://www.bioinformatics. babraham.ac.uk/projects/fastqc). The reads were aligned to the watermelon reference genome (Charleston Gray) with HISAT2 and annotated. All the uniquely mapped reads with a MAPQ ≥ 13 were retained for follow-up analyses .
Chipseeker software  was used for identifying m6A peaks in each anti-m6A IP sample, with the corresponding input sample serving as a control. A stringent cutoff threshold for the false discovery rate (FDR) < 0.05 was used to obtain high-confidence peaks. Only peaks consistently called in both two independent biological samples were considered as high-confidence peaks and used for subsequent analyses, and these peaks were annotated with PeakAnnotator (version 2.0) .
Differentially methylated peaks between samples were determined with MeTDiff software  using the following criteria: FoldChange ≥ 1.5 and P value < 0.05. The m6A-enriched motifs were identified with HOMER (version 4.7; http://homer.ucsd.edu/homer/) . All peaks mapped to mRNAs were used as target sequences, and exon sequences except for the peak-containing sequences were used as the background sequences. The motif length was restricted to six nucleotides and the E-value threshold was 1E-5. Visualization analysis of m6A peaks was conducted using Integrated Genome Browser (IGB, version 9.0.2) . Gene Ontology (GO) analysis of m6A-modified genes was performed with the Gene Ontology Consortium (http://www.geneontology.org/). GO terms with a Bonferroni corrected P value < 0.05 were considered to be significantly enriched for individual genes.
RNA-seq reads were aligned to the watermelon (Charleston Gray) reference genome, with default settings except for the maximum intron length, which was set to 10 kilobase (kb). Unique mapped reads were used as the input for Cufflinks v2.2.1 , which was used to normalize and estimate gene expression levels in terms of fragments per kilobase of transcript per million mapped reads (FPKM = counts of mapped fragments×109/ [transcript length×total count of the mapped fragments]). Differentially expressed genes (DEGs) were analyzed using the Cuffdiff program of Cufflinks . The watermelon (Charleston Gray) reference genome sequences and annotations were downloaded from the CuGenDB database (http://cucurbitgenomics.org/).
Validation of transcriptome data by m6A-IP-qPCR and qRT-PCR
To validate m6A methylation and expression levels revealed by the transcriptome data, watermelon leaves were collected at 48-h post-CGMMV inoculation and mock as described above and used for RNA extraction. The m6A-IP-qPCR was performed as previously described, with some modifications [23, 44]. The purified mRNA was fragmented into approximately 200 nucleotide sequences, and then stopped, and the mRNA fragments were purified by ethanol precipitation and resuspended. The fragmented mRNA was then resuspended and 100 μL of fragmented mRNAs used as the input sample (input mRNA). Another 100 μL of fragmented mRNAs were incubated with an anti-m6A polyclonal antibody. After washing twice with high-salt buffers and twice with IP buffers, the bound mRNAs were eluted from the beads, and immune-precipitated mRNA fragments were resuspended in 5 μL DEPC-treated water. Then, the immunoprecipitated mRNA and pre-immunoprecipitated mRNA (input mRNA) were reverse transcribed with random hexamers using M-MLV reverse transcriptase (Takara) to generate the cDNA for the IP and pre-IP samples. The cDNA was used as the template for the m6A-IP-qPCR and qRT-PCR. The enrichment of m6A methylation in specific gene regions was determined using the cycle threshold (CT) 2(-ΔCT) method . The β-actin gene was used as an internal control . The primers used for m6A-IP-qPCR and qRT-PCR are listed in Supplementary file 13.
Availability of data and materials
The raw sequencing data and processed peaks data for m6A-seq have been uploaded to the NCBI database under the SRA accession number: PRJNA597500 (https://www.ncbi.nlm.nih.gov/ bioproject/PRJNA597500). All the public databases used in this study are open. All the raw data generated in this study are included in the article and the supplemental files.
Cucumber green mottle mosaic virus
Wheat yellow mosaic virus
Rice black-streaked dwarf virus
Rice stripe virus
Alfalfa mosaic virus
Tobacco mosaic virus
Differentially methylated genes
Alkylated DNA repair proteins AlkB homologs
EVOLUTIONARILY CONSERVED C-TERMINAL REGION
Cleavage and Polyadenylation Specificity Factor 30
Alfalfa mosaic virus
Hidden Markov model
False discovery rate
Differentially expressed genes
Pearson correlation coefficient
- UTRs: DMEGs:
Untranslated regions: differentially methylated and expressed genes
Late embryogenesis abundant
Posttranscriptional gene silencing
Small interfering RNA
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:117.
Zhao X, Yang Y, Sun BF, Shi Y, Yang X, Xiao W, et al. FTO-dependent demethylation of N6-methyladenosine regulates mRNA splicing and is required for adipogenesis. Cell Res. 2014;24:1403–19.
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:315–28.
Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA modifications in gene expression regulation. Cell. 2017;169:1187–200.
Yue H, Nie X, Yan Z, Weining S. N6-methyladenosine regulatory machinery in plants: composition, function and evolution. Plant Biotechnol J. 2019;17:1194–208.
Bodi Z, Zhong S, Mehra S, Song J, Li H, Graham N, 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.
Růžička 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:157–72.
Lockhart J. A tale of three studies: uncovering the crucial roles of m6A readers. Plant Cell. 2018;30:947.
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:986–1005.
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:968–85.
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 Biol. 2014;11:1180–8.
Chen M, Urs MJ, Sanchez-Gonzalez I, Olayioye MA, Herde M, Witte CP. m6A RNA degradation products are catabolized by an evolutionarily conserved N6-methyl-AMP deaminase in plant and mammalian cells. Plant Cell. 2018;30:1511–22.
Martínez-Pérez M, Apariciob F, López-Gresa MP, Bellés JM, Sánchez-Navarro JA, Pallás V. 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–60.
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:87.
Zhang K, Zhuang X, Dong Z, Xu K, Chen X, Liu F, et al. The dynamics of N6-methyladenine RNA modification in interactions between rice and plant viruses. Genome Biol. 2021;22:189.
Zhang TY, Wang ZQ, Hu HC, Chen ZQ, Liu P, Gao SQ, et al. Transcriptome-WideN6-Methyladenosine (m6A) profiling of susceptible and resistant wheat varieties reveals the involvement of variety-specific m6A modification involved in virus-host interaction pathways. Front. Microbiol. 2021;12:656302.
Sun Y, Niu X, Di C, Fan M. High-throughput sequencing reveals vsiRNAs derived from Cucumber green mottle mosaic virus-infected watermelon. Hortic Plant J. 2017;3:17–22.
Sun Y, Niu X, Fan M. Genome-wide identification of cucumber green mottle mosaic virus-responsive microRNAs in watermelon. Arch Virol. 2017;162:1–12.
Li X, Bi X, An M, Xia Z, Wu Y. iTRAQ-based proteomic analysis of watermelon fruits in response to Cucumber green mottle mosaic virus infection. Int J Mol Sci. 2020;21:2541.
Dominissini D, Moshitch-Moshkovitz S, Salmon-Divon M, Amariglio N, Rechavi G. Transcriptome-wide mapping of N6-methyladenosine by m6A-seq based on immunocapturing and massively parallel sequencing. Nat Protoc. 2013;8:176.
Ke S, Alemu EA, Mertens C, Gantman EC, Fak JJ, Mele A. Haripal BA majority of m6A residues are in the last exons, allowing the potential for 3’ UTR regulation. Gene Dev. 2015;29:2037–53.
Luo G, MacQueen A, Zheng G, Duan H, Dore L, Lu Z, et al. Unique features of the m6A methylome in Arabidopsis thaliana. Nat Commun. 2014;5:5630.
Zhou L, Tian S, Qin G. RNA methylomes reveal the m6A-mediated regulation of DNA demethylase gene SlDML2 in tomato fruit ripening. Genome Biol. 2019;20:1–23.
Visvanathan A, Somasundaram K. mRNA traffic control reviewed: N6-methyladenosine (m6A) takes the driver’s seat. BioEssays. 2018;40:1700093.
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:345–60.
Bruggeman Q, Garmier M, De Bont LD, Soubigou-Taconnat L, Mazubert C, Benhamed M, et al. The polyadenylation factor subunit CLEAVAGE AND POLYADENYLATION SPECIFICITY FACTOR30: a key factor of programmed cell death and a regulator of immunity in arabidopsis. Plant Physiol. 2014;165:732–46.
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:935–40.
Li Y, Deng C, Shang Q, Zhao X, Liu X, Zhou Q. Characterization of siRNAs derived from Cucumber green mottle mosaic virus in infected cucumber plants. Arch Virol. 2016;161:455–8.
Sun Y, Fan M, He Y. Transcriptome analysis of watermelon leaves reveals candidate genes responsive to Cucumber green mottle mosaic virus infection. Int J Mol Sci. 2019;20:610.
Huang Y, Li M, Wu P, Xu Z, Que F, Wang F, et al. Members of WRKY group III transcription factors are important in TYLCV defense signaling pathway in tomato (Solanum lycopersicum). BMC Genomics. 2016;17:788.
Kumar M, Busch W, Birke H, Kemmerling B, Nurnberger T, Schöffl F. Heat shock factors HsfB1 and HsfB2b are involved in the regulation of Pdf1.2 expression and pathogen resistance in Arabidopsis. Mol Plant. 2009;2:152–65.
Kuhlmann M, Horvay K, Stathmann A, Heinekamp T, Fischer U, Bottner S, et al. The alpha-helical D1 domain of the bZIP transcription factor BZI-1 interacts with the ankyrin-repeat protein ANK1, and is essential for BZI-1 function, both in auxin signaling and pathogen response. J Biol Chem. 2003;278:8786–94.
Noman A, Fahad S, Aqeel M, Ali U, Amanullah, Anwar S, et al. miRNAs: Major modulators for crop growth and development under abiotic stresses. Biotechnol Lett. 2017;39:685–700.
Biemelt S, Sonnewald U. Plant-microbe interactions to probe regulation of plant carbon metabolism. J Plant Physiol. 2006;163:307–18.
Llave C. Dynamic cross-talk between host primary metabolism and viruses during infections in plants. Curr Opin Virol. 2016;19:50–5.
Underwood W, Somerville SC. Phosphorylation is required for the pathogen defense function of the Arabidopsis PEN3 ABC transporter. Plant Signal Behav. 2017;12:e1379644.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 2015;31:2382–3.
Salmon-Divon M, Dvinge H, Tammoja K, Bertone P. PeakAnalyzer: genome-wide annotation of chromatin binding and modification loci. BMC Bioinformatics. 2010;11:415.
Cui X, Zhang L, Meng J, Rao M, Chen Y, Huang Y. MeTDiff: a novel differential RNA methylation analysis for MeRIP-Seq data. IEEE ACM T COMPUT BI. 2018;15:526–34.
Heinz S, Benner C, Spann N, Bertolino E, Lin Y, Laslo P, et al. Simple combinations of lineage-determining transcription factors primecis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38:576–89.
Nicol J, Helt G, Blanchard S, Raja A, Loraine A. The integrated genome browser: free software for distribution and exploration of genome-scale datasets. Bioinformatics. 2009;25:2730–1.
Trapnell C, Hendrickson D, Sauvageau M, Goff L, Rinn J, Pachter L. Differential analysis of gene regulation at transcriptre solution with RNA-seq. Nat Biotechnol. 2013;31:46.
Xu K, Yang Y, Feng GH, Sun BF, Chen JQ, Li YF, et al. Mettl3-mediated m6A regulates spermatogonial differentiation and meiosis initiation. Cell Res. 2017;27:1100–14.
Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative (T) method. Nat Protoc. 2008;3:1101.
Kong Q, Yuan J, Gao L, Zhao S, Jiang W, Huang Y, et al. Identification of suitable reference genes for gene expression normalization in qRT-PCR analysis in watermelon. PLoS ONE. 2014;9:e90612.
This research was funded by the National Key Research & Development Program of China (2018YFD0100703), Natural Science Foundation of Zhejiang Province (LQ18C150003), National Natural Science Foundation of China (31572145). The supporters had no role in study design, data collection, data analysis, data interpretation, the writing of the manuscript or decision to publish.
Ethics approval and consent to participate
The watermelon germplasms R1288 (Citrullus mucosospermus) and S1252 (Citrullus lanatus) used in our study were procured from the collection of Zhejiang Academy of Agricultural Sciences, Institute of Vegetables, China. In addition, samples collection complied with the institutional, national and international guidelines.
Consent for publication
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
High Pearson correlation coefficient (PCC) for the abundance of confident peaks between biological replicates.
Gene ontology (GO) analysis (GO: level 2) of the m6A-containing transcripts identified in m6A -seq in watermelon.
Pie charts depicting the fraction of DMPs in four non-overlapping transcript segments.
The phenotype of two watermelon germplasm (R1288 and S1252) plants infected CGMMV 20 d.
Summary of m6A -seq information for CGMMV-infected and control samples.
High-confidence m6A peaks identified for control samples and CGMMV-infected samples.
Gene ontology (GO) analysis of the biological process, molecular function, and cellular component for the m6A -containing transcripts identified in m6A -seq in watermelon.
Characteristics of differentially methylated peaks (DMPs) in CGMMV-infected and control samples.
The m6A modification level of DEGs involved in immune responses and the production of siRNA involved in RNA interference in this study.
DEGs identified in RNA-seq in this study.
The response patterns of watermelon m6A elements to CGMMV infection in RNA-seq data.
The differentially m6A methylated and expressed genes (DMEGs) identified in this study.
Primers used for m6A-immunoprecipitation(IP)-qPCR (m6A -IP-qPCR) and Real-time quantitative PCR (qPCR).
About this article
Cite this article
He, Y., Li, L., Yao, Y. et al. Transcriptome-wide N6-methyladenosine (m6A) methylation in watermelon under CGMMV infection. BMC Plant Biol 21, 516 (2021). https://doi.org/10.1186/s12870-021-03289-8