N⁶-methyladenosine RNA modification regulates strawberry fruit ripening in an ABA-dependent manner

Abstract

Background

Epigenetic mark such as DNA methylation plays pivotal roles in regulating ripening of both climacteric and non-climacteric fruits. However, it remains unclear whether mRNA m⁶A methylation, which has been shown to regulate ripening of the tomato, a typical climacteric fruit, is functionally conserved for ripening control among different types of fruits.

Results

Here we show that m⁶A methylation displays a dramatic change at ripening onset of strawberry, a classical non-climacteric fruit. The m⁶A modification in coding sequence (CDS) regions appears to be ripening-specific and tends to stabilize the mRNAs, whereas m⁶A around the stop codons and within the 3′ untranslated regions is generally negatively correlated with the abundance of associated mRNAs. We identified thousands of transcripts with m⁶A hypermethylation in the CDS regions, including those of NCED5, ABAR, and AREB1 in the abscisic acid (ABA) biosynthesis and signaling pathway. We demonstrate that the methyltransferases MTA and MTB are indispensable for normal ripening of strawberry fruit, and MTA-mediated m⁶A modification promotes mRNA stability of NCED5 and AREB1, while facilitating translation of ABAR.

Conclusion

Our findings uncover that m⁶A methylation regulates ripening of the non-climacteric strawberry fruit by targeting the ABA pathway, which is distinct from that in the climacteric tomato fruit.

Introduction

As the most prevalent chemical modification in eukaryotic messenger RNAs (mRNAs), N⁶-methyladenosine (m⁶A) has been demonstrated to functionally modulate multiple biological processes through interfering mRNA metabolism [1,2,3,4]. In mammals, m⁶A methylation has been unveiled to play critical roles in regulating various physiological and pathological processes, such as embryonic and post-embryonic development, cell circadian rhythms, and cancer stem cell proliferation [5,6,7,8,9,10]. The m⁶A marks in mammals are dominantly installed by the methyltransferase complex composed of the stable catalytic core that is formed by methyltransferase-like 3 (METTL3) and METTL14 [11,12,13], the Wilms tumor 1-associating protein (WTAP) [14], and other concomitant functional elements [15,16,17]. The removal of m⁶A marks is executed by two m⁶A demethylases, fat mass and obesity-associated protein (FTO) and alkylated DNA repair protein AlkB homolog 5 (ALKBH5) [18, 19]. The m⁶A modification is recognized by the reader proteins, such as YTH-domain family proteins and specific RNA binding proteins (RBPs), which mediate the downstream effects of m⁶A methylation [4]. In plants, the m⁶A methylation machineries have been characterized in the model plant Arabidopsis thaliana to modulate development processes such as shoot stem cell proliferation, trichome branching, and floral transition [20,21,22,23]. Moreover, m⁶A has been demonstrated to play pivotal roles in mediating sporogenesis in rice [24] and regulating stress responses in maize [25]. By contrast, the m⁶A methylation machineries as well as the characteristics and functions of m⁶A in regulating physiological processes of horticultural crops remain largely unknown.

Fleshy fruits, which are enriched with nutrients, such as flavor compounds, fiber, vitamins and antioxidants, represent one of the commercially valuable structures of horticultural crops. As an important component of diets, fleshy fruits are indispensable for human health [26]. The ripening of fleshy fruits, which is characterized by dramatic changes in color, texture, flavor and aroma compounds [27], is a complex, genetically programmed process that impacts fruit nutritional quality and shelf life. Fruit ripening is regulated by both environmental and internal cues, including light, temperature, phytohormones, and developmental genes [28, 29]. Based on the different ripening mechanisms, fruits are classified into two groups: climacteric (e.g., tomato, apple, banana, and avocado) and non-climacteric (e.g., strawberry, grape, and citrus) [30]. Phytohormone ethylene is essential for the ripening of climacteric fruits [31, 32], and substantial insights have been made toward ethylene biosynthesis, ethylene perception and signal transduction, and downstream gene regulation [33]. In comparison, the ripening of non-climacteric fruits is thought to be abscisic acid (ABA)-dependent [32, 34], although the regulation of ABA pathway is poorly understood. A comprehensive understanding of the common regulatory mechanisms underlying ripening in climacteric and non-climacteric species has great potential for improving fruit quality and maintaining shelf-life.

Recently, it was shown that epigenetic marks, including DNA methylation and histone posttranslational modifications, play critical roles in the regulation of fruit ripening [35]. In a previous study, we uncovered that mRNA m⁶A methylation, which is considered as an mRNA “epitranscriptome”, exhibits dynamic changes during fruit ripening of tomato, a typical climacteric fruit [36]. Mutation of SlALKBH2, the m⁶A RNA demethylase gene, delays fruit ripening [36], indicating that m⁶A modification participates in the ripening control of climacteric fruit tomato. However, whether m⁶A is evolutionarily conserved among different types of fruits has not been defined. Moreover, the regulatory function of m⁶A in ripening of non-climacteric fruits remains elusive.

In the present study, we performed transcriptome-wide m⁶A methylation (m⁶A methylome) in strawberry, a classical non-climacteric fruit, and revealed that m⁶A represents a prevalent modification in the mRNAs of strawberry fruit. Compared to the dynamic changes in m⁶A modification around the stop codons or within the 3′ untranslated regions during the ripening of tomato fruit, a specific enrichment of m⁶A in the coding sequence (CDS) region, which tends to be positively correlated with the abundance of the transcripts, was observed in the ripe strawberry fruit. We demonstrated that, mediated by the methyltransferases MTA and MTB, m⁶A modification enhances mRNA stability or promotes translation efficiency of genes in the ABA biosynthesis and signaling pathway, thereby facilitating the ripening of strawberry fruit. Our study uncovers the regulatory effects of m⁶A methylation on non-climacteric strawberry fruit ripening and reveals a direct role for m⁶A methylation in the regulation of key elements in the ABA pathway.

Results

m⁶A methylation is a common feature of mRNAs in strawberry fruit

To investigate whether m⁶A methylation participates in modulating ripening of non-climacteric fruits, we performed m⁶A-seq [37] to characterize m⁶A 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 m⁶A methylome libraries were prepared with three independently biological replicates and subjected to massively parallel sequencing according to the standard m⁶A-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 m⁶A 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 m⁶A peaks, were used for subsequent analysis. We identified 9778, 10,853, and 10,095 confident m⁶A 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).

Fig. 1

Transcriptome-wide m⁶A methylomes in strawberry fruit. a The representative photographs of fruit at different developmental stages. S6, the growth stage 6; RS1, the ripening stage 1; RS3, the ripening stage 3. Scale bar = 1 cm. b Venn diagrams depicting the overlap of m⁶A peaks from three independent m⁶A-seq experiments on fruit at the three developmental stages. Rep, replicate. c Percentage of the m⁶A-containing transcripts containing various m⁶A peak numbers among samples

The m⁶A-seq results were validated by independent m⁶A immunoprecipitation followed by qPCR (m⁶A-IP-qPCR) analysis. Three m⁶A peak-containing transcripts, as well as three m⁶A peak-free transcripts, were randomly selected and examined (Additional file 1: Figure S3a). As expected, m⁶A enrichment was only observed in transcripts containing m⁶A peaks, but not in those without m⁶A peaks (Additional file 1: Figure S3b), indicating that our m⁶A-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 m⁶A 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 m⁶A modification is a common feature of mRNA in strawberry fruit. Most of the m⁶A-containing transcripts (> 86%) possess one m⁶A peak. Intriguingly, the percentage of transcripts harboring two or more m⁶A 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).

m⁶A distribution exhibits a dramatic change at the initiation stage of strawberry fruit ripening

We then evaluated the distribution of m⁶A 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, m⁶A 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 m⁶A 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 m⁶A 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 m⁶A peaks at the initiation stage of ripening [36].

Fig. 2

Dynamic distribution of m⁶A during strawberry fruit ripening. a Metagenomic profiles of m⁶A peak distribution along transcripts. UTR, untranslated region; CDS, coding sequence. The red arrows indicate the changes in m⁶A distribution in the CDS region adjacent to the start codon in the ripening process, while the green arrows show the changes in m⁶A distribution around the stop codon or within the 3′ UTR. S6, the growth stage 6; RS1, the ripening stage 1; RS3, the ripening stage 3. b, c The percentage (b) and relative enrichment (c) of m⁶A peaks in four non-overlapping transcript segments

The percentage of m⁶A 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 m⁶A 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 m⁶A 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 m⁶A 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 m⁶A 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 m⁶A, 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 m⁶A modification among various species.

m⁶A methylation overall affects mRNA abundance during the ripening of strawberry fruit

To gain insight into the potential roles of m⁶A in the regulation of strawberry fruit ripening, we next sought for transcripts showing differential m⁶A peaks, which exhibit a fold change ≥ 1.5 and P value < 0.05 in m⁶A enrichment between the samples, by comparing the m⁶A methylomes. A total of 1608 hypermethylated m⁶A peaks and 865 hypomethylated m⁶A 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 m⁶A peaks and 102 hypomethylated m⁶A 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 m⁶A 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 m⁶A levels during strawberry ripening as revealed by LC-MS/MS (Additional file 1: Figure S5), which might be caused by the different m⁶A enrichment values between hypermethylated and hypomethylated peaks (Additional file 7: Table S6; Additional file 8: Table S7).

Fig. 3

Correlation between m⁶A modification and mRNA abundance in strawberry fruit. a Volcano plot showing hypermethylated (red) and hypomethylated (blue) m⁶A peaks in fruit at RS1 stage compared to those at S6 stage. S6, the growth stage 6; RS1, the ripening stage 1. b Volcano plot depicting hypermethylated (red) and hypomethylated (blue) m⁶A peaks in fruit at RS3 stage compared to those at RS1 stage. RS3, the ripening stage 3. c, d Heat map of m⁶A enrichment ratios of differential m⁶A peaks shown in a (c) and b (d). e The distribution characteristics of the differential m⁶A peaks shown in a and b. UTR, untranslated region; CDS, coding sequence. f Volcano plots displaying the expression ratios of transcripts containing differential m⁶A peaks shown in a and b. Transcripts with significantly higher and lower mRNA levels (fold change ≥ 1.5; P value < 0.05) in fruit at RS1 stage compared to those at S6 stage, or in fruit at RS3 stage compared to those at RS1 stage, are highlighted in red and blue, respectively. g, h Cumulative distribution of gene expression changes (g) and box plot of gene expression ratios (h) between fruit at RS1 and S6 or fruit at RS3 and RS1. The m⁶A-modified transcripts were divided into three categories, including hypermethylated, hypomethylated, and non-differential transcripts. i Box plot showing gene expression ratios based on the distributions of differential m⁶A peaks between fruit at RS1 and S6 stages. Asterisks indicate significant differences (***P < 0.001; Wilcoxon test)

The 1608 hypermethylated m⁶A peaks were highly enriched (83.04%) in the CDS region, whereas the 865 hypomethylated m⁶A 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 m⁶A distribution pattern (Fig. 2), showing that the percentage of m⁶A 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 m⁶A peaks, 1424 (88.56%) fell into the newly generated peaks, which represents ripening-specific peaks (Additional file 11: Table S10). For the differential m⁶A peaks identified after ripening initiation, both the hypermethylated and hypomethylated m⁶A 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 m⁶A modification, we performed Gene Ontology (GO) analysis on genes containing differential, non-differential, and ripening-specific m⁶A peaks. In line with the progression of fruit development and ripening, genes harboring ripening-specific m⁶A 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 m⁶A peaks at the initiation stage of ripening (from S6 to RS1) (Additional file 1: Figure S7b). In contrast, genes with non-differential m⁶A peaks were enriched in multiple biological processes in addition to developmental pathways (Additional file 1: Figure S7c). These results suggest that dynamic changes in m⁶A modification are responsible for those genes to exert their functions during fruit development and ripening.

Accumulating evidences have confirmed that m⁶A deposition affects mRNA abundance [1, 22, 41, 44]. To assess the potential correlation between m⁶A modification and mRNA abundance in strawberry fruit, we compared the list of transcripts harboring altered m⁶A 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 m⁶A modification indicated that among the 1398 transcripts with hypermethylated m⁶A 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 m⁶A 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 m⁶A-modified genes. Indeed, genes bearing hypermethylated or hypomethylated m⁶A peaks in fruit at RS1 compared to S6 predominantly exhibited higher expression (Fig. 3g, h). This trend did not exist for non-differential m⁶A-modified genes (Fig. 3g, h). Considering the distribution characteristics of hypermethylated m⁶A (highly enriched in the CDS region) and hypomethylated m⁶A (mainly distributed around the stop codon or within the 3′ UTR) (Fig. 3e), it was proposed that m⁶A deposition in the CDS region overall exhibits a positive effect on mRNA abundance, while m⁶A 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 m⁶A peaks in different segments against their expression levels confirmed the complex relationship between m⁶A modification and gene expression levels (Fig. 3i). The negative role of m⁶A modification on mRNA abundance was also observed in transcripts with differential m⁶A 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 m⁶A 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 m⁶A methylation in the regulation of strawberry fruit ripening.

Genes in ABA biosynthesis and signaling pathway exhibit differential m⁶A 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 m⁶A-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 m⁶A modifications were confirmed by m⁶A-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 m⁶A depositions and mRNA abundances. It is noteworthy that similar results of m⁶A 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).

Fig. 4

m⁶A modification facilitates mRNA stability or translation of genes in the ABA pathway. a A brief model of the two core ABA signaling pathways in plants. b Integrated Genome Browser (IGB) tracks showing the distribution of m⁶A reads in transcripts of the 9-cis-epoxycarotenoid dioxygenase 5 (NCED5), putative ABA receptor (ABAR), and ABA-responsive element-binding protein 1 (AREB1). The hypermethylated m⁶A peaks (fold change ≥ 1.5; P value < 0.05) in fruit at the RS1 stage compared to those at the S6 stage are indicated by shadow boxes. S6, the growth stage 6; RS1, the ripening stage 1. Rep, replicate. c m⁶A enrichment for NCED5, ABAR, and AREB1 from m⁶A-seq data. d Validations of the m⁶A enrichment by m⁶A-immunoprecipitation (IP)-qPCR. e, f Transcript levels of NCED5, ABAR, and AREB1 determined by RNA-seq (e) and quantitative RT-PCR (f). FPKM, fragments per kilobase of exon per million mapped fragments. The ACTIN gene served as an internal control in f. g Schematic diagram of the expression cassettes used for mRNA stability assay. The intact (WT) or mutated (MU) CDS of NCED5, ABAR, and AREB1 were separately cloned into the pCambia2300 vector driven by the CaMV 35S promoter. The specific m⁶A sites (highlighted in red) identified in m⁶A-seq and validated by SELECT were mutated from adenosine (A) to guanine (G) using site-directed mutagenesis kit. h Determination of mRNA stability for NCED5, ABAR, and AREB1. The intact (WT) or mutated (MU) CDS were expressed in Nicotiana benthamiana leaves. After actinomycin D treatment at an indicated time point, the total RNAs were extracted and submitted to quantitative RT-PCR assay with the N. benthamiana ACTIN gene serving as an internal control. i m⁶A-IP-qPCR assay showing the relative m⁶A enrichment in the intact or mutated transcripts. j Brief workflow for the analysis of translation efficiency. k Translation efficiency of NCED5, ABAR, and AREB1. Translation efficiency was expressed as the abundance ratio of mRNA in the polysomal RNA versus the total RNA. Data are presented as mean ± standard deviation (n = 3). Asterisks indicate significant differences (*P < 0.05, **P < 0.01, ***P < 0.001; Student’s t test)

To better understand how m⁶A 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 m⁶A sites identified in m⁶A-seq and validated by SELECT analysis (Additional file 1: Figure S9; Additional file 18: Table S17), a method for detection of single m⁶A 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 m⁶A depositions (Fig. 4i), indicating that site-specific m⁶A modification stabilizes mRNAs of NCED5 and AREB1.

The m⁶A modification in the CDS region has been demonstrated to affect translation efficiency beyond mRNA stability [3]. We next investigated whether m⁶A 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 m⁶A site was mutated (Fig. 4k), demonstrating that m⁶A methylation facilitates translation of ABAR. Collectively, these data suggest that critical genes in ABA pathway undergo m⁶A-mediated post-transcriptional regulation, which promotes mRNA stability or facilitates translation.

Characterizations of m⁶A methyltransferases in strawberry fruit

Having observed the changes in m⁶A 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 m⁶A methyltransferases. In mammals, METTL3 and METTL14 form a stable heterodimer wherein METTL3 serves as the m⁶A catalytic subunit and METTL14 facilitates RNA binding [12, 13] (Fig. 5a). The m⁶A 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).

Fig. 5

The m⁶A methyltransferase MTA interacts with MTB in strawberry. a The working model for m⁶A installations mediated by the methyltransferase complex in mammals. The m⁶A methyltransferases METTL3 and METTL14 interact and function as the stable catalytic core for internal m⁶A installations in the form of heterodimer. b Phylogenetic analysis of eukaryotic m⁶A methyltransferases. The phylogenetic tree was generated by MEGA (version 5.2). Bootstrap values from 1000 replications for each branch are presented. Species names are abbreviated as follows: Hs, Homo sapiens; Ms, Mus musculus; At, Arabidopsis thaliana; Os, Oryza sativa; Zm, Zea mays; Sl, Solanum lycopersicum; Nb, Nicotiana benthamiana; Fa, Fragaria × ananassa; Fve, Fragaria vesca. c Transcript levels of the m⁶A methyltransferase genes MTA and MTB in diploid woodland strawberry at different developmental stages revealed by RNA-seq. S6, the growth stage 6; RS1, the ripening stage 1; RS3, the ripening stage 3. d Representative images of octoploid cultivated strawberry fruit at various developmental stages. SG, small green; Wt, white; IR, initial red; FR, full red. Scale bar = 1 cm. e Transcript levels of MTA and MTB in octoploid strawberry fruit determined by quantitative RT-PCR. The ACTIN gene was used as an internal control. Data are presented as mean ± standard deviation (n = 3). Asterisks indicate significant differences (*P < 0.05, **P < 0.01, ***P < 0.001; Student’s t test). f Y2H assay revealing the interactions between MTA and MTB. The MTA fused with the activation domain (AD) of GAL4 (AD-MTA) and the MTB fused with the binding domain (BD) of GAL4 (BD-MTB) were co-expressed in yeast. The transformants were grew on SD/-Leu/-Trp (-LW), and further selected on SD/-Leu/-Trp/-His (-LWH) and SD/-Leu/-Trp/-His/-Ade (-LWHA) with or without X-α-gal. g LCI assay revealing the interactions between MTA and MTB. The MTA fused with the N-terminus of LUC (MTA-nLUC) was co-expressed with the MTB or its MT-A70 domain fused with the C-terminus of LUC (cLUC-MTB or cLUC-MTB^D) in N. benthamiana leaves. h, i Subcellular localization (h) and colocalization (i) of MTA and MTB. The MTA-mCherry or/and MTB-eGFP fusion proteins were transiently expressed into N. benthamiana leaves. The N. benthamiana leaves expressing eGFP or/and mCherry were used as the negative control. Scale bar = 10 μm

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).

m⁶A 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.

Fig. 6

m⁶A methyltransferases positively regulate the ripening of strawberry fruit. a Ripening phenotypes of MTA/MTB RNA interference (RNAi-MTA/RNAi-MTB) and overexpression (OE-MTA/OE-MTB) fruits. Strawberry fruit agroinfiltrated with empty plasmid were used as controls. The experiments were performed with more than three highly reproducible biological replicates and the representative results are presented. Scale bar = 1 cm. b Transcript levels of MTA and MTB, as well as the two important ripening genes CHS and PG1, in the RNAi and overexpression fruits determined by quantitative RT-PCR. The ACTIN gene served as an internal control. c LC-MS/MS assay revealing the changes in global m⁶A methylation levels in the MTA-silenced (RNAi-MTA) or MTA-overexpressed (OE-MTA) fruits. d RNA immunoprecipitation (RIP) assay revealing the binding of MTA protein to the transcripts of NCED5, ABAR, and AREB1. The protein-RNA complexes were extracted from strawberry fruit expressing the MTA-eGFP fusion protein and subjected to immunoprecipitation with anti-GFP monoclonal antibody or mouse IgG (negative control). e–h The changes in relative m⁶A enrichment (e), gene expression (f), mRNA stability (g), and translation efficiency (h) of NCED5, ABAR, and AREB1 in the MTA-silenced (RNAi-MTA) or MTA-overexpressed (OE-MTA) fruits. The relative m⁶A enrichment and gene expression were determined by m⁶A-IP-qPCR and quantitative RT-PCR, respectively. For mRNA stability assay, the total RNAs were extracted after actinomycin D treatment at an indicated time point and submitted to quantitative RT-PCR assay. Translation efficiency was expressed as the abundance ratio of mRNA in the polysomal RNA versus the total RNA. Data are presented as mean ± standard deviation (n = 3). Asterisks indicate significant differences (*P < 0.05, **P < 0.01, ***P < 0.001; Student’s t test). NS, no significance

We next explored how the ripening of strawberry was regulated by MTA, the core component of the methyltransferase complexes with m⁶A catalytic activity. The MTA RNAi fruits showed lower m⁶A levels, while the MTA-overexpressed fruits displayed higher m⁶A 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 m⁶A 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 m⁶A 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 m⁶A deposition because the transcripts of these genes are not m⁶A-modified according to our m⁶A-seq datasets. We speculate that MTA may regulate translation efficiency of numerous transcripts beyond direct m⁶A installation. As expected, MTA repression or overexpression also altered the translation efficiency of a number of ripening-related genes without m⁶A 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 m⁶A-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 m⁶A 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 m⁶A 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 m⁶A modification stabilizes mRNAs of these genes. We propose that, in addition to direct modulation of translation efficiency via m⁶A installation in the target transcripts, MTA may indirectly modulate translation efficiency via m⁶A-mediated regulation of the translation initiation factors or elongation factors.

Fig. 7

Influence of m⁶A on translation initiation factors or elongation factors in strawberry. a Integrated Genome Browser (IGB) tracks displaying the distribution of m⁶A reads in transcripts of genes encoding translation initiation factors (EIF2, EIF2B, EIF3A, and EIF3C) and elongation factors (EF1A). The hypermethylated m⁶A peaks (fold change ≥ 1.5; P value < 0.05) in fruit at RS1 stage compared to those at S6 stage are indicated by shadow boxes. S6, the growth stage 6; RS1, the ripening stage 1. Rep, replicate. b m⁶A enrichment for EIF2, EIF2B, EIF3A, EIF3C, and EF1A from m⁶A-seq data. c Transcript levels of EIF2, EIF2B, EIF3A, EIF3C, and EF1A determined by RNA-seq. FPKM, fragments per kilobase of exon per million mapped fragments. d RNA immunoprecipitation (RIP) assay revealing the binding of MTA protein to the transcripts of EIF2, EIF2B, EIF3A, EIF3C, and EF1A. e–i The changes in relative m⁶A enrichment (e, g), gene expression (f, h), and mRNA stability (i), for EIF2, EIF2B, EIF3A, EIF3C, and EF1A in the MTA-silenced (RNAi-MTA) or MTA-overexpressed (OE-MTA) fruits. The relative m⁶A enrichment and gene expression were determined by m⁶A-IP-qPCR and quantitative RT-PCR, respectively. The ACTIN gene served as an internal control. For mRNA stability assay, the total RNAs were extracted after actinomycin D treatment at an indicated time point and subjected to quantitative RT-PCR assay. Data are presented as mean ± standard deviation (n = 3). Asterisks indicate significant differences (*P < 0.05, **P < 0.01, ***P < 0.001; Student’s t test)

MTA is indispensable for global m⁶A mRNA methylation

To gain a more detailed insight into the MTA-mediated m⁶A methylation, we performed m⁶A-seq on fruits of the MTA RNAi and the control. A total of 8693 and 8704 confident m⁶A 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 m⁶A 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 m⁶A peaks covering 1128 gene transcripts and 1510 hypermethylated m⁶A 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 m⁶A hypomethylation when the MTA was silenced. The hypomethylated m⁶A peaks were mainly distributed in the CDS region (76.19%), while the hypermethylated m⁶A 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 m⁶A peaks identified during fruit ripening (Fig. 3e).

Fig. 8

Changes in m⁶A mRNA methylome in MTA-silenced strawberry fruits. a Venn diagrams depicting the overlap of m⁶A peaks from three independent m⁶A-seq experiments of fruits of the MTA RNAi (RNAi-MTA) and the control. b Volcano plot showing hypermethylated (red) and hypomethylated (blue) m⁶A peaks in the MTA-silenced fruits compared to those in the control group. c Heat map of m⁶A enrichment ratios of differential m⁶A peaks shown in b. d The distribution characteristics of differential m⁶A peaks shown in b. UTR, untranslated region; CDS, coding sequence. e, f Box plot of gene expression ratios (e) and cumulative distribution of gene expression changes (f) between the MTA-silenced fruits and the controls. Expression of transcripts containing hypomethylated and non-differential m⁶A peaks in the MTA-silenced fruits compared to the controls was analyzed based on our RNA-seq data. Asterisks indicate significant differences (***P < 0.001; Wilcoxon test). g Gene ontology (GO) enrichment for the transcripts containing hypomethylated m⁶A peaks in the MTA-silenced fruits compared to the control. The GO enrichment was analyzed on the agriGO database (http://systemsbiology.cau.edu.cn/agriGOv2/), and the most impacted GO biological process terms with a statistical significance of Yekutieli-corrected P value < 0.05 were shown

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 m⁶A-modified transcripts based on our RNA-seq data, we found that genes harboring hypomethylated m⁶A modification in the MTA RNAi fruits were more likely downregulated, and this trend did not exist for non-differential m⁶A-modified genes (Fig. 8e, f), suggesting a positive correlation between MTA-mediated m⁶A methylation and mRNA abundance. GO analysis revealed that genes containing hypomethylated m⁶A 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 m⁶A mRNA methylation, which ensures the proper progression of fruit ripening of strawberry by targeting ripening-related genes among various biological processes.

Discussion

m⁶A methylation participates in the regulation of strawberry fruit ripening

Compared to the understanding of the molecular basis underlying fruit ripening in climacteric fruits such as tomato, our knowledge regarding the regulation of ripening in non-climacteric fruits, e.g., strawberry, is still limited. Recently, it was elucidated that the ripening of strawberry involves the remodeling of DNA methylation [55]. However, it is unclear whether m⁶A methylation, which has been revealed to modulate the ripening of tomato fruit [36], is involved in the regulation of strawberry fruit ripening. In the present study, we found that m⁶A methylation represents a widespread mRNA modification in strawberry fruit and exhibits a dramatic change at the initiation stage of ripening. Overexpression of MTA or MAB, the m⁶A methyltransferase genes, promotes fruit ripening, whereas repression of either gene delays ripening (Fig. 6a, b), demonstrating that m⁶A methylation participates in the regulation of strawberry fruit ripening.

The transcript levels of MTA and MTB increased significantly upon ripening initiation of strawberry fruit (Fig. 5c–e), which may account for the m⁶A hypermethylation in the CDS region. It should be noted that the m⁶A methyltransferase genes in tomato express stably and no obviously global m⁶A hypermethylation was observed during fruit ripening [36]. By contrast, it is the m⁶A demethylase SlALKBH2 that positively regulates tomato fruit ripening through mediating the mRNA stability of SlDML2, a key DNA demethylase gene determining the DNA methylation patterns during ripening [36]. Due to the pivotal role of DNA methylation in the regulation of fruit ripening in both climacteric and non-climacteric fruits, it is reasonable to assume that m⁶A modification may modulate strawberry fruit ripening by modulating the DNA methylation machinery.

Strawberry undergoes an overall loss of DNA methylation during ripening [55]. There are four SlDML2 homologs, FveDME1, FveROS1.1, FveROS1.2, and FveROS1.3, in the strawberry genome (Additional file 1: Figure S15a). Our m⁶A-seq analysis revealed that FveDME1 and FveROS1.1 contain differential m⁶A peaks upon ripening initiation (from S6 to RS1), while FveROS1.2 and FveROS1.3 are not m⁶A-modified (Additional file 1: Figure S15b). However, none of these genes showed significantly increased expression from S6 to RS1 (Additional file 1: Figure S15c), suggesting that the homologs of SlDML2 might be dispensable for m⁶A-mediated regulation of ripening in strawberry fruit. Previous study has shown that the reprogramming of DNA methylation during the ripening of strawberry fruit is governed by components in the RNA-directed DNA methylation (RdDM) pathway rather than those in the demethylation pathway [55]. Our m⁶A-seq data indicated that there was no differential m⁶A modification in the transcripts of DNA methyltransferase genes in the RdDM pathway during strawberry fruit ripening (Additional file 1: Figure S16). Together, these data suggest that distinct mechanisms underlie the m⁶A-mediated ripening regulation in strawberry fruit.

m⁶A methylation regulates strawberry fruit ripening by targeting ABA pathway

The phytohormone ABA plays crucial roles in plant growth, development, and stress responses [56]. The ripening of non-climacteric fruits, such as strawberry, has been revealed to be ABA-dependent [34]. Identification and functional definition of the core elements in ABA biosynthesis and signal transduction pathway have expanded our understanding of the mechanistic roles of ABA-mediated regulation of non-climacteric fruit ripening.

The ABA biosynthesis rate-limiting enzyme NCED was reported to play an essential role in the ripening of strawberry fruit [34, 38]. Moreover, several critical constituents of ABA signaling, including the ABA receptor FaPYR1 and FaABAR, the type 2C protein phosphatase FaABI1, and the SNF1-related kinase FaSnRK2.6, have been revealed to be indispensable for normal fruit ripening of strawberry [57,58,59]. Nevertheless, the regulatory mechanisms underlying ABA biosynthesis and signaling pathway remain largely unknown.

In Arabidopsis, ABA receptor PYLs could be regulated by post-translational modification, such as phosphorylation [60], tyrosine nitration [61], and ubiquitination [62], which synergistically modulate the abundance and activity of the receptors [56]. In contrast, the downstream signal molecules PP2Cs and SnRK2s are mainly controlled by protein phosphorylation and dephosphorylation for maintaining the appropriate ABA signal transduction [63,64,65,66]. The transcription of genes in ABA biosynthesis and signaling pathway display dynamic changes in various development processes or in response to environmental stresses [56, 67,68,69], implying that they undergo precise regulation at transcriptional level. However, little is known about the regulation of genes in the ABA pathway at the post-transcriptional level.

In this study, we found that NCED5, ABAR, and AREB1, the genes in ABA biosynthesis and signaling pathway undergo m⁶A-mediated post-transcriptional regulation (Fig. 4). The m⁶A modifications promote the mRNA stability of NCED5 and AREB1, while enhancing the translation efficiency of ABAR. These findings identify a novel layer of gene regulation in ABA biosynthesis and signaling pathway and establish a link between m⁶A-mediated ABA pathway and strawberry fruit ripening. Given the essential roles of ABA in plant development and stress resistance, it is interesting and necessary to explore the regulation of m⁶A methylation on these physiological processes.

m⁶A modification exhibits diverse effects on mRNAs in strawberry

As the most prominent modification in mRNA, m⁶A was initially thought to accelerate the degradation of mRNAs through promoting their transfer from the translatable pool to the decay sites [1]. Later findings have indicated that m⁶A modification also harbors the capacity to stabilize the mRNAs or promotes the translation in various biological processes [3, 22, 70, 71], although the underlying mechanisms are poorly understood, suggesting that m⁶A methylation possesses different regulatory roles on mRNAs [4]. This functional diversity was proposed to be heavily dominated by m⁶A distribution and the local sequence contexts within transcripts. One likely explanation is that the specific sequence contexts around m⁶A marks determine the recruitment of diverse m⁶A readers or other RNA binding proteins (RBPs) that carry distinguishing and even opposite molecular functions [4].

In this study, we observed a dramatic change in m⁶A modification at the initiation of strawberry fruit ripening (Fig. 2). The m⁶A depositions in the CDS region tend to stabilize the mRNAs, whereas those in the 3′ UTR or around the stop codon exhibit the opposite effects (Fig. 3). These data suggest that, depending on their distribution, the m⁶A modification plays distinct roles on mRNAs in strawberry. Notably, the differential m⁶A modification in the ripening process of tomato fruit mainly appeared around the stop codon or within the 3′ UTR, and these m⁶A depositions were generally negatively correlated with the abundance of the mRNAs. Combining the observation in strawberry and tomato, we propose that m⁶A modification around the stop codon or within the 3′ UTR tends to negatively correlated with mRNA abundance in the ripening process of both climacteric and non-climacteric fruits, whereas m⁶A deposition in the CDS region generally stabilizes the mRNAs upon ripening initiation of non-climacteric fruits. Interestingly, besides modulating mRNA stability, m⁶A modification also affected translation efficiency of some transcripts. The changes in translation efficiency might be directly regulated by m⁶A deposition on the transcripts, or indirectly by m⁶A-mediated regulation of translation initiation factors and elongation factors.

In conclusion, our work reveals a regulatory role of m⁶A methylation on the ripening of the non-climacteric strawberry fruit. The molecular basis, which involves the m⁶A-mediated regulation of genes in the ABA biosynthesis and signaling pathway, is distinct from what was described in climacteric tomato fruit (Additional file 1: Figure S17). These findings provide new insights into understanding the regulatory networks controlling fruit ripening.

Methods

Plant materials

The diploid woodland strawberry (Fragaria vesca cv. “Hawaii-4”) and the octoploid cultivated strawberry (Fragaria × ananassa cv. “Benihoppe”) were planted in a greenhouse with standard culture conditions [59]. To accurately determine fruit ages through development, flowers were tagged at the anthesis. Fruits of “Hawaii-4” at the growth stage 6 (S6), the ripening stage 1 (RS1), and the ripening stage 3 (RS3) [38], which were on average 15, 21, and 27 days post-anthesis (DPA), respectively, were harvested, and then frozen in liquid nitrogen. The fruits with the removal of attached achenes were subsequently stored at − 80 °C until use. The ‘Benihoppe’ fruits were harvested at the small green (SG), large green (LG), white (Wt), initial red (IR), and full red (FR) stages, respectively, based on the size, weight, shape, and color [72], and then maintained at fresh status for further studies or directly frozen and stored as the “Hawaii-4.”

m⁶A-seq and data analysis

The m⁶A-seq was performed according to the method described by Dominissini et al. (2013) [37]. Briefly, total RNAs were extracted from the woodland strawberry fruit at S6, RS1, and RS3 stages, and the RNAi-MTA and control fruit by the plant RNA extraction kit (Magen, R4165-02), and then 300 μg of intact total RNAs were used for mRNA isolation by the Dynabeads mRNA purification kit (Life Technologies, 61006). The purified mRNAs were randomly fragmented into ~ 100 nucleotide-long fragments by incubation at 94 °C for 5 min in the RNA fragmentation buffer (10 mM Tris-HCl, pH 7.0, and 10 mM ZnCl₂). The reaction was terminated with 50 mM EDTA, and then the fragmented mRNAs were purified by standard phenol-chloroform extraction and ethanol precipitation. For immunoprecipitation (IP), 5 μg of fragmented mRNAs was mixed with 10 μg of anti-m⁶A polyclonal antibody (Synaptic Systems, 202003) and incubated at 4 °C for 2 h in 450 μL of IP buffer consisting of 10 mM Tris-HCl, pH 7.4, 150 mM NaCl, 0.1% NP-40 (v/v), and 300 U mL^-1 RNase inhibitor (Promega, N2112S). After the addition of 50 μL Dynabeads protein-A (Life Technologies, 10002A), the mixture was incubated at 4 °C for another 2 h. The beads were subsequently washed twice with high-salt buffer containing 50 mM Tris-HCl, pH 7.4, 1 M NaCl, 1 mM EDTA, 1% NP-40 (v/v), and 0.1% SDS (w/v) and twice with IP buffer. The m⁶A-containing fragments were eluted from the beads by incubation with 6.7 mM N⁶-methyladenosine (TargetMol, T6599) in IP buffer at 4 °C for 2 h, followed by standard phenol-chloroform extraction and ethanol precipitation. Then, 50 ng of m⁶A-containing mRNAs or pre-immunoprecipitated mRNAs (the input) were used for library construction by the NEBNext ultra RNA library preparation kit (NEB, E7530). High-throughput sequencing was conducted on the Illumina HiSeq X sequencer with a paired-end read length of 150 bp following the standard protocols. The sequencing was performed with three independent biological replicates, and each RNA sample was prepared from the mix of at least 60 strawberry fruits to avoid individual differences among fruits.

For data analysis, the quality of raw sequencing reads was initially assessed by FastQC tool (version 0.11.7; http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Adaptors and low-quality bases with a score < 20 were trimmed using Cutadapt (version 1.16) [73], and then short reads with a length < 18 nucleotides or reads containing ambiguous nucleotides were filtered out by Trimmomatic (version 0.30) [74]. The remaining reads were aligned to the strawberry reference genome v2.0.a1 (ftp://ftp.bioinfo.wsu.edu/species/Fragaria_vesca/Fvesca-genome.v2.0.a1) by Burrows-Wheeler Aligner (BWA; version 0.30) [75]. Mapping quality (MAPQ) of all aligned reads was concurrently evaluated, and only uniquely mapped reads with a MAPQ ≥ 13 were retained for subsequent analysis [37].

The identification of m⁶A peaks was carried out by MACS software (version 2.0.10) [76], using the corresponding input as a control. High-confidence peaks were obtained by a stringent cutoff threshold for MACS-assigned false discovery rate (FDR) < 0.05. PeakAnalyzer (version 2.0) [77] was applied to annotate the identified peaks to the strawberry genome annotation file (ftp.bioinfo.wsu.edu/species/Fragaria_vesca/Fvesca-genome.v2.0.a2/genes/). To identify the differentially methylated peaks between samples, the m⁶A site differential algorithm [78] was applied with a criterion of fold change in m⁶A enrichment ≥ 1.5 and P value < 0.05. HOMER (version 4.7; http://homer.ucsd.edu/homer/) [42] was employed to identify the m⁶A motifs with a restricted length of six nucleotides. The differential m⁶A peaks identified between S6 and RS1 stages were used as the target sequences, and the exon sequences without m⁶A peaks were used as the background sequences. Integrated Genome Browser (IGB, version 9.0.2) [79] was used for the visualization of m⁶A peaks. Gene Ontology (GO) enrichment was analyzed on the agriGO database (version 2.0; http://systemsbiology.cau.edu.cn/agriGOv2/) and only statistically significant terms with a Yekutieli-corrected P value < 0.05 were remained [80].

RNA-seq and data analysis

The sequencing reads from the input samples in m⁶A-seq were used for RNA-seq analysis as previously reported [21]. In brief, the uniquely mapped reads with a MAPQ ≥ 13 were assembled by Cufflinks [81]. Gene expression was presented as fragments per kilobase of exon per million mapped fragments (FPKM) by using Cuffdiff [81], which concurrently provides statistical routines for capturing differentially expressed genes. The Benjamini and Hochberg’s approach [82] was used to adjust the resulting P values for controlling the false discovery rate (FDR). Differential gene expression was defined on basis of a cutoff criterion of FPKM fold change ≥ 1.5 and P value < 0.05.

Quantitative RT-PCR analysis

Total RNAs were extracted from the strawberry fruits or N. benthamiana leaves using the plant RNA extraction kit (Magen, R4165-02). The extracted RNAs were reverse transcribed into cDNAs by the HiScript® III RT SuperMix for qPCR kit (Vazyme, R323-01). The synthesized cDNAs were then employed as templates for PCR amplification using the ChemQ Universal SYBR qPCR Master Mix (Vazyme, Q711-02-AA) and a StepOne Plus Real-Time PCR System (Applied Biosystems) in the following program: 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 30 s. Relative quantification of gene transcription levels were performed by the cycle threshold (C_T) 2^(−ΔCT) method [83]. Strawberry ACTIN (gene22626) or N. benthamiana ACTIN (Niben101Scf03410g03002) was applied to normalize the expression values. The primers for PCR amplification are listed in Additional file 18: Table S17. The experiment was conducted with three biological replicates, and each contained three technical repeats.

m⁶A-IP-qPCR

m⁶A-IP-qPCR was carried out as previously described with minor modifications [84]. Briefly, 5 μg of intact mRNAs were fragmented into ~ 300 nucleotide-long fragments by an incubation at 94 °C for 30 s in the RNA fragmentation buffer (10 mM Tris-HCl, pH 7.0, and 10 mM ZnCl₂), followed by the addition of 50 mM EDTA to terminate the reaction. The fragmented mRNAs were purified by standard ethanol precipitation and resuspended in 250 μL of DEPC-treated water. Then, 5 μL of fragmented mRNAs was taken as the input control and 100 μL were incubated with 5 μg of anti-m⁶A polyclonal antibody (Synaptic Systems, 202003) at 4 °C for 2 h in 450 μL of IP buffer containing 10 mM Tris-HCl, pH 7.4, 150 mM NaCl, 0.1% NP-40 (v/v), and 300 U mL^-1 RNase inhibitor (Promega, N2112S). The mixture was subsequently incubated with 20 μL of Dynabeads protein-A (Life Technologies, 10002A) at 4 °C for another 2 h. After washing twice with high-salt buffer containing 50 mM Tris-HCl, pH 7.4, 1 M NaCl, 1 mM EDTA, 1% NP-40 (v/v), and 0.1% SDS (w/v) and twice with IP buffer, the m⁶A-containing fragments were eluted by an incubation with 6.7 mM N⁶-methyladenosine (m⁶A; TargetMol, T6599) in 200 μL of IP buffer at 4 °C for 2 h, followed by ethanol precipitation. The immunoprecipitated mRNA fragments were finally resuspended in 5 μL DEPC-treated water. Then, both the m⁶A-containing mRNAs and the input mRNAs were submitted to quantitative RT-PCR by using the primers listed in Additional file 18: Table S17. Relative m⁶A enrichment in a specific region of a transcript was calculated using the cycle threshold (C_T) 2^(−ΔCT) method [83]. The value for the immunoprecipitated sample was normalized against that for the input. The experiment was performed with three biological replicates, and each contained three technical repeats.

SELECT analysis

To validate the specific m⁶A sites in mRNAs, the SELECT analysis was performed according to the reported method [49]. Briefly, 20 ng of mRNAs were mixed with 40 nM Up-Primer, 40 nM Down-Primer, and 5 μM dNTP in 17 μL 1× CutSmart buffer containing 20 mM Tris-HAc, pH 7.9, 50 mM KAc, 10 mM MgAc2, and 100 μg ml^-1 BSA. The mixture was annealed by incubating at a temperature gradient: 90 °C for 1 min, 80 °C for 1 min, 70 °C for 1 min, 60 °C for 1 min, 50 °C for 1 min, and 40 °C for 6 min. Then, 0.01 U Bst 2.0 DNA polymerase (NEB, M0537S), 0.5 U SplintR ligase (NEB, M0375S), and 10 nmol ATP were added to the mixture to a final volume of 20 μL. The final mixture was incubated at 40 °C for 20 min, denatured at 80 °C for 20 min, and kept at 4 °C. Subsequently, the mixture was diluted 1:50 in DEPC-treated water and submitted to quantitative RT-PCR as the template. Data was analyzed by the StepOne Plus Real-Time PCR Software (version 2.3). All used primers are listed in Additional file 18: Table S17.

Identification of strawberry m⁶A methyltransferases and phylogenetic analysis

To identify m⁶A methyltransferases in strawberry, the amino acid sequence of conserved MT-A70 domain (PF05063) was downloaded from the Pfam (http://pfam.xfam.org/) and then utilized to search the potential homologs against the strawberry protein dataset (ftp.bioinfo.wsu.edu/species/Fragaria_vesca/Fvesca-genome.v2.0.a2/genes/) using HMMER 3.1 with default parameters [85]. The resulting protein sequences were analyzed on the CDD database (https://www.ncbi.nlm.nih.gov/cdd/) [86], and only that containing a MT-A70 domain remained as the m⁶A methyltransferase candidates. The known m⁶A methyltransferases METTL3 and METTL14 in mammals and MTA and MTB in Arabidopsis were then employed to perform BLASTP-algorithm in NCBI (https://www.ncbi.nlm.nih.gov/) to further confirm the identified homologs. For phylogenetic analysis, the sequences of the identified strawberry m⁶A methyltransferase were aligned with the sequences of m⁶A methyltransferase in mouse (Mus musculus), rice (Oryza sativa), maize (Zea mays), tomato (Solanum lycopersicum), tobacco (Nicotiana benthamiana), and Arabidopsis, using ClustalX 2.1 with standard parameters [87]. The alignment result was imported into MEGA software (version 5.2) to create the phylogenetic tree by using Neighbor-Joining method with 1000 bootstrap replicates [88].

Y2H analysis

Y2H analysis was performed as the previously described [89]. In brief, the coding sequence of MTA and MTB was amplified from the cDNAs of diploid woodland strawberry and then ligated into the pGADT7 (AD) and pGBKT7 (BD) vectors, respectively. The resulting plasmids were co-transformed into Saccharomyces cerevisiae strain AH109 according to the protocols in the Matchmaker GAL4 Two-Hybrid System 3 (Clontech). The yeast cells were cultured on SD/-Leu-Trp (-LW) medium, and then transferred onto the SD/-Leu-Trp-His (-LWH) or SD/-Leu-Trp-His-Ade medium (-LWHA) with or without X-α-gal. Transformants carrying empty pGADT7 (AD) or pGBKT7 (BD) vectors were used as negative controls, and that concurrently carrying pGBKT7-53 and pGADT7-T vectors was used as the positive control. The primers used for vector constructions are listed in Additional file 18: Table S17.

LCI assay

LCI assay was performed as previously reported [90]. Briefly, the coding sequence of MTA was cloned into the pCambia1300-nLUC plasmid, and the coding sequence of MTB or its MT-A70 domain sequence (MTB^D) was cloned into the pCambia1300-cLUC plasmid. The resulting constructs were separately transformed into Agrobacterium tumefaciens strain GV3101. The agrobacteria were cultured at 28 °C for 18 h in Luria-Bertani (LB) liquid medium containing 50 μg mL^− 1 kanamycin, 50 μg mL^− 1 gentamycin, and 50 μg mL^− 1 rifampicin. After being pelleted by centrifugation at 5,000 g for 5 min, the agrobacteria were resuspended in the infiltration buffer (10 mM MES, pH 5.6, 10 mM MgCl₂, and 100 μM acetosyringone) to a final OD₆₀₀ of 0.5. Then, the suspension was infiltrated into N. benthamiana leaves for co-expression of fusion protein nLUC-MTA with MTB-cLUC or MTB^D-cLUC. After culture for 36 h, the leaves were incubated with 1 mM luciferin dissolved in ddH₂O supplemented with 0.01% Triton X-100 at room temperature for 5 min, and then observed under a chemiluminescence imaging system (Tanon). Empty vectors expressing cLUC or nLUC were co-transformed as the negative controls. The primers used for vector constructions are listed in Additional file 18: Table S17.

Subcellular localization

For subcellular localization analysis, the coding sequence of MTA and MTB was amplified from the cDNAs of diploid woodland strawberry and then inserted into the pCambia2300-mCherry and pCambia2300-eGFP plasmids to generate 35S::MTA-mCherry and 35S::MTB-eGFP vectors, respectively. The resulting constructs were separately transformed into A. tumefaciens strain GV3101. The agrobacteria were subsequently infiltrated into N. benthamiana leaves for the individual expression of mCherry-tagged MTA (MTA-mCherry) and eGFP-tagged MTB (MTB-eGFP) or the co-expression of the two fusion proteins. After culture for 36 h, the mesophyll protoplasts were isolated from N. benthamiana leaves as previously reported [91] and observed under a Leica confocal microscope (Leica DMI600CS). Protoplasts expressing eGFP or mCherry were used as negative controls. The primers used for vector constructions are listed in Additional file 18: Table S17.

Agroinfiltration-mediated transient transformation in strawberry fruit

Transient transformation of strawberry fruit mediated by agroinfiltration was performed as previously described [54]. To construct the RNA interference (RNAi) vectors, a ~ 300-bp fragment targeting the coding sequence region of MTA or MTB was cloned and inserted into the pCR8 plasmid, and then restructured into the pK7GWIWGD (II) plasmid by using the Gateway LR Clonase^TM Enzyme Mix (Invitrogen, 11791-020). To construct the overexpression (OE) vectors, the coding sequence of MTA and MTB was amplified and ligated into the pCambia2300-eGFP plasmid to generate 35S::MTA-eGFP and 35S::MTB-eGFP vector, respectively. The resulting constructs were separately transformed into the A. tumefaciens strain GV3101. The agrobacteria were cultured at 28 °C overnight in LB liquid medium supplemented with 50 μg mL⁻¹ kanamycin, 50 μg mL⁻¹ gentamycin, and 50 μg mL⁻¹ rifampicin, and then diluted 1:100 in 100 mL of fresh LB medium to continue culturing for approximately 8 h. The agrobacteria cells were subsequently collected by centrifugation at 5,000 g for 5 min and resuspended in the infiltration buffer (10 mM MES, pH 5.6, 10 mM MgCl₂, and 100 μM acetosyringone) to a final OD₆₀₀ of 0.8. After being kept at room temperature for 2 h without shaking, the suspensions were injected into the octoploid strawberry fruit at large green (LG) stage by using a 1 mL syringe. The infiltrated fruits were cultured for 5–7 days in a growth room with the following conditions: 23 °C, 80 % relative humidity, and a 16/8-h light/dark photoperiod with a light intensity of 100 μmol m⁻² s⁻¹. The experiment was performed with more than three independent biological replicates, and each group contained at least fifteen fruits. The primers used for vector constructions are listed in Additional file 18: Table S17.

mRNA stability assay

For mRNA stability assay in N. benthamiana leaves, the coding sequence of NCED5, ABAR, and AREB1 was amplified from the cDNAs of diploid woodland strawberry. The mutated forms of the amplified sequence in which the potential m⁶A sites were replaced by cytidine (C) or guanine (G) were constructed using the QuikChange II XL Site-directed Mutagenesis Kit (Agilent Technologies, 200518). The fragments were separately inserted into the pCambia2300 vector, which were subsequently transformed into A. tumefaciens strain GV3101. Then, the agrobacteria were cultured at 28 °C for 18 h in LB liquid medium containing 50 μg mL⁻¹ kanamycin, 50 μg mL⁻¹ gentamycin, and 50 μg mL⁻¹ rifampicin. After being collected by centrifugation at 5,000 g for 5 min, the agrobacteria were diluted to an OD₆₀₀ of 0.5 in the infiltration buffer (10 mM MES, pH 5.6, 10 mM MgCl₂, 100 μM acetosyringone), and then infiltrated into the N. benthamiana leaves. After 36 h of incubation, the infiltration parts in the leaves were injected with 20 μg mL^− 1 actinomycin D (Sigma, A4262) dissolved in ddH₂O. After culture for 30 min, leaf discs were taken and considered as time 0 controls, and subsequent samples were harvested every 3 h in triplicate.

For mRNA stability assay in strawberry, the MTA RNAi fruit, the MTA-overexpressed fruit, and the controls were sliced into ~ 3 mm slices and then transferred onto plates containing 20 μg mL^− 1 actinomycin D dissolved in ddH₂O. After incubation for 30 min, the slices were collected as the N. benthamiana leaf discs. The mRNA levels of genes were subsequently examined by quantitative RT-PCR as described above. All primers used for PCR amplification were listed in Additional file 18: Table S17.

Translation efficiency assay

Translation efficiency was assayed according to the method described by Merchante et al. (2015) [50]. Briefly, 6 g of strawberry fruits or 3 g of N. benthamiana leaves were ground into fine powder in liquid nitrogen. One gram of sample was used for total RNA extraction, and the rest was suspended in 15 mL of polysome extraction buffer (200 mM Tris-HCl, pH 9.0, 35 mM MgCl₂, 200 mM KCl, 25 mM EGTA, 1% Triton X-100 (v/v), 1% IGEPAL CA-630 (v/v), 5 mM DTT, 1 mM PMSF, 50 μg mL⁻¹ Chloramphenicol, and 100 μg mL⁻¹ Cycloheximide) at 4 °C for 20 min with slight shaking. The mixture was centrifuged twice at 16,000g for 20 min at 4 °C. Then, 12.5 mL of supernatant were slowly transferred onto 13.5 mL of sucrose buffer (1.75 M Sucrose, 400 mM Tris-HCl, pH 9.0, 35 mM MgCl₂, 5 mM EGTA, 200 mM KCl, 5 mM DTT, 50 μg mL⁻¹ Chloramphenicol, and 50 μg ml⁻¹ Cycloheximide). After centrifugation at 200,000g for 4 h at 4 °C, the supernatant was carefully removed, and the polysomes in the bottom were resuspended in 300 μL of DEPC-treated water. The polysomal RNAs, as well as the total RNAs, were isolated by using the plant RNA extraction kit (Magen, R4165-02), and used for quantitative RT-PCR analysis as described above. Translation efficiency was calculated by the abundance ratio of mRNA in the polysomal RNA versus the total RNA using the cycle threshold (CT) 2^(−ΔCT) method [83] with the GADPH2 gene as an internal reference. The primers used for PCR amplification are listed in Additional file 18: Table S17.

Quantitative analysis of m⁶A level by LC-MS/MS

The global m⁶A levels in strawberry fruit were detected by LC-MS/MS as previously described [36] with minor modifications. In brief, 200 ng of mRNAs was digested at 37 °C for 6 h with 1 unit of Nuclease P1 (Wako, 145-08221) in 50 μL of reaction buffer containing 10 mM ammonium acetate, pH 5.3, 25 mM NaCl, and 2.5 mM ZnCl₂. Then, 1 unit of alkaline phosphatase (Sigma-Aldrich, P6774) and 5.5 μL of 1 M fresh NH₄HCO₃ were added, followed by incubation at 37 °C for another 6 h. After centrifugation at 15,000 g for 5 min, the supernatant was used for LC-MS/MS analysis. The digested nucleosides were separated by UPLC (Waters, ACQUITY) equipped with a ACQUITY UPLC HSS T3 column (Waters), and then detected by a Triple Quad Xevo TQ-S (Waters) mass spectrometer in positive ion mode with multiple reaction monitoring. The mobile phase was composed of buffer A (0.1% formic acid in ultrapure water) and buffer B (100% acetonitrile). Nucleosides were accurately quantified depending on the nucleoside-to-base ion mass transitions of m/z 268.0 to 136.0 (A) and m/z 282.0 to 150.1 (m⁶A). The pure commercial adenosine (A; TargetMol, T0853) and N⁶-methyladenosine (m⁶A; TargetMol, T6599) were used to generate standard curves, which were subsequently employed to calculate the contents of A and m⁶A in each sample. The global m⁶A levels were presented in the form of m⁶A/A ratio. The experiment was repeated with three independent biological replicates.

RNA immunoprecipitation

RNA immunoprecipitation (RIP) was carried out following the method of Wei et al. (2018) [22] with minor modifications. Briefly, the octoploid strawberry fruit expressing the MTA-eGFP protein were sliced into ~ 2 mm slices, and then fixed with 1% formaldehyde under a vacuum for 30 min on ice. The fixation was terminated by the addition of 150 mM glycine, followed by incubation on ice for 5 min. The fixed fruit tissues (2 g) were homogenized in 5 mL of lysis buffer (50 mM HEPES, pH 7.5, 2 mM EDTA, 150 mM KCl, 0.5% NP-40 (v/v), 2 mM EDTA, 0.5 mM DTT, 1× cocktail protease inhibitor (Sigma, 04693132001), and 300 U mL⁻¹ RNase Inhibitor). After incubation at 4 °C for 1 h, the mixture was centrifuged at 15,000g for 30 min. Then, 200 μL of the supernatant were taken as the input control, and the remainder was subjected to immunoprecipitation (IP) with anti-GFP monoclonal antibody or rabbit lgG at 4 °C overnight. Fifty microliters of Dynabeads Protein-A (Life Technologies, 10002A) were added to the mixture and incubated at 4 °C for 2 h. After washing four times with PBS buffer, the RNA-protein mix was catalyzed by proteinase K (Takara, 9034) at 55 °C for 1 h. The immunoprecipitated mRNAs and input mRNAs were subsequently isolated by the plant RNA extraction kit (Magen, R4165-02). Relative enrichment of individual transcript was determined by quantitative RT-PCR analysis as described above. The primers used for PCR amplification are listed in Additional file 18: Table S17. The analysis was performed with three biological replicates, and each contained three technical repeats.