m6A methylation is a common feature of mRNAs in strawberry fruit
To investigate whether m6A methylation participates in modulating ripening of non-climacteric fruits, we performed m6A-seq [37] to characterize m6A methylomes on diploid woodland strawberry (Fragaria vesca) at three developmental stages, i.e. S6 (the growth stage 6, approximately 15 days post-anthesis (DPA)), RS1 (the ripening stage 1, 21 DPA), and RS3 (the ripening stage 3, 27 DPA) (Fig. 1a) [38]. The transition from S6 to RS1 represents the initiation of ripening, and that from RS1 to RS3 indicates the phase after ripening initiation. The m6A methylome libraries were prepared with three independently biological replicates and subjected to massively parallel sequencing according to the standard m6A-seq protocols [37]. High Pearson correlation coefficients were observed between biological replicates, indicating reliable repeatability (Additional file 1: Figure S1, S2). A total of 24–37 million raw reads were generated for each library (Additional file 2: Table S1), and this sequencing depth is comparative to that observed in mammal (11–24 million reads) [39], rice (23–47 million reads) [40], and tomato (20–30 million reads) [36]. After adaptor trimming and reads filtration, 24–37 million clean reads remained at each library, and almost 95% of these reads were uniquely aligned to the strawberry genome v2.0.a1, representing high mapping quality (Additional file 2: Table S1). The peak-calling algorithm was used to identify m6A peaks with an estimated false discovery rate (FDR) < 0.05 [37], and only those consistently detected in all three biological replicates, which we called confident m6A peaks, were used for subsequent analysis. We identified 9778, 10,853, and 10,095 confident m6A peaks within 8934, 8990, and 8374 gene transcripts, in fruit at S6, RS1, and RS3, respectively (Fig. 1b; Additional file 3: Table S2; Additional file 4: Table S3; Additional file 5: Table S4).
The m6A-seq results were validated by independent m6A immunoprecipitation followed by qPCR (m6A-IP-qPCR) analysis. Three m6A peak-containing transcripts, as well as three m6A peak-free transcripts, were randomly selected and examined (Additional file 1: Figure S3a). As expected, m6A enrichment was only observed in transcripts containing m6A peaks, but not in those without m6A peaks (Additional file 1: Figure S3b), indicating that our m6A-seq data were accurate and robust.
Based on the parallel RNA-seq analyses, we estimated that the transcriptome of diploid woodland strawberry contains 0.6–0.7 m6A peaks per actively expressed transcript, which shows FPKM (fragments per kilobase of transcript per million fragments mapped) ≥ 1 (Additional file 6: Table S5). These levels are comparable with those observed in Arabidopsis or tomato [36, 41], demonstrating that m6A modification is a common feature of mRNA in strawberry fruit. Most of the m6A-containing transcripts (> 86%) possess one m6A peak. Intriguingly, the percentage of transcripts harboring two or more m6A peaks increases dramatically when the strawberry fruit turn to ripen, which changes from 2.95% at S6 to 13.07% at RS1, and 13.63% at RS3 (Fig. 1c).
m6A distribution exhibits a dramatic change at the initiation stage of strawberry fruit ripening
We then evaluated the distribution of m6A peaks in the whole transcriptome of strawberry fruit. The transcript was divided into four non-overlapping segments: 5′ untranslated region (UTR), coding sequence (CDS), stop codon (100-nucleotide window centered on the stop codon), and 3′ UTR. As shown in Fig. 2a, m6A modifications in all three samples (fruit at S6, RS1, and RS3 stages) were highly enriched around the stop codon and within the 3′ UTR, but the proportion of m6A peaks within these regions declined markedly in the ripening process (from S6 to RS1 or RS3; indicated by green arrowheads). Surprisingly, a substantial increase in the proportion of m6A peaks was concurrently observed in the CDS region adjacent to the start codon in fruit at RS1 or RS3 compared to those at S6 (indicated by red arrowheads). This is different from that observed in tomato, which shows no dramatic changes in proportion of m6A peaks at the initiation stage of ripening [36].
The percentage of m6A peaks locating in the CDS region increased from 5.18 to 23.26% from S6 to RS1, but displayed no distinct changes from RS1 to RS3 (Fig. 2b). In comparison, the percentage of m6A peaks felling into the stop codon region and the 3′ UTR decreased from 39.95% and 53.81% to 34.22% and 38.99%, respectively, from S6 to RS1, with no obvious changes observed from RS1 to RS3 (Fig. 2b). After segment normalization by the relative fraction that each segment occupied in the transcriptome, the m6A enrichment in fruit at RS1 or RS3 consistently exhibits a preferential localization in the CDS region, besides enriched around the stop codon and within the 3′ UTR (Fig. 2c). Thus, the transcriptome-wide m6A distribution displays a dramatic change at the initiation stage of ripening, but not after ripening initiation.
Several high-confidence sequence motifs were identified within the m6A peaks (Additional file 1: Figure S4), by using the hypergeometric optimization of motif enrichment (HOMER; http://homer.ucsd.edu/homer/) [42]. The conserved RRACH consensus sequence observed in various organisms [1, 21, 41, 43], where R represents adenosine (A) or guanosine (G), underlined A indicates m6A, and H represents A, cytidine (C), or uridine (U), appears in the list, whereas the UGUA sequence motif found in Arabidopsis [22], tomato [36], and maize [25] was not identified. These data suggest the complexity of m6A modification among various species.
m6A methylation overall affects mRNA abundance during the ripening of strawberry fruit
To gain insight into the potential roles of m6A in the regulation of strawberry fruit ripening, we next sought for transcripts showing differential m6A peaks, which exhibit a fold change ≥ 1.5 and P value < 0.05 in m6A enrichment between the samples, by comparing the m6A methylomes. A total of 1608 hypermethylated m6A peaks and 865 hypomethylated m6A peaks, corresponding to 1398 and 790 transcripts, respectively, were identified in fruit at RS1 compared to those at S6 (Fig. 3a, c; Additional file 7: Table S6; Additional file 8: Table S7). By contrast, only 113 hypermethylated m6A peaks and 102 hypomethylated m6A peaks, which were distributed in 107 and 90 transcripts, respectively, were identified in fruit at RS3 compared to those at RS1 (Fig. 3b, d; Additional file 9: Table S8; Additional file 10: Table S9). These results confirmed that substantial changes in m6A methylome occurred at the initiation stage of ripening (from S6 to RS1), but not after ripening initiation (from RS1 to RS3). We did not observe an apparent increase in total m6A levels during strawberry ripening as revealed by LC-MS/MS (Additional file 1: Figure S5), which might be caused by the different m6A enrichment values between hypermethylated and hypomethylated peaks (Additional file 7: Table S6; Additional file 8: Table S7).
The 1608 hypermethylated m6A peaks were highly enriched (83.04%) in the CDS region, whereas the 865 hypomethylated m6A peaks were mainly distributed around the stop codon (40.66%) or within the 3′ UTR (58.59%) (Fig. 3e; Additional file 1: Figure S6). This is in accordance with the dynamic m6A distribution pattern (Fig. 2), showing that the percentage of m6A peaks locating in the CDS region increased sharply, while that around the stop codon or within the 3′ UTR declined, from S6 to RS1. Interestingly, of the 1608 hypermethylated m6A peaks, 1424 (88.56%) fell into the newly generated peaks, which represents ripening-specific peaks (Additional file 11: Table S10). For the differential m6A peaks identified after ripening initiation, both the hypermethylated and hypomethylated m6A peaks were highly enriched (over 90%) around the stop codon or within the 3′ UTR (Fig. 3e). To explore the biological significance of genes with dynamic m6A modification, we performed Gene Ontology (GO) analysis on genes containing differential, non-differential, and ripening-specific m6A peaks. In line with the progression of fruit development and ripening, genes harboring ripening-specific m6A peaks were mostly annotated to developmental pathways, including response to hormone stimulus, developmental process, histone modification, small molecular biosynthetic process, and protein transport (Additional file 1: Figure S7a). Similar enrichment was observed for genes covering differential m6A peaks at the initiation stage of ripening (from S6 to RS1) (Additional file 1: Figure S7b). In contrast, genes with non-differential m6A peaks were enriched in multiple biological processes in addition to developmental pathways (Additional file 1: Figure S7c). These results suggest that dynamic changes in m6A modification are responsible for those genes to exert their functions during fruit development and ripening.
Accumulating evidences have confirmed that m6A deposition affects mRNA abundance [1, 22, 41, 44]. To assess the potential correlation between m6A modification and mRNA abundance in strawberry fruit, we compared the list of transcripts harboring altered m6A methylation with the differentially expressed genes (fold change ≥ 1.5 and P value < 0.05) obtained from our parallel RNA-seq analyses (Additional file 12: Table S11; Additional file 13: Table S12). There were 3562 upregulated and 3583 downregulated genes, respectively, in RS1 compared to S6 fruits. Comparison of these differentially expressed genes with the transcripts containing differential m6A modification indicated that among the 1398 transcripts with hypermethylated m6A peaks in fruit at RS1 compared to those at S6, 440 and 147 transcripts displayed higher and lower expression levels, respectively (Fig. 3f; Additional file 14: Table S13). Accordingly, among the 790 transcripts showing hypomethylated m6A peaks in fruit at RS1 compared to those at S6, 267 transcripts showed higher expression levels, whereas only 51 transcripts exhibited lower expression levels (Fig. 3f; Additional file 15: Table S14). We subsequently extended this analysis to the whole transcriptome with all m6A-modified genes. Indeed, genes bearing hypermethylated or hypomethylated m6A peaks in fruit at RS1 compared to S6 predominantly exhibited higher expression (Fig. 3g, h). This trend did not exist for non-differential m6A-modified genes (Fig. 3g, h). Considering the distribution characteristics of hypermethylated m6A (highly enriched in the CDS region) and hypomethylated m6A (mainly distributed around the stop codon or within the 3′ UTR) (Fig. 3e), it was proposed that m6A deposition in the CDS region overall exhibits a positive effect on mRNA abundance, while m6A modification around the stop codon or within the 3′ UTR is generally negatively correlated with the abundance of the transcripts. Plotting the genes harboring differential m6A peaks in different segments against their expression levels confirmed the complex relationship between m6A modification and gene expression levels (Fig. 3i). The negative role of m6A modification on mRNA abundance was also observed in transcripts with differential m6A peaks identified after ripening initiation (Fig. 3g, h; Additional file 16: Table S15; Additional file 17: Table S16), which were mainly distributed around the stop codon or within the 3′ UTR (Fig. 3e).
Notably, hundreds of ripening-induced and ripening-repressed genes, which display significantly higher or lower expression in RS1 compared to S6 (Additional file 12: Table S11), exhibit differential m6A modification (Additional file 14: Table S13; Additional file 15: Table S14). GO analysis revealed that these genes were enriched in processes such as multicellular organismal development, developmental process, nucleocytoplasmic transport, and anatomical structure development (Additional file 1: Figure S7d), implicating the involvement of m6A methylation in the regulation of strawberry fruit ripening.
Genes in ABA biosynthesis and signaling pathway exhibit differential m6A methylation upon ripening initiation
The plant hormone ABA has been elucidated to be essential for strawberry fruit ripening [32]. Two core ABA signal transduction pathways, the “ABA-PYR/PYL-PP2C-SnRK2-AREB/ABF” pathway [45, 46] and the “ABA-ABAR-WRKY40-ABI5” pathway [47], have been proposed in Arabidopsis, respectively (Fig. 4a). In the m6A-seq analysis, we found that transcripts of key genes in ABA biosynthesis and signaling pathway, including 9-cis-epoxycarotenoid dioxygenase 5 (NCED5) [38], putative ABA receptor (ABAR) [34], and ABA-responsive element-binding protein 1 (AREB1) [46], exhibit hypermethylation in the CDS region at the initiation stage of strawberry fruit ripening (from S6 to RS1) (Fig. 4b, c). NCED5 encodes the rate-limiting enzyme for ABA biosynthesis, while ABAR and AREB1 encode an ABA receptor and an ABA-responsive element, respectively. All three genes have been proposed to regulate strawberry fruit ripening [34, 38, 48]. The differential m6A modifications were confirmed by m6A-IP-qPCR analysis (Fig. 4d). The transcript levels of NCED5 and AREB1, but not ABAR, increased significantly in fruit at RS1 compared to those at S6 as revealed by both RNA-seq (Fig. 4e) and quantitative RT-PCR analyses (Fig. 4f), indicating a positive correlation between m6A depositions and mRNA abundances. It is noteworthy that similar results of m6A enrichment as well as transcript levels of these three genes were found upon ripening initiation of the octoploid cultivated strawberry (Fragaria × ananassa) (Additional file 1: Figure S8a, b).
To better understand how m6A methylation affects the abundance of the transcripts, we examined the mRNA stability of NCED5, AREB1, and ABAR by monitoring the degradation rate of mRNAs in the presence of transcription inhibitor actinomycin D. The intact CDS sequences of the three genes and their mutated forms in which the specific m6A sites identified in m6A-seq and validated by SELECT analysis (Additional file 1: Figure S9; Additional file 18: Table S17), a method for detection of single m6A locus [49], were mutated from A to guanosine (G) were separately inserted into the vector for transient expression in the Nicotiana benthamiana leaves (Fig. 4 g). As shown in Fig. 4 h, the mRNAs of NCED5, AREB1, and ABAR degraded obviously after actinomycin D treatment. The degradation rate was substantially increased for mRNAs of NCED5 and AREB1, but not ABAR, in the mutated form (Fig. 4h), concomitant with significantly diminished m6A depositions (Fig. 4i), indicating that site-specific m6A modification stabilizes mRNAs of NCED5 and AREB1.
The m6A modification in the CDS region has been demonstrated to affect translation efficiency beyond mRNA stability [3]. We next investigated whether m6A methylation modulates translation efficiency of NCED5, AREB1, and ABAR, which was determined by calculating the abundance ratio of mRNA in the polysomal RNA versus the total RNA [50] (Fig. 4j). The translation efficiency of ABAR rather than NCED5 and AREB1 was significantly decreased when the specific m6A site was mutated (Fig. 4k), demonstrating that m6A methylation facilitates translation of ABAR. Collectively, these data suggest that critical genes in ABA pathway undergo m6A-mediated post-transcriptional regulation, which promotes mRNA stability or facilitates translation.
Characterizations of m6A methyltransferases in strawberry fruit
Having observed the changes in m6A methylation in a large number of transcripts (Fig. 3), including those of ABA biosynthesis and signaling genes (Fig. 4), during the ripening of strawberry fruit, we next explored the mechanistic basis. We speculate that the ripening-specific hypermethylation in the CDS region is regulated by mRNA m6A methyltransferases. In mammals, METTL3 and METTL14 form a stable heterodimer wherein METTL3 serves as the m6A catalytic subunit and METTL14 facilitates RNA binding [12, 13] (Fig. 5a). The m6A methyltransferase complexes are conserved between mammals and plants [51], and the Arabidopsis mRNA adenosine methyltransferase (MTA) and MTB appear to be the homologs of METTL3 and METTL14, respectively [52, 53]. We identified the single homologs of Arabidopsis MTA and MTB in the genome of the diploid woodland strawberry. Phylogenetic analysis indicated that MTA and MTB exhibit high similarity among plant species and are evolutionarily conserved with mammal METTL3 and METTL14 (Fig. 5b). Both MTA and MTB in strawberry contain the highly conserved MT-A70 domain (Additional file 1: Figure S10), which displays extremely similar protein sequences to those observed in mammals and Arabidopsis (Additional file 1: Figure S11).
Transcriptome analysis showed that MTA and MTB increased significantly at the initiation stage of ripening (from S6 to RS1) (Fig. 5c). The homolog genes of MTA and MTB in the octoploid cultivated strawberry also exhibited increased expression from white (Wt) stage to initial red (IR) stage (Fig. 5d, e), suggesting that the two methyltransferases may play important roles in modulating strawberry fruit ripening.
To explore the possibility that MTA and MTB in strawberry function in the form of heterodimer as METTL3 and METTL14 in mammals [12, 13], we subsequently analyzed the interactions between MTA and MTB using the yeast two-hybrid (Y2H) system. As shown in Fig. 5f, the yeast cells co-expressing AD-MTA and BD-MTB, but not the negative controls, displayed normal growth on the selective SD/-Leu-Trp-His (-LWH) and SD/-Leu-Trp-His-Ade (-LWHA) solid medium and turn to blue with the addition of X-α-gal, indicating that MTA interacts with MTB. The interactions between MTA and MTB were further verified by the split luciferase complementation imaging (LCI) assay, in which the luciferase activity was detected when MTA-nLUC and cLUC-MTB were co-expressed in N. benthamiana leaves (Fig. 5g). It should be noted that, compared with the MT-A70 domain of MTB, the full-length MTB protein exhibit relatively weaker combining capacity with MTA (Fig. 5g). Subcellular localization analysis by transiently expressing MTA-mCherry and MTB-eGFP fusion proteins in N. benthamiana leaves showed that MTA is present in both the nucleus and cytoplasm, while the MTB protein is specifically localized in the nucleus (Fig. 5h). Interestingly, when MTA-mCherry was co-expressed with MTB-eGFP, the two proteins tend to colocalize in the nucleus (Fig. 5i).
m6A methyltransferases positively regulate strawberry fruit ripening
We subsequently examined the function of MTA and MTB in strawberry fruit ripening. The RNA interference (RNAi) and overexpression constructs of MTA or MTB under the control of a 35S cauliflower mosaic virus promoter were agroinfiltrated into the octoploid strawberry fruit according to the reported procedures [54]. By comparing the fruits of RNAi with the control, we found that suppression of either MTA or MTB delayed fruit ripening (Fig. 6a). A visible color change was observed at the seventh day after agroinfiltration in the control, whereas the MTA or MTB RNAi fruits were almost green at this stage (Fig. 6a). Conversely, overexpression of MTA or MTB accelerated fruit ripening (Fig. 6a). Gene expression analysis indicated that MTA and MTB were successfully silenced in the RNAi fruits while enhanced in the overexpressed fruits (Fig. 6b). The ripening genes chalcone synthase (CHS) and polygalacturonase 1 (PG1) displayed a noticeable decrease in the MTA or MTB RNAi fruits, but were dramatically enhanced in the overexpressed fruits (Fig. 6b). These results suggest that MTA and MTB are necessary for normal fruit ripening of strawberry.
We next explored how the ripening of strawberry was regulated by MTA, the core component of the methyltransferase complexes with m6A catalytic activity. The MTA RNAi fruits showed lower m6A levels, while the MTA-overexpressed fruits displayed higher m6A levels, than the control as revealed by LC-MS/MS (Fig. 6c). RNA immunoprecipitation (RIP) analysis indicated that the MTA directly binds to the transcripts of the ABA biosynthesis or signaling genes NCED5, ABAR, and AREB1 (Fig. 6d). Accordingly, the m6A enrichments in the transcripts of NCED5, ABAR, and AREB1 (Fig. 6e) and the corresponding mRNA levels (Fig. 6f) were significantly decreased in the RNAi fruits, but increased in the overexpressed fruits, compared to the control. The mRNA stability assay demonstrated that mRNA of NCED5 and AREB1, but not that of ABAR, degraded more quickly in the MTA RNAi fruits, but more slowly in the MTA-overexpressed fruits, compared to the control (Fig. 6g; Additional file 1: Figure S12a). By contrast, MTA overexpression or repression markedly enhanced or decreased the translation efficiency of the ABAR mRNA, but showed no significant effects on that of NCED5 and AREB1 (Fig. 6h). Together, these results suggest that MTA may regulate strawberry fruit ripening by targeting genes in the ABA pathway, leading to the increase in mRNA stability or translation efficiency of these genes.
MTA-mediated m6A methylation modulates genes encoding translation initiation factors and elongation factors
Due to the critical role of ABA in the regulation of fruit ripening of strawberry, we evaluated whether MTA affects translation efficiency of other genes in the ABA signaling pathway. We found that genes, such as WRKY DNA-binding protein 40 (WRKY40), exhibited significant changes in translation efficiency when the MTA was silenced or overexpressed (Additional file 1: Figure S13). This could not be reasonably explained by m6A deposition because the transcripts of these genes are not m6A-modified according to our m6A-seq datasets. We speculate that MTA may regulate translation efficiency of numerous transcripts beyond direct m6A installation. As expected, MTA repression or overexpression also altered the translation efficiency of a number of ripening-related genes without m6A modification, such as PG1 relevant to firmness and dihydroflavonol 4-reductase (DFR) associated with anthocyanin biosynthesis (Additional file 1: Figure S13).
To investigated the possible mechanisms, we mined our m6A-seq and RNA-seq data and found that the transcripts of genes encoding translation initiation factors (EIF2, EIF2B, EIF3A, and EIF3C) and elongation factors (EF1A), which play pivotal roles in facilitating protein synthesis by promoting the initiation and elongation of mRNA translation, respectively, exhibited m6A hypermethylation in the CDS region upon ripening initiation of strawberry fruit (Fig. 7a, b), concomitant with an increase in the transcript levels of these genes (Fig. 7c). Similar results were observed in the ripening process of the octoploid strawberry fruit (Additional file 1: Figure S8c, d). RIP analysis showed direct interactions between MTA and the transcripts of EIF2, EIF2B, EIF3A, and EIF3C, and EF1A (Fig. 7d). More importantly, the m6A enrichment in these gene transcripts (Fig. 7e, g) and the mRNA abundance (Fig. 7f, h) declined or increased noticeably when MTA was silenced or overexpressed. The mRNA degradation rates were obviously increased in the MTA RNAi fruits, but reduced in the MTA-overexpressed fruits, compared to the control (Fig. 7i; Additional file 1: Figure 12b), implying that MTA-mediated m6A modification stabilizes mRNAs of these genes. We propose that, in addition to direct modulation of translation efficiency via m6A installation in the target transcripts, MTA may indirectly modulate translation efficiency via m6A-mediated regulation of the translation initiation factors or elongation factors.
MTA is indispensable for global m6A mRNA methylation
To gain a more detailed insight into the MTA-mediated m6A methylation, we performed m6A-seq on fruits of the MTA RNAi and the control. A total of 8693 and 8704 confident m6A peaks within 7933 and 5180 gene transcripts in the MTA RNAi fruits and the control, respectively, were identified (Fig. 8a; Additional file 19: Table S18; Additional file 20: Table S19). Sequence preference analysis indicated that the conserved RRACH motif was identified among these m6A peaks, which harbour high Pearson correlation coefficients between biological replicates, implicating high reliability (Additional file 1: Figure S14). Compared with the control, we identified 2520 hypomethylated m6A peaks covering 1128 gene transcripts and 1510 hypermethylated m6A peaks within 1239 gene transcripts in the MTA RNAi fruits (Fig. 8b, c; Additional file 21: Table S20; Additional file 22: Table S21), indicating an overall m6A hypomethylation when the MTA was silenced. The hypomethylated m6A peaks were mainly distributed in the CDS region (76.19%), while the hypermethylated m6A peaks were highly enriched around the stop codon (31.26%) or within the 3′ UTR (61.46%) (Fig. 8d). This is similar to the distribution of differential m6A peaks identified during fruit ripening (Fig. 3e).
RNA-seq analysis indicated that there were 4488 downregulated and 4754 upregulated genes (fold change ≥ 1.5 and P value < 0.05), respectively, in the MTA RNAi fruits (Additional file 23: Table S22). By analyzing the overall expression of m6A-modified transcripts based on our RNA-seq data, we found that genes harboring hypomethylated m6A modification in the MTA RNAi fruits were more likely downregulated, and this trend did not exist for non-differential m6A-modified genes (Fig. 8e, f), suggesting a positive correlation between MTA-mediated m6A methylation and mRNA abundance. GO analysis revealed that genes containing hypomethylated m6A modification were significantly annotated to biological processes including response to hormone stimulus, developmental process, multicellular organismal development, cell wall organization, and secondary metabolic process (Fig. 8g), which are prerequisites for fruit development and ripening. Importantly, we found that a number of ripening-related genes, including those involved in anthocyanin biosynthesis, sugar metabolism, cell wall organization and degradation, and hormone signaling, displayed hypomethylation in the MTA RNAi fruits (Additional file 24: Table S23). Together, these results demonstrated that MTA is indispensable for global m6A mRNA methylation, which ensures the proper progression of fruit ripening of strawberry by targeting ripening-related genes among various biological processes.