Skip to main content

Transcriptome-wide high-throughput deep m6A-seq reveals unique differential m6A methylation patterns between three organs in Arabidopsis thaliana

Abstract

Background

m6A is a ubiquitous RNA modification in eukaryotes. Transcriptome-wide m6A patterns in Arabidopsis have been assayed recently. However, differential m6A patterns between organs have not been well characterized.

Results

Over two-third of the transcripts in Arabidopsis are modified by m6A. In contrast to a recent observation of m6A enrichment in 5′ mRNA, we find that m6A is distributed predominantly near stop codons. Interestingly, 85 % of the modified transcripts show high m6A methylation extent compared to their transcript level. The 290 highly methylated transcripts are mainly associated with transporters, stress responses, redox, regulation factors, and some non-coding RNAs. On average, the proportion of transcripts showing differential methylation between two plant organs is higher than that showing differential transcript levels. The transcripts with extensively higher m6A methylation in an organ are associated with the unique biological processes of this organ, suggesting that m6A may be another important contributor to organ differentiation in Arabidopsis. Highly expressed genes are relatively less methylated and vice versa, and different RNAs have distinct m6A patterns, which hint at mRNA fate. Intriguingly, most of the transposable element transcripts maintained a fragmented form with a relatively low transcript level and high m6A methylation in the cells.

Conclusions

This is the first study to comprehensively analyze m6A patterns in a variety of RNAs, the relationship between transcript level and m6A methylation extent, and differential m6A patterns across organs in Arabidopsis.

Background

Over 100 types of chemical modifications have been discovered in RNAs from all of the living species [1, 2]. The most diverse modifications were present in ribosomal RNA (rRNA) and transfer RNA (tRNA). Chemical modifications are also prevalent in messenger RNA (mRNA) and other non-coding RNA (ncRNA) in eukaryotes [1, 2]. Among those, the most important is modified by N 6-methyladenosine (m6A) [2–4]. m6A has been found ubiquitously distributed in rRNA, tRNA, mRNA, and some snRNA of eukaryotes, such as yeast [5], mammals [4, 6], insects [2], and plants [7]. Recently, transcriptome-wide analyses showed that one-third of the transcribed genes (mRNA) were modified by m6A in human and mouse [4, 6, 8]. The m6A enriched sites were found near stop codons, in 3′UTRs and mRNA segments derived from large exons [4, 6, 8]. These studies also showed that this modification was highly conserved in eukaryotes [2, 6], suggesting that a delicate regulatory mechanism may be responsible for this selective modification, and provided clues of the important metabolisms that this modification involved in or was responsible for, for example, RNA splicing [6, 8], RNA export [4], and RNA stability [4, 6].

The availability of antibody that specifically binds the m6A modified sites and efficiently enriches RNAs containing m6A modification facilitates the transcriptome-wide analysis of the patterns of this RNA modification through the biotechnologies of RNA sequencing (RNA-seq), RNA immunoprecipitation (RIP), and m6A-seq [4, 6, 8]. RIP was primarily used to analyze RNA-protein interaction [9]. However, the aim of the RIP experiment for the m6A-seq study was to pull down the RNA of interest containing m6A modification through application of m6A antibody to the randomly fragmented RNA pool. m6A-seq is a recently reported technology integrating the powers of both RIP and high-throughput RNA sequencing for transcriptome-wide analysis of m6A methylation patterns in eukaryotes [6, 10].

Transcriptome-wide analysis of m6A in mammals and plants provided insights into topological patterns and facilitated discovery of some functions of this RNA modification [4, 6, 10–13]. However, the differential m6A methylation among plant organs, for example, leaves, flowers, and roots, has not been well characterized. In this study, we significantly improved biotechnologies for RNA isolation and RIP, thus deep and high quality m6A-seq and massively m6A-mapped datasets in Arabidopsis are now available. This study aimed to: (1) comprehensively and transcriptome-wide characterize the m6A distributing patterns in numerous types of RNAs in Arabidopsis; (2) analyze the relationship between the transcript level and the m6A modification extent in the Arabidopsis transcriptome; (3) characterize differential patterns of the m6A methylation among leaves, flowers, and roots; and (4) discuss new functions of m6A modification in the transcripts extensively modified by m6A from the clues of the potential biological functions in these transcripts. This is the first study to comprehensively analyze m6A differential patterns across organs in plants. This study opens up a new avenue to greatly understand the transcriptome-wide patterns of m6A modification in different RNAs, relationship between m6A methylation extent and gene transcript level, and m6A differential patterns across organs in plants.

Results

Quality and depth of the RNA sequencing in this study

Commercial m6A antibody has proved to specifically bind to m6A RNA and has been successfully used for the m6A RNA immunoprecipitation experiments in the previous studies [4, 6, 10–13]. In this study, we collected samples from three organs of Arabidopsis: leaves, flowers, and roots; and performed m6A-seq, mRNA-seq, and input RNA-seq (total fragmented RNA without RIP experiments for sequencing and as the control for m6A-seq) with two replicates for each sample (Additional file 1). A total of 90 to 156 million reads were generated for each m6A-seq sample; 48 to 92 million reads for each mRNA-seq sample and 25 to 53 million reads for each input RNA-seq sample (Additional file 1). The proportion of the cleanly mapped reads and transcripts in m6A-seq were around 65–70 % (Additional file 1). Compared to m6A-seq data in the mammalian (11 to 24 million reads for each sample) [6] and in the rice (23 to 47 million reads for each sample) [12], the depth of the m6A-seq in this study (Additional file 1) was greatly high. The HPLC-MS/MS results indicated that RIP efficiency for m6A enrichment was high in this study (Additional file 2), and the non-specific immunoprecipitation rate was extensively low (lower than 1 %) in this study (Additional file 2), suggesting that the experimental error caused by non-precipitation was also low in the m6A-seq in this study.

General features and extent of m6A methylation in Arabidopsis

We identified that 16,688 to 19,305 transcripts were modified by m6A in the three Arabidopsis organs. For all three organs, at least 83 % agreement was found between two m6A-seq replicates in this study (Fig. 1). This agreement proportion between replicates was the highest compared to the previous reports [6, 13]. We found that 70.6 %, 73.7 %, and 76.7 % of the transcribed genes (transcripts) were chemically modified by m6A in the leaves, flowers and roots of Arabidopsis, respectively (Additional file 3). This estimation was greatly higher than previous reports (over one-third) in human, mouse [6] and (approximately 50 %) plant [13]. The estimation differences may be due to the different criteria used for calling of the ‘m6A modified transcripts’ in m6A-seq [6, 13]. In the previous studies, input data were used as the control for calling the m6A modified transcripts due to relatively high non-specific immuno-precipitation rate in their experiments. This could result in underestimation of both total m6A peaks and proportion of the m6A transcripts in the transcriptome [6, 13]. However, all the mapped reads after removal of PCR duplicates in the m6A-seq were counted and considered to be derived from RNA fragments containing m6A modification in this study. And thus the transcripts with m6A mapped reads were considered as the modified transcripts due to very low non-specific immuno-precipitation rate (lower than 1 %) in this study (Additional file 2). But our estimation may better reflect the m6A modification extent in the cells (see details in the Discussion section of this paper).

Fig. 1
figure 1

Number of overlapped m6A transcripts in the two m6A-seq replicates of three organs of leaves (a), flowers (b), and roots (c) in Arabidposis. The analysis indicated that over 83 % of m6A transcripts were overlapped in the two m6A-seq replicates

On average of two replicates, over 80 % of the m6A modified transcripts were common among the leaves, flowers, and roots (Fig. 2). On average, around 32,300 m6A sites from the leaves, approximately 43,400 sites from the flowers, and approximately 48,100 from the roots were successfully mapped to the Arabidopsis genome with an estimation of approximately 2.0 to 2.6 m6A sites per m6A transcript and approximately 1.4 to 2.0 m6A sites per transcript in the whole transcriptome (Additional file 4). This estimation was comparable to that in mammals (approximately 1.5 m6A sites per transcript) [6], but higher than the estimated in the plant (0.7 to 1.0 m6A site per transcript) [13]. However our observation was closer to the earlier reports, for example, approximately 3-5 m6A sites per transcript [1, 2]. The ratio of m6A/A in the m6A modified transcripts was in the range of 0.44 % to 0.61 % in the three organs, and this ratio was in a range of 0.35 % to 0.50 % in the whole transcriptome of three organs in Arabidopsis (Additional file 5). This ratio estimation is also comparable to the recently reported in the plant [13].

Fig. 2
figure 2

Number of the overlapped transcripts in the two (a, b) m6A-seq replicates

m6A topological patterns in Arabidopsis

Over two-third of the methylated transcripts exhibited one or two m6A sites (Fig. 3a, details in Additional file 6). And over 17 % contained four or more sites (Fig. 3a, details in Additional file 6), which was much higher than previously reported (only 5.5 %) in human [6].

Fig. 3
figure 3

The general m6A patterns in Arabidopsis. a Proportion of the modified transcripts containing different m6A sites. b Sequence logo representing the most common consensus motif (RRm6ACH) in the m6A peaks in Arabidopsis

The consensus sequence of m6A modification has been identified as ‘RRm6ACH’, where R is A/G and H is A/C/U [4, 6, 14]. Our data showed that over 75 % of the RIP fragments in the m6A-seq contained the consensus sequence RRm6ACH in Arabidopsis (Fig. 3b). The most two frequent motifs were AAm6ACU (19.3 %) and AAm6ACA (19.0 %) (Fig. 3b). This observation is consistent with the recent reports in plants [12, 13].

The m6A topology in mRNA was categorized into two types according to the m6A patterns distributing in the whole transcript. One type was characterized by dominant m6A enrichment observed at stop codon or 3′UTR. And 73.0 % to 76.3 % of the mRNA in three organs was modified by this type of the m6A pattern (Fig. 4a, Additional file 7). Thus, most of the methylated mRNA was characterized by this typical m6A topology in Arabidopsis (Type 1, Fig. 4a): one or two high peaks at stop codon or at 3′UTR with extremely low m6A signals observed in the coding regions. In most cases, the peak height in 3′UTR or at stop codon was two to a dozen of folds of the signal heights in the coding regions (Fig. 4a, Fig. 5). This dominant m6A enrichment was not found in the remaining mRNA (Type 2, Fig. 4b). The overall m6A signals were also relatively low in Type 2 (Fig. 4b). Transcriptome-wide analysis showed that the overall m6A patterns distributing within genes were highly close with each other among three organs (Fig. 5). Statistic analysis indicated that the normalized read depth representing the overall m6A patterns had non-significant differences among three organs (P = 0.716, Additional file 8), suggesting that recognition of motif for m6A methylation was extensively conserved among plant organs.

Fig. 4
figure 4

Schematic screen shots of two typical types of m6A topologies in mRNA in Arabidopsis. a Type 1 (representative gene, ‘AT2G01910’; trace files of two organs (flowers and roots) were presented), one or two predominant peaks in 3′UTR or at stop codon with several much lower signals in the codon regions. The peaks in 3′UTR or at stop codon were two to tens of folds of the signals in the codon regions. Most of the messenger RNA (over 70 %) presented this type of m6A topology. b Type 2 (representative gene, ‘At3g29030’; trace files of three organs, leaves (the upper), flowers (in the middle), and roots (the lower) were presented). Several m6A sites distributing in the transcripts with low m6A signals (in the middle and lower parts of the figure, representative of flowers and roots) or without peaks in 3′UTR or at stop codon (in the upper part, representative of leaves)

Fig. 5
figure 5

The overall m6A distributing pattern from 5′ (left) to 3′ (right) in the m6A-seq datasets and the overall transcriptional pattern from 5′ (left) to 3′ (right) in the mRNA-seq datasets in the three organs of Arabidopsis. The patterns with the caption of ‘CK’ in the image was deduced from the mRNA-seq datasets. The other threes were from the m6A-seq dataset. The number ‘0’ on the left refers to TSS. The numbers from 49 to 59 represents stop codon or the proximate 3′UTR. As shown in this figure, only one dominant peak of m6A enrichment was detected around 3′UTR or stop codon in the overall dataset of the Arabidopsis transcriptome in this study

Two types of m6A patterns were observed in rRNAs: one was modified by one m6A site, and the other was methylated by several m6A sites (Fig. 6a and b). The m6A topology in tRNAs was also categorized into two major types: approximately 10 % of tRNAs were slightly modified by m6A (Fig. 6c) and m6A methylation was not observed in the remaining tRNAs. Both snRNA and snoRNA were highly methylated by m6A, but only a single m6A site was found in these two types of RNAs (Fig. 6d). Therefore, different types of RNAs were endowed with distinct m6A topologies.

Fig. 6
figure 6

Schematic screen shots of the m6A patterns in rRNA, tRNA and sn(o)RNA. a Numerous m6A sites distributed in a rRNA transcript (reprehensive ‘AT2g01010’). b A single m6A in a rRNA transcript (representative ‘AT2g01020’). c Low m6A methylation in most of tRNA (approximately 90 %, reprehensive ‘AT1g06480’). d sn(o) RNA was highly methylated by a single m6A site (representative ‘AT3G56705’)

m6A methylation extent versus gene transcript level

To compare m6A methylation extent in m6A-seq with gene transcript level (TL) in mRNA-seq in three organs of Arabidopsis, we categorized the m6A methylation extent into three groupings based on the comparison of modified Fragments Per Kilobase of Transcript Per Million Fragments Mapped (MFPKM) of the transcript in the m6A-seq with the FPKM of the counterpart in the mRNA-seq using χ 2 test (to test whether ratio of MFPKM to FPKM in each transcript fits 1:1 in each organ, P <0.05) (Table 1). The ratio of MFPKM to FPKM is higher (or lower) than 1 by χ 2 test (P <0.05), representing that the m6A methylation extent in m6A-seq was relatively high (or low) compared to the transcript level in mRNA-seq. When the ratio of MFPKM to FPKM in each transcript fits 1:1 (P <0.05), the m6A methylation extent was considered relatively ‘equivalent’ to the transcript level. This comparison was applied to estimate the m6A methylation extent of a transcript compared to its transcript level. All the m6A methylated transcripts were used for this estimation and the transcripts without m6A methylation detected were excluded for this comparison.

Table 1 Three groupings of the m6A methlylation extent compared to the transcript level in three organs of the Arabidopsis m6A transcriptome

Interestingly, 83.4 % to 89.1 % of the methylated transcripts showed a high m6A modification level, while <4 % of the m6A modified transcripts had a low modification extent versus their transcript level in mRNA-seq in the three organs (Table 1). Thus, most of the m6A modified transcripts showed a relatively high m6A modification extent versus their gene transcript level in the three organs of Arabidopsis.

To further analyze relationship between the m6A methylation extent in m6A-seq and the transcript level in mRNA-seq in three organs of Arabidopsis, the transcript level was categorized into three groupings: high, moderate and low. Each category contained one-third of the m6A modified transcripts from the highest to the lowest FPKM in mRNA-seq. The comparison of ratio of the average MFPKM in m6A-seq to the average FPKM in mRNA-seq between three groupings using t-test (Table 2) showed that most of the highly expressed transcripts were relatively less modified by m6A, and most transcripts with a low expression level were more likely modified by m6A (P <0.005). The moderately expressed transcripts tended to be moderately methylated in Arabidopsis (P <0.005). This observation somewhat differed from the report in human and mouse: the moderately expressed transcripts were more likely to be methylated, and those expressed at the two extremes were less methylated by m6A [6]. Both of our two groups used the similar methods to analyze the relationship between the m6A methylation extent and the gene expression level [6]. The reasons causing the differences may not be due to different methodologies or different biological species because a number of studies have proved that m6A methylation is highly conserved among eukaryotes [2, 6, 8, 12]. However, our observation may authentically reflect m6A methylation in the cells: the genes with a lower transcript level may require a relatively higher m6A modification extent to maintain RNA stability in the cells [12, 13] and vice versa.

Table 2 Relationship between the m6A methlylation extent and the transcript level in three organs of the Arabidopsis m6A transcriptome

We found that approximately 5.5 % of the m6A modified transcripts were extensively methylated by m6A in each organ (with fold change (ratio of MFPKM in the m6A-seq to FPKM in the mRNA-seq) ≥10, false discovery rate (FDR) <10−12, and read number per transcript ≥30). And 290 (2.1 %) extensively m6A methylated transcripts were found common in all of the three organs of leaves, flowers, and roots in Arabidopsis. These transcripts extensively modified by m6A were mainly associated with transporter, defenses or stress response, redox, nucleic acid binding, signal transduction, regulation of transcription, DNA/RNA/protein modification, cell proliferation or cycle, transposable element gene, pseudogenes, and snRNA (Table 3) [15–63]. Gene ontology (GO) analysis showed that the major molecular functions in these transcripts were responsible for molecule binding, transferase, hydrolase, stress responses, transporter, and kinase activity (Fig. 7).

Table 3 Potential functions in the 290 transcripts presenting extensively high methylation in all of the three organs of leaves, flowers and roots in Arabidopsis
Fig. 7
figure 7

GO analysis of the biological process (a) and molecular functions (b) for the 290 transcripts extensively methylated by m6A in the Arabidopsis transcriptome. Proportions in the figures indicated the level of the hits for each classification in the blast. As the majority of mRNA was methylated by m6A, methylation occurs on a functionally distinct subset of transcripts. However, most of the gene transcripts extensively methylated by m6A were related to the proteins that they were a direct transporter for energy molecules, for example, ATP or GTP, or accomplishment of their biological functions was dependent on these energy chaperones

Differential m6A methylation among leaves, flowers, and roots

High m6A sequencing depth (reads) of a transcript in the m6A-seq dataset may not suggest that this transcript was highly methylated by m6A in the cells because the gene transcript level of the counterpart may also be relatively high in the mRNA-seq dataset. To minimize effect of gene transcript level on estimation of m6A methylation extent of the transcripts and to more precisely evaluate differential m6A methylation among organs, we applied an algorithm ‘NFPKM’ to each transcript to estimate differential m6A methylation among three organs (see details in the Methods section of this paper). Two fold change and chi-square were applied for estimation of differential m6A methylation between organs using our algorithm. Accordingly, these two criteria were also used to estimate differential gene transcripts in mRNA-seq and to have a parallel comparison with the results of differential m6A methylation among organs in this study.

On average, 26.6 % of the transcripts presented differential in mRNA-seq (fold change of FPKM between two organs >2 or <0.5, and FDR <0.05), while 33.5 % showed differential methylation between two organs (fold change of NFPKM between two organs >2 or <0.5, and FDR <0.05) (Table 4). A paired analysis indicated that ratio of transcripts showing differential m6A methylation was significantly higher than that showing the differential gene transcripts in the three Arabidopsis organs (P <0.00035). The comparison also showed that the leaves had the highest extent of m6A methylation among three organs followed by the flowers. And the roots were less likely methylated among three organs (Table 4).

Table 4 The gene transcripts presenting differential transcript level and differential m6A methylation across three organs in Arabidopsis (fold change >2 or <0.5, FDR <0.05)

Eleven genes were randomly chosen for validation of our analysis above (Additional file 9). As products of qRT-PCR cover a short span in the transcriptome (50 to 150 bp) [13], two flanks of the amplicon containing one m6A peak in IGV program and showing differential m6A methylation were chosen to design primers (Additional file 9). The correlation coefficient between the qRT-PCR relative abundance results and the RIP-seq expected abundances was significant (r = 0.8632 (r between 0.7738-0.9526 for 95 % % CI), n = 33 genes, and P <10−4), indicating that our qRT-PCR data were consistent with the data estimated by m6A-seq and mRNA-seq using the IGV program (Additional file 10).

The heatmap representing the overall patterns of both differential transcript level and m6A methylation extent revealed completely distinct patterns: (1) between transcript level and m6A methylation; (2) among three organs; (3) among different genes; (4) among different regions on the same or different chromosomes; and (5) among five chromosomes (Fig. 8). This suggested that both differential gene transcripts and differential m6A methylation were highly heterogeneous in the Arabidopsis transcriptome and that regulation of the gene transcript level and extent of RNA m6A methylation may be relatively independent events.

Fig. 8
figure 8

Heatmap by the gene alphabetic order presenting the overall differential patterns of both gene transcript and m6A methylation level in the transcripts in Arabidopsis. a Overview of the differential patterns of m6A methylation between leaves and flowers. b Overview of the differential patterns of m6A methylation between leaves and roots. c Overview of the differential patterns of m6A methylation between flowers and roots. d Overview of the differential patterns of gene transcript level between leaves and flowers. e Overview of the differential patterns of gene transcript level between leaves and roots. f Overview of the differential patterns of gene transcript level between flowers and roots. The patterns in all of the six comparisons (a to f) above were based on the alphabetic order of the gene ID representing Chromosomes 1 to 5 (up to down)

Analysis of common elements between two replicate datasets showed that 2,628 (18.1 %) m6A modified transcripts in leaves, 1,920 (13.5 %) in flowers, and 1,166 (8.0 %) in roots showed a higher extent of methylation than the other two organs (fold change of NFPKM between two organs >2 or <0.5, FDR <0.05). Based on GO analysis, the major molecular functions in these transcripts were relevant to binding activity, transferase, hydrolase, kinase, transporter, and transcription factor (Fig. 9). Based on KEGG pathway analysis, certain transcripts presenting higher methylation in leaves than the other two organs were related to pathways of photosynthesis, carbohydrate, and nitrogen metabolism (Table 5); transcripts with higher methylation in flowers were concerning metabolic pathways of RNA degradation, DNA replication, and protein synthesis metabolisms (Table 5); transcripts presenting higher methylation in roots involved in biosynthesis of alkaloids, and carbonate metabolism (Table 5) (fold change of NFPKM between two organs >2 or <0.5, FDR <0.005). We found that 43 transcripts in leaves, 41 in flowers, and 23 in roots showed an extensively higher methylation level than the other two organs (fold change of NFPKM between two organs ≥10 or <0.1, FDR <10−10, and read number per transcript ≥20). The transcripts extensively methylated in leaves were mainly expressed for proteins located in mitochondria or chloroplast, photosynthesis, regulation of transcription, or stress response (Table 6) [64–71]; those in flowers (Table 7) were differentially or specifically expressed in the reproductive organs during flowering, or for cell proliferation, circadian rhythm, protein metabolism, transporter, or defense response [72–88]; and those in roots were mainly expressed for defense or stress response, transporters, redox process, or signaling transduction (Table 8) [89–100]. Therefore, the transcripts representing extensively higher m6A methylation in an organ revealed an intriguing phenomenon that functions of these transcripts were required for or highly related with unique biological roles of this organ (Tables 6, 7, and 8).

Fig. 9
figure 9

GO analysis of molecular functions of transcripts presenting a higher extent of m6A methylation in leaves (a), flowers (b), or roots (c) than other two organs in Arabidopsis

Table 5 KEGG pathways of the transcripts representing a higher methylation extent in one organ than the other two organs in Arabidopsis
Table 6 Potential functions of the 43 transcripts presenting extensively higher m6A methylation in leaves than that in the other two organs in Arabidopsis (fold change >10 or <0.1, FDR <10−10)
Table 7 Potential functions of the 41 transcripts presenting extensively higher m6A methylation in flowers than that in leaves and roots in Arabidopsis (fold change >10 or <0.1, FDR <10−10)
Table 8 Potential functions of the 23 transcripts presenting extensively higher m6A methylation in roots than that in leaves and flowers in Arabidopsis

Gene transcript level and m6A RNA methylation patterns in the transposable element genes

Interestingly, more than 97 % of the transcripts from the transposable element gene (TE) exhibited relatively high extent of m6A modification compared to gene transcript level (fold change of NFPKM between two organs >2 or <0.5, FDR <0.05). Another intriguing phenomenon was that most (>75 %) of the TE transcripts presented a ‘fragmented’ form in the both m6A-seq and mRNA-seq data (Fig. 10). The transcript level was distinct among the fragmented TE transcripts, though the transcript fragments were derived from the same TE in the genome (Fig. 10a). In some cases, the 5′ region of TEs was highly transcribed, but the other regions were much less expressed (Fig. 10a). Therefore, most of the TE transcripts remained a ‘fragmented’ form in the cells with a relatively high m6A extent and a low transcript level.

Fig. 10
figure 10

The transcriptional pattern (a) and m6A patterns (b) in the TE transcripts. The ‘fragmented form’ was observed in most of the TE transcript (the representative gene, AT5G35835, and trace files of two organs (leaves and flowers))

Discussion

Improvement of biotechnology is prerequisite for successful m6A-seq

RNA isolation from plants is sometimes difficult due to a thick cell wall which may be difficult to disrupt during RNA extraction [101, 102]. In addition, plant cells also produce secondary molecules (for example, polysaccharides), which are collected during RNA precipitation and can be problematic for RNA purification [101]. We modified the CTAB protocol for RNA isolation. The improved buffer showed a potential to efficiently disrupt cell walls, and specifically precipitate RNA, thus a large quantity of high purity RNA from Arabidopsis was obtained for mRNA-seq, m6A-seq, and input RNA-seq in this study (Additional file 11).

Another significant improvement in this study is that we optimized the protocol for m6A RNA immunoprecipitation to improve RIP efficiency and to minimize non-specific precipitation. We found that procedures including pre-precipitation through adding beads and IP buffer in the RNA solution before adding m6A antibody to exclude any non-specific binding RNA, and a vigorously washing of the binding beads after incubation with m6A antibody three to four times using IP buffer can significantly reduce the background (non-specific precipitation) in the RIP experiment. Incubation of the beads binding with RIP RNA in the dilution buffer at 50 Â°C for 90 min facilitated washing down of m6A RNA from the beads to improve m6A RNA yield.

m6A patterns between plant and mammal

The m6A patterns between plant (Figs. 4a and 5) and mammal were similar, for example, both rich near stop codons and 3′UTR [6, 8, 12, 13]. Both kingdoms had a similar consensus m6A methylation motif ‘RRm6ACH’ and a close frequency of m6A site per transcript in the transcriptome [6, 12, 13]. These phenomena indicated that m6A RNA methylation may be conserved between plant and mammal.

Some differences, however, were found between Arabidopsis and mammals as aforementioned (for example, Fig. 3a). These differences more likely represented the different m6A methylation extent and unique patterns between two species, or may result from the distinct depths of m6A-seq between two studies. Resolution in the m6A-seq with a low depth in the previous studies may be insufficient to discern the weak m6A signals in the coding regions in most of mRNA (Figs. 4a and 5). The previous report also suggested under-estimation of m6A sites in their studies [6, 8].

Dominant m6A enrichment in 5′UTR in Arabidopsis was reported by Luo et al. [13]. However, both our study from six Arabidopsis m6A-seq datasets (Fig. 5) and a recent report from rice m6A-seq of two different tissues [12] did not observer this dominant m6A pattern as the authors claimed a unique m6A pattern in the plant in their paper [13]. The m6A modified motif was confirmed by numerous studies highly conserved among the eukaryotes, including between two kingdoms of plant and animal [2, 6–8, 103]. A number of studies did not observe the dominant m6A enrichment in 5′UTR in mammals either [6, 8, 11]. The dominant m6A near the start codon reported by Luo et al. [13] may be an experimental artifact, for example probably caused by contamination from mRNA as a dominant peak was also found near the start codon in the overall pattern of mRNA-seq (Fig. 5). Nevertheless the difference in m6A enrichment pattern needs further confirmation.

The proportion of the m6A modified transcripts (over two-third) estimated in this study is relatively much higher than the previously reported [6, 12]. The previous studies with the robust experiments confirmed that the m6A modification is required for RNA stability and RNA transport from nucleus to cytoplasm [4, 6]. Otherwise, RNA will be degraded in the cells [3, 4, 6]. This suggested that the major proportion of mRNA in the cells, for example, over 70 % as estimated in our study, may rationally be modified by m6A [3]. Therefore, our estimation may reflect the genuine phenomenon in the cells.

Extensively high m6A methylation in certain transcripts may be suitable or required for the biological functions of these transcripts

Eukaryotic cells are highly compartmentalized and functionally differentiated [104–106]. mRNA translated into proteins for energy carriers (for example, ATP or GTP), transporter (for example, ion transporter), stress response, redox, protein post-translational modification, and protein located in mitochondria or chloroplast were found highly methylated by m6A in this study (Table 3). The common feature of these proteins was that they were a direct transporter for energy molecules, for example, ATP or GTP, or accomplishment of their biological functions was dependent on these energy chaperones [106, 107]. These transcripts may be transported for their protein translation near compartments for energy metabolism in the cells (for example, near mitochondria or chloroplast) or in the cells which were highly differentiated for energy process or stress response [109–111]. A long evolutionary history of m6A modification may imprint these transcripts with extensively high m6A methylation for their molecule stability [4, 112].

Extensive high m6A methylation was also found in the transcripts for some signaling factors, regulation factors (for example, regulation of DNA replication, RNA transcript, and protein synthesis and post-transcriptional process) and in certain sn(o)RNA (Table 3). These transcripts were highly or delicately regulated to maintain a low TL in the cells in most cases [107, 108, 113]. High methylation in these transcripts may endow RNA stability to these transcripts with a low TL, or confer signaling recognition to these transcripts [1, 3, 6, 112].

The transcripts of certain pseudogenes were extensively modified by m6A (Table 3). Accordingly, the transcripts of their relative functional genes had high methylation, for example, gene for NADH dehydrogenase (‘AT2G07709’) (Table 3). m6A modification proved a highly selective and regulative event [6, 8, 112]. If the transcript of a parent functional gene was extensively methylated by a regulative system, the transcript of the relative pseudogene may also be highly methylated by the same regulative system.

Dysfunction of m6A modification in the human cells resulted in an abnormal circadian rhythm [11]. Several transcripts responsible for regulation of circadian rhythm were highly methylated by m6A in Arabidopsis (Tables 3 and 7), suggesting that m6A may also play an important role in regulation of circadian rhythm in plants.

The m6A topological patterns and their potential functions in different types of RNAs

A frequency of 15 consensus (RRm6ACH) sequences per transcript should occur in the transcriptome [2, 14]. In deed, the occurrence of m6A methylation in the transcriptome was much lower than this expected frequency [2, 6, 114], suggesting that the majority of the consensus sequences were either not modified or many of them were demethylated promptly after the accomplishment of the biological metabolisms due to m6A modification [2]. As a result, methylation and demethylation may maintain a dynamic oscillation in the cells in response to environmental stimuli [2, 115, 116].

The extensively lower m6A signals in the coding regions including the splicing sites in most of mRNA discovered in this study (Figs. 4a and 5) could come from the sequencing bias. But the technical bias resulting in a large number (over 70 %) of the m6A modified mRNA containing this feature is impossible (Fig. 4a). The m6A methylation is required for a correct RNA splicing event [6, 116, 117]. Nevertheless, the most possible reason for the low signal detected in the coding regions in most of mRNA in this study may result from a prompt m6A demethylation event after accomplishment of the RNA splicing event [115, 116]. The consequent m6A demethylation in the coding regions may facilitate a speed of movement of ribosomes through mRNA chain thus may confer a high efficiency protein synthesis, while reservation of a great higher extent of methylation at the stop codon or 3′UTR may be responsible for RNA stability, signaling for transport and translocation, or as regulatory elements for protein translation through the recruitment of specific factors onto these m6A sites for RNA transport or protein synthesis [4, 112, 118, 119].

Dominant m6A enrichment near stop codons and 3′UTR (as shown in Fig. 4a) was observed in most (over 70 %) of mRNA in this study as previously reported [6, 8, 12, 13]. This m6A distributing type (Fig. 4a) may represent the typical m6A topological pattern in most of the mature mRNA. However, a small proportion of mRNA did not present this m6A enrichment in 3′UTR (as shown in Fig. 4b). m6A methylation confers RNA stability [4]. Non-dominant m6A enrichment in 3′UTR in this small proportion of mRNA may reflect a phenomenon of ‘subsequent-demethylation’ to the m6A dominant peaks in 3′UTR or stop condon, further suggesting that these RNA may undergo a process of degradation due to m6A demethylation [4]. Therefore, different m6A topological patterns in mRNA may reflect RNA status or fate in the cells, for example, pre-mature RNA, mature RNA, or RNA being in degradation.

Both rRNA and tRNA were mainly modified by numerous types of cytosine methylation [1, 3]. However, a relatively high m6A methylation was found in all rRNAs (Fig. 6a and b) and in a certain proportion of tRNAs (approximately 10 %) (Fig. 6c) in Arabidopsis in this study. Mature rRNA and tRNA were derived from pre-rRNA and pre-tRNA through a series of complex biological and molecular processes, including splicing and folding. The splicing machinery in rRNA and tRNA was similar to that in mRNA in eukaryotes [3]. m6A may guide a correct splicing event in these types of RNAs as the role of the m6A methylation required for the correct splicing events in mRNA [6]. m6A may not undergo a subsequent demethylation event after splicing for the mature rRNA, but demethylation of m6A after splicing may be required for the mature tRNA as hypothetically aforementioned for most of mRNA. This would be the reason that m6A methylation was reserved and observed in all rRNAs (Fig. 6a and b), but it was detected in a small proportion of tRNAs (Fig. 6c).

Therefore, diverse m6A patterns and topologies may be unique or required for the miscellaneous functions in the different types of RNAs [3, 112].

Potential roles of differential m6A methylation among plant organs

Differential gene expression among plant organs has proved responsible for organ differentiation and development [120, 121]. Differential level of m6A methylaytion among three organs was much higher than that in gene transcript level (Table 4), suggesting that m6A modification may be another important contributor for organ differentiation or maintenance of differential status among the organs in Arabidopsis.

The transcripts presenting extensively higher m6A modification in an organ than other organs showed a connection of the functions of these transcripts required for or specific to this organ. For example, the transcripts presenting an extensively higher level of m6A methylation in leaves were related to photosynthesis metabolism or proteins located in mitochondria or chloroplast, and those in roots showed response to stresses, redox process, and transporters (Tables 6 and 8). However, most of the differential and extensively methylated transcripts in flowers were related to regulation of reproductive organ development, stress response, cell proliferation, differentiation, or circadian rhythm (Table 7) [18, 122]. Dysfunction of METTL3, the gene responsible for m6A modification, can result in an arrest of the early development in embryo at the globular stage in Arabidopsis [7], suggesting that m6A methylation may play an important role in differentiation and development of the reproductive organs in plants.

The proportion of the m6A modified transcripts was the highest in the roots among three organs (Additional file 3). The surrounding environment for root growth is more complex than that for leaves and flowers [123, 124], which may require a higher proportion of the m6A transcripts in roots to adapt to more diverse conditions [2–4]. However, the overall extent for m6A RNA methylation showed the highest in the leaves among three organs (Table 1). The leaves are the major organ in plants responsible for photosynthesis metabolism and energy transition [123]. This may require a higher extent of m6A RNA methylation in the leaves to fit for these metabolisms that were processed in severe conditions in most of cases if our hypothesis aforementioned is rational.

Potential biological significance of gene transcript and m6A patterns in TE

Most of the TE transcripts (>85 % in this study) remained a relatively low level in the cells. This may repress a neo-transposition of TEs into the genome [125]. 17 % to 85 % of the genomic sequences are composed of TE in the higher plant species [126, 127]. In most of cases, the neo-transposition of TEs was highly repressed because it is often detrimental to the host due to a prompt expansion of the genomic size or dysfunction of the functional genes by this event [128]. The fragmented TE transcripts may not be reversely transcribed into a complete cDNA, thus the malignant TE transposition may be consequently avoided. In addition, the distribution of these unconnected transcripts in the cells may play a role similar to the small interfering RNA (siRNA), which may in turn repress the transcription or the neo-transposition of these TEs [125, 129]. The remnant TE transcript fragments may also be induced in response to stress, which may further trigger a serial of complex reactions [127]. Some of the fragmented TE transcripts may play a role as ‘a regulator’ for other gene expression or as a direct regulator for the nearby genes [125, 127, 130, 131], or may work as powerful regulators of the immune response as the functions recently discovered in some miRNA [132]. Because the fragmented transcripts were more likely to be degraded than the intact transcripts [3], relatively high m6A methylation in the fragmented TE transcripts may prevent these segmented transcripts from further degradation so that a relatively low level of these transcripts can be maintained in the cells [3].

Conclusions

Thanks to significant improvements of technologies for RIP experiments in this study, high resolution of transcriptome-wide mapping of m6A was available in Arabidopsis. This is the first study for comprehensive characterization of m6A patterns of different types of RNAs, relationship between m6A methylation extent and gene transcript level, and differential features of m6A methylation among three plant organs. Two-third of the transcripts were modified by m6A in Arabidopsis. 35,000 to 48,000 m6A sites and approximately 1.4 to 2.0 m6A sites per transcript were mapped to the Arabidopsis genome.

Over 85 % of the m6A modified transcripts had a relatively high m6A methylation level (P <0.05), while <4 % had a low m6A extent compared with their transcript level (P <0.05). Approximately 5.5 % of the methylated transcripts presented extensively high methylation (fold change >10 compared to TL, and FDR <10−12). The highly expressed transcripts were relatively less methylated by m6A and vice versa (P <0.001). The 290 (2.1 %) highly methylated transcripts were mainly expressed for stress response, redox, signaling factors, regulation factors, and some ncRNA. Most of the biological functions in these transcripts were involved in molecule binding, transferase, transporter, and kinase activity.

Most of mRNA (over 70 %) was characterized by a typical m6A topology, that is, one or two predominant peaks at the stop codon or 3′UTR accompanying with extensively low m6A signals in the coding regions. Unlike a recent observation of another predominant m6A enrichment in the 5′ mRNA in Arabidopsis [13], we found that m6A predominantly distributed only at the stop codon or 3′UTR. Some sn(o) RNAs were also highly methylated with a single m6A site in these transcripts. All rRNAs was relatively highly methylated by one or several m6A sites. Non- or slight-m6A methylation was observed in most of tRNAs (approximately 90 %), and the remaining tRNAs were relatively hyper-methylated by a single m6A site. Interestingly, most (over 75 %) of the transcribed TEs maintained a relatively high m6A methylation. Therefore, the topologies of m6A in different RNAs not only confer diverse m6A patterns in the cell, the unique m6A pattern in a specific RNA may endow special functions to this RNA.

The similar m6A patterns between plant and mammal suggested that m6A methylation may be conserved between two living kingdoms. Differences were also found between plant and mammal, which may represent their unique m6A patterns in the two living kingdoms.

Proportion of transcripts (33.5 %) showing differential m6A methylation among three organs was greatly higher than that (22.6 %) presenting differential transcript level in Arabidopsis (P <0.00035). Function of the transcripts with extensively higher methylation in an organ than others were required or suitable for unique biological roles of this organ. Therefore, m6A methylation may be an important contributor to the organ differentiation and may confer unique functions to this organ.

Methods

Plant materials

Wild type, Columbia ecotype (Col-0), of the mouseear cress (Arabidopsis thaliana) was used in this study. The plants were grown in the greenhouse under a photoperiod of 16 h light/8 h dark at 22–24 Â°C. When the plants were at the blooming period (5 weeks after seed germination), the plant materials of the flowers, rosette leaves, and roots were separately collected, treated with liquid nitrogen and stored at −80 Â°C until use.

RNA isolation and purification

All centrifuge tubes and pipette-tips are RNase-free or must be treated with DEPC. All buffers were RNase free or prepared using DEPC-treated dd H2O. The modified CTAB buffer was used for RNA isolation [133]. In brief, approximately 10 g of the frozen plant materials were ground into fine powder in liquid nitrogen. The plant powder was promptly transferred into four 50 mL tubes with 25 mL CTAB buffer (2 % CTAB, 1.0 % PVP-40, 2.0 M NaCl, 100 mM Tris, 25 mM EDTA-Na2, and 1.0 % β-mercaptoethanol) in each tube. Incubate the tubes at 65 Â°C for 10 min and invert tubes for several seconds every minute during incubation. Add 200 Î¼L chloroform to each tube, invert quickly 150–200 times. Centrifuge tubes at 4,500 g for 2 min at 4 Â°C, and discard most of the chloroform solution using a pipette-tip. Centrifuge tubes at 15,000 g for 5 min at 4 Â°C. Transfer the supernatant to new tubes. Add 8.0 mL LiCl (8.0 M) into each tube and mix well. Store the mixture at −20 Â°C for 30 min, then centrifuge at 17,000 g for 30 min at 4 Â°C. Use 2 mL 80 % ethanol to rinse RNA pellet. Centrifuge at 17,000 g for 2 min, discard the ethanol, and then dry tubes in the laminar-flow hood. Dissolve the dried RNA pellet using 400 Î¼L RNase-free water and treat with DNase (Promega, Madison, WI, USA) to remove DNA contamination. Add 300 Î¼L chloroform-phenol (1:1, v/v) into the tube, invert quickly 200 times. Centrifuge tubes at 13,000 g for 1 min at 4 Â°C, and discard most of the chloroform-phenol solution in the hood. Centrifuge tubes at 15,000 g for 5 min at 4 Â°C. Transfer the supernatant to new tubes. Add 2.8 volume of ethanol and 0.10–0.15 volume of NaAc (pH 5.6), and mix well. Store the mixture at −20 Â°C for 30 min, then centrifuge at 15,000 g for 30 min at 4 Â°C. Store the RNA pellet in the ethanol-NaAc solution at −80 Â°C until use.

RNA fragmentation

The purified total RNA was diluted in the fragmentation buffer (10 mM ZnCl2 and 10 mM Tris–HCl, pH 7) with the RNA final concentration of approximately 1.0 Î¼g Î¼L−1 [6, 10]. The diluted RNA was fragmented into approximately 100-nucleotide-long by incubation at 95 Â°C for 5 min. The chemical fragmentation reaction was stopped with 0.05 M EDTA. The ethanol-NaAc solution was used to precipitate RNA. The fragmented RNA was prepared for RIP and m6A-seq.

RNA immunoprecipitation (RIP)

Approximately 2.5 mg fragmented total RNA was resuspended in the 1,000 Î¼L IP buffer (150 mM NaCl, 0.1 % NP-40, 10 mM Tris–HCl, pH 7.4, 2 mM RVC (Sigma-Aldrich, St. Louis, MO, USA), 200 U RNasin (Promega, Madison, WI, USA), and 0.5 mg mL− 1 BSA). Non-specific binding RNA was pre-precipitated by adding 3.0 Î¼L Protein A bead (Life Technology, Grand Island, NY, USA) and incubation on a rotating wheel for 2 h. The magnet rack was used to precipitate the beads. The supernatant was transferred to a new tube. Add 5.0 Î¼g m6A-specific monoclonal antibody (Merck Millipore, Billerica, MA, USA) into the tube containing RNA and IP buffer, and incubate on a rotating wheel for 3–4 h. Then add 5.0 Î¼L Protein A bead (Life Technology, Grand Island, NY, USA) for an additional rotation for 2 h. The magnet rack was used to precipitate the beads. The beads were vigorously washed using 1,000 Î¼L IP buffer three to four times. Discard the IP buffer. Add 300 Î¼L dilution buffer (10 mM Tris–HCl pH 7.5) into the bead tube and incubate at 50 Â°C for 90 min. Precipitate the beads using the magnet rack and transfer the supernatant to a new tube. The ethanol-NaAc solution and glycogen were used to precipitate m6A RNA. The m6A RNA pellet was washed using 80 % ethanol and then resuspended into 15 Î¼L dd-H2O for m6A-seq and HPLC-MS/MS analysis.

mRNA-seq, m6A-seq, and RNA-seq from the input samples

The quality control (QC) tests for the RNA samples were performed using Agilent Technologies (Santa Clara, CA, USA). Library for mRNA-seq was generated using mRNA sequencing kit plus random primers. High throughput m6A-seq, mRNA-seq, and input RNA-seq of three samples of leaves, flowers, and roots was performed on HiSeq 2000 (Illumina Inc., San Diego, CA, USA) at Purdue University Genomics Core Facility (http://www.genomics.purdue.edu/services/core.shtml).

Approximately 10.0 Î¼g purified total RNA was reserved for mRNA-seq before the RNA samples were used for fragmentation. And approximately 2.5 Î¼g fragmented total RNA was reserved and used for the input RNA-seq before the RIP experiments. Thus, mRNA-seq, m6A-seq, and input RNA-seq were parallel and their data were mutually comparable in this study [6]. RNA integrity number (RIN) was estimated using a Nanodrop 2000 UV vis (Thermo Fisher Scientific, Wilmington, NC, USA). The QC tests were done by Agilent Technologies (Santa Clara, CA, USA). All RNA sequencing of three samples of leaves, flowers, and roots was performed on the same sequencer at the same batch.

HPLC-MS/MS

A total of 100–200 ng of input total RNA or m6A RNA from RIP experiment was digested by 2 U nuclease P1 (US Biological Life Science, Salem, MA, USA) at 37 Â°C for 2 h, and 0.5 U alkaline phosphatase (Promega, Madison, WI, USA) at 37 Â°C for an additional 2 h. A total of 5 Î¼L of the digested and purified solutions were assayed by HPLC-MS/MS. Nucleosides were separated by reversed-phase high performance liquid chromatography, using a Waters Xterra C18 column (2.1 × 150 mm, 3.5 Î¼m) with a water/acetonitrile gradient. Mass spectrometry was performed using an Agilent 6460 triple quadrupole mass spectrometer in positive electrospray ionization mode. Mass transitions were 268.1 -- > 136.1 for A and 282.1 -- > 150.1 for m6A (Additional file 2).

Alignment of reads and visualization of m6A peaks

Both RNA-seq and m6A-seq datasets were mapped to the Arabidopsis genome (TAIR10) using TopHat2 with a parameter of ‘-b2-fast’ [134]. The potential PCR duplicates were removed by the parameter ‘rmdup’ rooted in SAMtools [135]. The fragment numbers for each transcript were estimated using the featureCounts with a parameter of ‘-p’ [136].

The peaks and distributing patterns of m6A in the Arabidopsis transcriptome were visualized using free software, Integrative Genomics Viewer (IGV2.3, Boston, MA, USA [137]).

We tried several previously published protocols, for example, moving-window, to call m6A peaks based on comparison of m6A-seq with input RNA-seq counterpart. And we found all of these methods failed recognition of most (>92 %) of m6A sites in CDS (Additional file 12). We finally did not apply these methods in this study. Because of the extensively low non-specific immunoprecipitation rate (<1 %) in this study (Additional file 2), all the mapped reads in the m6A-seq were considered to be sequenced and derived from the specific immunoprecipitation of the m6A RNA fragments. Thus, an estimation of m6A peak number of a m6A modified transcript was calculated by this formula: total mapped absolute length covered by m6A fragments within the transcript/150, considered that library construction for m6A-seq was created from a m6A RNA pool with an average RNA length of approximately 106 nucleotides (Additional file 11: Fig. S3b) and average coverage of a peak in the m6A-seq data base was approximately 150-nucleotide long in this study as visualized by IGV 2.3.

Discernment of m6A topological patterns

Distribution of m6A sites in the different regions of the transcripts was estimated by Dominissini et al.’s method [6]. The consensus m6A motif sequences were figured out by Luo et al.’s protocol [13] with modification: approximately 1,000 the highest m6A peaks and approximately 100 nt length around each m6A peak were used for deduction of the consensus m6A motif sequences. The typical m6A patterns of different types of RNA were drawn from the m6A mapping analysis based on visualization using the IGV2.3 program [13].

The overall m6A distributing patterns were discerned by this method: a gene was splitted in to 60 bins. The read depth of each bin was normalized by per 1 kb per 1 Mb data, then the normalized depth was used to plot the patterns.

Comparison of m6A methylation extent versus transcript level

The sequenced fragment number of each transcript in mRNA-seq was normalized using the algorithm of Fragments Per Kilobase of Transcript Per Million Fragments Mapped (FPKM = Counts of mapped fragments × 109)/(Length of transcript × Total count of the mapped fragments); ‘fragment’ refers to mapped reads after removal of PCR duplicates [138]. While the sequenced fragment number of each transcript in m6A-seq was normalized using a modified FPKM (MFPKM = Counts of mapped fragments × 109)/(Total absolute mapped length on the chromosome covered by m6A fragments within the transcript × Total count of the mapped m6A fragments)), considered that the library in m6A-seq was derived only from RNA fragments containing m6A sites, and not from the entire transcript [138].

The m6A methylation extent of a transcript were categorized into three groupings based on comparison of MFPKM of the transcript in the m6A-seq with the FPKM of the same transcript in the mRNA-seq using χ 2 test: (1) m6A methylation extent ‘equivalent’ to transcript level (‘equivalent’, ratio of FPKM to MFPKM fits 1:1 (P <0.05)); (2) methylation extent higher than transcript level (‘Hi’, ratio of FPKM to MFPKM < 1 (P <0.05)); and (3) methylation extent lower than transcript level (‘Low’, ratio of FPKM to MFPKM >1 (P <0.05)). The chi-square value was determined by this formula: χ 2 = (MFPKM*100/(MFPKM + FPKM) - 50)^2/50 + (FPKM*100/(MFPKM + FPKM) - 50)^2/50.

Analysis of both differential transcript level and differential m6A methylation among plant organs

mRNA-seq data were normalized by FPKM as described above. χ 2 tests were used to estimate whether FPKM was significantly different between two organs using R 3.1 (http://cran.r-project.org/bin/windows/base/). The transcripts with fold change in FPKM >2.0 or <0.5, and FDR <0.02 were considered differentially expressed between two organs [6, 10].

To minimize influence of transcript level on estimation of differentiation of m6A extent, m6A-seq data were normalized by a specific algorithm, NFPKM (NFPKM = MFPKM in m6A-seq/LOG (FPKM in mRNA-seq, 2)). χ 2 tests were also used to estimate whether NFPKM of a m6A modified transcript was significantly different between two organs using R 3.1. The transcripts with fold change of NFPKM >2.0 or <0.5, and FDR <0.005 were considered differentially methylated between two organs [6, 10]. The transcripts commonly presenting higher m6A extent in an organ in the two replicates were used for gene ontology (GO) analysis. P values of the chi-square tests to all gene transcripts of a replicate were input in R 3.1.2 (‘Pumpkin Helmet’) to estimate FDR of each transcript under a parameter of ‘p.adjust (P value,method = ‘fdr’)’. The fold change is a fixed threshold for all classifications, thus FDR thresholds may vary with different analyses.

Visualization of the overall differential patterns

The overall patterns of differential transcript level of all transcripts were discerned using ratio of FPKM of mRNA-seq between two organs. Similarly, the overall patterns of differential m6A methylation of all transcripts were depicted using ratio of NFPKM derived from both m6A-seq and mRNA-seq (see above) between two organs. The heatmap representing the overall differential patterns of gene transcript and m6A methylation was created using R3.1 (http://cran.r-project.org/bin/windows/base/).

GO and KEGG pathway analysis

The GO results and potential molecular functions in the transcripts were deducted from the online tool released in the TAIR website (http://www.arabidopsis.org/) [139]. And molecular functions in some transcripts were inferred from the recent publications.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in certain transcripts of interests were figured out using the online tool released in the DAVID bioinformatics resources (http://david.abcc.ncifcrf.gov/).

qRT-PCR

Quantitative real-time PCR (qRT-PCR) was performed to assess the relative abundance of m6A RNA in the RIP samples. All purified RNA templates were transferred into cDNA using Quanta qScriptâ„¢ cDNA Synthesis Kits (Quanta BioSciences, Inc., Gaithersburg, MD, USA). Eleven genes were randomly chosen for this test (Additional file 9). The qRT-PCR primers were designed to span exon-exon junctions in order to eliminate the potential amplification of genomic DNA and non-spliced mRNA. qRT-PCR was performed on C1000 Thermal Cycler (Bio-RAD) using SYBR Green SuperMix buffer (Bio-RAD) and 300 ng total cDNA template for amplification. Because the qRT-PCR amplicon spanned an exon-exon junction with a length of 80–150 bp and the m6A enrichment in the region around splicing sites was usually low in most of mRNA in this study (Fig. 4a), cDNA from total RNA of the Actin2 gene was used for housekeeping gene to estimate the relative abundance (RA) of m6A RNA in the qRT-PCR amplicons based on this algorithm: RA = 100 × 2−ΔC. The expected abundance (EA) of m6A RNA in the m6A-seq dataset was estimated by this algorithm: EA = 100× (the mapped m6A RNA reads of the test gene in m6A-seq and in the region for qRT-PCR test/ the mapped RNA reads of the Actin2 gene in mRNA-seq and in the region for qRT-PCR test). The consistency between the AR and ER patterns among three organs was compared (Additional file 10).

Ethics statement

The plant materials used in this study are freely used and available to all researchers without any protection for intellectual property right. This research meets all applicable standards for the ethics of experimentation and research integrity from all five institutes that provide supports to this study.

Data access

The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE72706.

Abbreviations

ATP:

adenosine triphosphate

BSA:

bovine serum albumin

cDNA:

complementary DNA

CDS:

coding sequence

CI:

confidence intervals

CTAB:

cetyl trimethylammonium bromide

DEPC:

diethylpyrocarbonate

EA:

expected abundance

EDTA:

ethylenediaminetetraacetic acid

FDR:

false discovery rate

FPKM:

fragments per kilobase of transcript per million fragments mapped of the transcript in mRNA-seq

GO:

gene ontology

GTP:

guanosine triphosphate

HPLC:

High-performance liquid chromatography

KEGG:

the Kyoto Encyclopedia of Genes and Genomes

m6A:

N 6-methyladenosine

m6A-seq:

RNA sequencing based on m6A RNA immunoprecipitation

MFPKM:

modified fragments per kilobase of transcript per million fragments mapped of the transcript in m6A-seq

mRNA-seq:

RNA sequencing from oligo-dT RNA

MS:

mass spectrometry

NADH:

nicotinamide adenine dinucleotide

NFPKM:

MFPKM in m6A-seq divided by LOG (FPKM in mRNA-seq, 2)

PVP:

polyvinylpyrrolidone

QC:

quality control test

qRT-PCR:

quantitative real-time polymerase chain reaction

RA:

relative abundance

RIN:

RNA integrity number

RIP:

RNA immunoprecipitation

sn(o)RNA:

small nuclear RNA and small nucleolar RNA

TAIR:

the Arabidopsis Information Resource (https://www.arabidopsis.org/)

TE:

transposable element gene

TL:

transcript level

Tris:

tris(hydroxymethyl)aminomethane

TSS:

transcription start site

U:

unit

UTR:

untranslated region

References

  1. Matera AG, Terns RM, Terns MP. Non-coding RNAs: lessons from the small nuclear and small nucleolar RNAs. Nat Rev Mol Cell Biol. 2007;8:209–20.

    Article  CAS  PubMed  Google Scholar 

  2. Pan T. N 6-methyl-adenosine modification in messenger and long non-coding RNA. Trends Biochem Sci. 2013;4:204–9.

    Article  CAS  Google Scholar 

  3. Grosjean H. Fine-tuning of RNA functions by modification and editing. Berlin: Springer-Verlag; 2005.

    Google Scholar 

  4. Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, et al. N 6-methyladenosine-dependent regulation of messenger RNA stability. Nature. 2013;7481:117–20.

    Article  CAS  Google Scholar 

  5. Bodi Z, Bottley A, Archer N, May ST, Fray RG. Yeast m6A methylated mRNAs are enriched on translating ribosomes during meiosis, and under rapamycin treatment. PLoS One. 2015;10(7):e0132090.

    Article  PubMed Central  PubMed  CAS  Google Scholar 

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

    Article  CAS  Google Scholar 

  7. Zhong S, Hongying L, Bodi Z, Button JD, Vespa L, Herzog M, et al. MTA is an Arabidopsis messenger RNA adenosine methylase and interacts with a homolog of a sex-specific splicing factor. Plant Cell. 2008;20:1278–88.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  8. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3' UTRs and near stop codons. Cell. 2012;7:1635–46.

    Article  CAS  Google Scholar 

  9. Keene JD, Komisarow JM, Friedersdorf MB. RIP-Chip: the isolation and identification of mRNAs, microRNAs and protein components of ribonucleoprotein complexes from cell extracts. Nat Protoc. 2006;1:302–7.

    Article  CAS  PubMed  Google Scholar 

  10. Dominissini D, Moshitch-Moshkovitz S, Salmon-Divon M, Amariglio N, Rechavi G. Transcriptome-wide mapping of N 6-methyladenosine by m6A-seq based on immunocapturing and massively parallel sequencing. Nat Protoc. 2013;1:176–89.

    Article  CAS  Google Scholar 

  11. Fustin JM, Doi M, Yamaguchi Y, Hida H, Nishimura S, Yoshida M, et al. RNA-methylation-dependent RNA processing controls the speed of the circadian clock. Cell. 2013;155:793–806.

    Article  CAS  PubMed  Google Scholar 

  12. Li YL, Wang XL, Li CP, Hu SN, Yu J, Song SH. Transcriptome-wide N6-methyladenosine profiling of rice callus and leaf reveals the presence of tissue-specific competitors involved in selective mRNA modification. RNA Biol. 2014;11(9):1180–8.

    Article  PubMed  Google Scholar 

  13. Luo GZ, MacQueen A, Zheng G, Duan H, Dore LC, Lu Z, et al. Unique features of the m6A methylome in Arabidopsis thaliana. Nat Commun. 2014;5:5630.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  14. Csepany T, Lin A, Baldick Jr CJ, Beemon K. Sequence specificity of mRNA N6-adenosine methyltransferase. J Biol Chem. 1990;265:20117–22.

    CAS  PubMed  Google Scholar 

  15. Schledzewski K, Brinkmann H, Mendel RR. Phylogenetic analysis of components of the eukaryotic vesicle transport system reveals a common origin of adaptor protein complexes 1, 2, and 3 and the F subcomplex of the coatomer COPI. J Mol Evol. 1999;6:770–8.

    Article  Google Scholar 

  16. Footitt S, Slocombe SP, Larner V, Kurup S, Wu Y, Larson T, et al. Control of germination and lipid mobilization by COMATOSE, the Arabidopsis homologue of human ALDP. EMBO J. 2002;12:2912–22.

    Article  Google Scholar 

  17. Bedhomme M, Hoffmann M, McCarthy EA, Gambonnet B, Moran RG, Rébeillé F, et al. Folate metabolism in plants: an Arabidopsis homolog of the mammalian mitochondrial folate transporter mediates folate import into chloroplasts. J Biol Chem. 2005;41:34823–31.

    Article  CAS  Google Scholar 

  18. Heyndrickx KS, Vandepoele K. Systematic identification of functional plant modules through the integration of complementary data sources. Plant Physiol. 2012;159:884–901.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  19. Jammes F, Song C, Shin D, Munemasa S, Takeda K, Gu D, et al. MAP kinases MPK9 and MPK12 are preferentially expressed in guard cells and positively regulate ROS-mediated ABA signaling. Proc Natl Acad Sci U S A. 2009;48:20520–5.

    Article  Google Scholar 

  20. Aranda-Sicilia MN, Cagnac O, Chanroj S, Sze H, Rodríguez-Rosales MP, Venema K. Arabidopsis KEA2, a homolog of bacterial KefC, encodes a K(+)/H(+) antiporter with a chloroplast transit peptide. Biochim Biophys Acta. 2012;9:2362–71.

    Article  CAS  Google Scholar 

  21. Tsugeki R, Ditengou FA, Sumi Y, Teale W, Palme K, Okada K. NO VEIN mediates auxin-dependent specification and patterning in the Arabidopsis embryo, shoot, and root. Plant Cell. 2009;10:3133–51.

    Article  CAS  Google Scholar 

  22. Charrier B, Champion A, Henry Y, Kreis M. Expression profiling of the whole Arabidopsis shaggy-like kinase multigene family by real-time reverse transcriptase-polymerase chain reaction. Plant Physiol. 2002;2:577–90.

    Article  CAS  Google Scholar 

  23. Lellis AD, Kasschau KD, Whitham SA, Carrington JC. Loss-of-susceptibility mutants of Arabidopsis thaliana reveal an essential role for eIF(iso)4E during potyvirus infection. Curr Biol. 2002;12:1046–51.

    Article  CAS  PubMed  Google Scholar 

  24. Kawamura Y, Uemura M. Mass spectrometric approach for identifying putative plasma membrane proteins of Arabidopsis leaves associated with cold acclimation. Plant J. 2003;2:141–54.

    Article  CAS  Google Scholar 

  25. Sarry JE, Kuhn L, Ducruix C, Lafaye A, Junot C, Hugouvieux V, et al. The early responses of Arabidopsis thaliana cells to cadmium exposure explored by protein and metabolite profiling analyses. Proteomics. 2006;7:2180–98.

    Article  CAS  Google Scholar 

  26. Teotia S, Lamb RS. The paralogous genes RADICAL-INDUCED CELL DEATH1 and SIMILAR TO RCD ONE1 have partially redundant functions during Arabidopsis development. Plant J. 2008;1:29–41.

    Google Scholar 

  27. Chen LT, Luo M, Wang YY, Wu K. Involvement of Arabidopsis histone deacetylase HDA6 in ABA and salt stress response. J Exp Bot. 2010;12:3345–53.

    Article  CAS  Google Scholar 

  28. Kindgren P, Kremnev D, Blanco NE, de Dios Barajas López J, Fernández AP, Tellgren-Roth C, et al. The plastid redox insensitive 2 mutant of Arabidopsis is impaired in PEP activity and high light-dependent plastid redox signalling to the nucleus. Plant J. 2012;2:279–91.

    Article  CAS  Google Scholar 

  29. Yang S, Zeng X, Li T, Liu M, Zhang S, Gao S, et al. AtACDO1, an ABC1-like kinase gene, is involved in chlorophyll degradation and the response to photooxidative stress in Arabidopsis. J Exp Bot. 2012;63(10):3959–73.

    Article  CAS  PubMed  Google Scholar 

  30. Zhang X, Wang C, Zhang Y, Sun Y, Mou Z. The Arabidopsis mediator complex subunit positively regulates salicylate-mediated systemic acquired resistance and jasmonate/ethylene-induced defense pathways. Plant Cell. 2012;10:4294–309.

    Article  CAS  Google Scholar 

  31. Pu L, Liu MS, Kim SY, Chen LF, Fletcher JC, Sung ZR. EMBRYONIC FLOWER1 and ULTRAPETALA1 act antagonistically on Arabidopsis development and stress response. Plant Physiol. 2013;2:812–30.

    Article  CAS  Google Scholar 

  32. Shrestha B, Guragain B, Sridhar VV. Involvement of co-repressor LUH and the adapter proteins SLK1 and SLK2 in the regulation of abiotic stress response genes in Arabidopsis. BMC Plant Biol. 2014;14:54.

    Article  PubMed Central  PubMed  CAS  Google Scholar 

  33. Lee SJ, Lee MH, Kim JI, Kim SY. Arabidopsis putative MAP kinase kinase kinases Raf10 and Raf11 are positive regulators of seed dormancy and ABA response. Plant Cell Physiol. 2015;1:84–97.

    Article  Google Scholar 

  34. Zolman BK, Silva ID, Bartel B. The Arabidopsis pxa1 mutant is defective in an ATP-binding cassette transporter-like protein required for peroxisomal fatty acid beta-oxidation. Plant Physiol. 2001;3:1266–78.

    Article  Google Scholar 

  35. Nito K, Kamigaki A, Kondo M, Hayashi M, Nishimura M. Functional classification of Arabidopsis peroxisome biogenesis factors proposed from analyses of knockdown mutants. Plant Cell Physiol. 2007;6:763–74.

    Article  CAS  Google Scholar 

  36. Attallah CV, Welchen E, Martin AP, Spinelli SV, Bonnard G, Palatnik JF, et al. Plants contain two SCO proteins that are differentially involved in cytochrome c oxidase function and copper and redox homeostasis. J Exp Bot. 2011;12:4281–94.

    Article  CAS  Google Scholar 

  37. Hare PD, Moller SG, Huang LF, Chua NH. LAF3, a novel factor required for normal phytochrome A signaling. Plant Physiol. 2003;4:1592–604.

    Article  CAS  Google Scholar 

  38. Parker G, Schofield R, Sundberg B, Turner S. Isolation of COV1, a gene involved in the regulation of vascular patterning in the stem of Arabidopsis. Development. 2003;10:2139–48.

    Article  CAS  Google Scholar 

  39. van Leeuwen W, Okrész L, Bögre L, Munnik T. Learning the lipid language of plant signalling. Trends Plant Sci. 2004;8:378–84.

    Article  CAS  Google Scholar 

  40. Caryl AP, Armstrong SJ, Jones GH, Franklin FC. A homologue of the yeast HOP1 gene is inactivated in the Arabidopsis meiotic mutant asy1. Chromosoma. 2000;109(1–2):62–71.

    Article  CAS  PubMed  Google Scholar 

  41. Kanno T, Bucher E, Daxinger L, Huettel B, Böhmdorfer G, Gregor W, et al. A structural-maintenance-of-chromosomes hinge domain-containing protein is required for RNA-directed DNA methylation. Nat Genet. 2008;40(5):670–5.

    Article  CAS  PubMed  Google Scholar 

  42. Campi M, D'Andrea L, Emiliani J, Casati P. Participation of chromatin-remodeling proteins in the repair of ultraviolet-B-damaged DNA. Plant Physiol. 2012;158(2):981–95.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  43. Liu Y, Deng Y, Li G, Zhao J. Replication factor C1 (RFC1) is required for double-strand break repair during meiotic homologous recombination in Arabidopsis. Plant J. 2013;73:154–65.

    Article  CAS  PubMed  Google Scholar 

  44. Hanada K, Zhang X, Borevitz JO, Li WH, Shiu SH. A large number of novel coding small open reading frames in the intergenic regions of the Arabidopsis thaliana genome are transcribed and/or under purifying selection. Genome Res. 2007;5:632–40.

    Article  CAS  Google Scholar 

  45. Habu Y, Mathieu O, Tariq M, Probst AV, Smathajitt C, Zhu T, et al. Epigenetic regulation of transcription in intermediate heterochromatin. EMBO Rep. 2006;12:1279–84.

    Article  CAS  Google Scholar 

  46. Yaish MW, Peng M, Rothstein SJ. AtMBD9 modulates Arabidopsis development through the dual epigenetic pathways of DNA methylation and histone acetylation. Plant J. 2009;1:123–35.

    Article  CAS  Google Scholar 

  47. Joukov V, Groen AC, Prokhorova T, Gerson R, White E, Rodriguez A, et al. The BRCA1/BARD1 heterodimer modulates ran-dependent mitotic spindle assembly. Cell. 2006;3:539–52.

    Article  CAS  Google Scholar 

  48. Liu Y, Wang F, Zhang H, He H, Ma L, Deng XW. Functional characterization of the Arabidopsis ubiquitin-specific protease gene family reveals specific role and redundancy of individual members in development. Plant J. 2008;5:844–56.

    Article  CAS  Google Scholar 

  49. Silady RA, Ehrhardt DW, Jackson K, Faulkner C, Oparka K, Somerville CR. The GRV2/RME-8 protein of Arabidopsis functions in the late endocytic pathway and is required for vacuolar membrane flow. Plant J. 2008;1:29–41.

    Article  CAS  Google Scholar 

  50. Li G, Siddiqui H, Teng Y, Lin R, Wan XY, Li J, et al. Coordinated transcriptional regulation underlying the circadian clock in Arabidopsis. Nat Cell Biol. 2011;5:616–22.

    Article  CAS  Google Scholar 

  51. Lei L, Li S, Du J, Bashline L, Gu Y. Cellulose synthase INTERACTIVE3 regulates cellulose biosynthesis in both a microtubule-dependent and microtubule-independent manner in Arabidopsis. Plant Cell. 2013;12:4912–23.

    Article  CAS  Google Scholar 

  52. Fujimori K, Ohta D. Isolation and characterization of a histidine biosynthetic gene in Arabidopsis encoding a polypeptide with two separate domains for phosphoribosyl-ATP pyrophosphohydrolase and phosphoribosyl-AMP cyclohydrolase. Plant Physiol. 1998;1:275–83.

    Article  Google Scholar 

  53. Han SK, Song JD, Noh YS, Noh B. Role of plant CBP/p300-like genes in the regulation of flowering time. Plant J. 2007;1:103–14.

    Google Scholar 

  54. Li MW, AuYeung WK, Lam HM. The GCN2 homologue in Arabidopsis thaliana interacts with uncharged tRNA and uses Arabidopsis eIF2α molecules as direct substrates. Plant Biol (Stuttg). 2013;1:13–8.

    Article  CAS  Google Scholar 

  55. Lu DP, Christopher DA. Endoplasmic reticulum stress activates the expression of a sub-group of protein disulfide isomerase genes and AtbZIP60 modulates the response in Arabidopsis thaliana. Mol Genet Genomics. 2008;3:199–210.

    Article  CAS  Google Scholar 

  56. Brooks MA, Dziembowski A, Quevillon-Cheruel S. Structure of the yeast Pml1 splicing factor and its integration into the RES complex. Nucleic Acids Res. 2009;37:129–43.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  57. Chung T, Wang D, Kim CS, Yadegari R, Larkins BA. Plant SMU-1 and SMU-2 homologues regulate pre-mRNA splicing and multiple aspects of development. Plant Physiol. 2009;3:1498–512.

    Article  CAS  Google Scholar 

  58. Hamada T, Igarashi H, Taguchi R, Fujiwara M, Fukao Y, Shimmen T, et al. The putative RNA-processing protein, THO2, is a microtubule-associated protein in tobacco. Plant Cell Physiol. 2009;4:801–11.

    Article  CAS  Google Scholar 

  59. Kang CH, Feng Y, Vikram M, Jeong IS, Lee JR, Bahk JD, et al. Arabidopsis thaliana PRP40s are RNA polymerase II C-terminal domain-associating proteins. Arch Biochem Biophys. 2009;1:30–8.

    Article  CAS  Google Scholar 

  60. Heazlewood JL, Tonti-Filippini JS, Gout AM, Day DA, Whelan J, Millar AH. Experimental analysis of the Arabidopsis mitochondrial proteome highlights signaling and regulatory components, provides assessment of targeting prediction programs, and indicates plant-specific mitochondrial proteins. Plant Cell. 2004;1:241–56.

    Article  CAS  Google Scholar 

  61. Lee CP, Eubel H, O'Toole N, Millar AH. Heterogeneity of the mitochondrial proteome for photosynthetic and non-photosynthetic Arabidopsis metabolism. Mol Cell Proteomics. 2008;7:1297–316.

    Article  CAS  PubMed  Google Scholar 

  62. Zybailov B, Rutschow H, Friso GL, Rudella A, Emanuelsson O, Sun Q, et al. Sorting Signals, N-Terminal Modifications and Abundance of the Chloroplast Proteome. PLoS One. 2008;3(4), e1994.

    Article  PubMed Central  PubMed  CAS  Google Scholar 

  63. Nikolovski N, Rubtsov D, Segura MP, Miles GP, Stevens TJ, Dunkley TP, et al. Putative glycosyltransferases and other plant Golgi apparatus proteins are revealed by LOPIT proteomics. Plant Physiol. 2012;2:1037–51.

    Article  CAS  Google Scholar 

  64. Rosgaard L, Zygadlo A, Scheller HV, Mant A, Jensen PE. Insertion of the plant photosystem I subunit G into the thylakoid membrane. FEBS J. 2005;15:4002–10.

    Article  CAS  Google Scholar 

  65. de Bianchi S, Dall'Osto L, Tognon G, Morosinotto T, Bassi R. Minor antenna proteins CP24 and CP26 affect the interactions between photosystem II subunits and the electron transport rate in grana membranes of Arabidopsis. Plant Cell. 2008;20(4):1012–28.

    Article  PubMed Central  PubMed  CAS  Google Scholar 

  66. Sawchuk MG, Donner TJ, Head P, Scarpella E. Unique and overlapping expression patterns among members of photosynthesis-associated nuclear gene families in Arabidopsis. Plant Physiol. 2008;4:1908–24.

    Article  CAS  Google Scholar 

  67. Woo HR, Kim JH, Nam HG, Lim PO. The delayed leaf senescence mutants of Arabidopsis, ore1, ore3, and ore9 are tolerant to oxidative stress. Plant Cell Physiol. 2004;7:923–32.

    Article  Google Scholar 

  68. Mikkelsen MD, Thomashow MF. A role for circadian evening elements in cold-regulated gene expression in Arabidopsis. Plant J. 2009;2:328–39.

    Article  CAS  Google Scholar 

  69. Zygadlo A, Jensen PE, Leister D, Scheller HV. Photosystem I lacking the PSI-G subunit has a higher affinity for plastocyanin and is sensitive to photodamage. Biochim Biophys Acta. 2005;2:154–63.

    Article  CAS  Google Scholar 

  70. Ferro M, Brugière S, Salvi D, Seigneurin-Berny D, Court M, Moyet L, et al. AT_CHLORO, a comprehensive chloroplast proteome database with subplastidial localization and curated information on envelope proteins. Mol Cell Proteomics. 2010;6:1063–84.

    Article  CAS  Google Scholar 

  71. Kwak JM, Mori IC, Pei ZM, Leonhardt N, Torres MA, Dangl JL, et al. NADPH oxidase AtrbohD and AtrbohF genes function in ROS-dependent ABA signaling in Arabidopsis. EMBO J. 2003;11:2623–33.

    Article  Google Scholar 

  72. Gillissen B, Bürkle L, André B, Kühn C, Rentsch D, Brandl B, et al. A new family of high-affinity transporters for adenine, cytosine, and purine derivatives in Arabidopsis. Plant Cell. 2000;2:291–300.

    Article  Google Scholar 

  73. Diener AC, Gaxiola RA, Fink GR. Arabidopsis ALF5, a multidrug efflux transporter gene family member, confers resistance to toxins. Plant Cell. 2001;7:1625–38.

    Article  Google Scholar 

  74. Mäser P, Thomine S, Schroeder JI, Ward JM, Hirschi K, Sze H, et al. Phylogenetic relationships within cation transporter families of Arabidopsis. Plant Physiol. 2001;4:1646–67.

    Article  Google Scholar 

  75. Pertry I, Václavíková K, Depuydt S, Galuszka P, Spíchal L, Temmerman W, et al. Identification of Rhodococcus fascians cytokinins and their modus operandi to reshape the plant. Proc Natl Acad Sci U S A. 2009;3:929–34.

    Article  CAS  Google Scholar 

  76. Bae MS, Cho EJ, Choi EY, Park OK. Analysis of the Arabidopsis nuclear proteome and its response to cold stress. Plant J. 2003;5:652–63.

    Article  CAS  Google Scholar 

  77. Sels J, Mathys J, De Coninck BM, Cammue BP, De Bolle MF. Plant pathogenesis-related (PR) proteins: a focus on PR peptides. Plant Physiol Biochem. 2008;11:941–50.

    Article  CAS  Google Scholar 

  78. Wormit A, Butt SM, Chairam I, McKenna JF, Nunes-Nesi A, Kjaer L, et al. Osmosensitive changes of carbohydrate metabolism in response to cellulose biosynthesis inhibition. Plant Physiol. 2012;1:105–17.

    Article  CAS  Google Scholar 

  79. Yang Y, Xu R, Ma CJ, Vlot AC, Klessig DF, Pichersky E. Inactive methyl indole-3-acetic acid ester can be hydrolyzed and activated by several esterases belonging to the AtMES esterase family of Arabidopsis. Plant Physiol. 2008;3:1034–45.

    Article  CAS  Google Scholar 

  80. Tanaka H, Ishikawa M, Kitamura S, Takahashi Y, Soyano T, Machida C, et al. The AtNACK1/HINKEL and STUD/TETRASPORE/AtNACK2 genes, which encode functionally redundant kinesins, are essential for cytokinesis in Arabidopsis. Genes Cells. 2004;12:1199–211.

    Article  CAS  Google Scholar 

  81. Schmid M, Davison TS, Henz SR, Pape UJ, Demar M, Vingron M, et al. A gene expression map of Arabidopsis thaliana development. Nat Genet. 2005;5:501–6.

    Article  CAS  Google Scholar 

  82. Borg M, Brownfield L, Khatab H, Sidorova A, Lingaya M, Twell D. The R2R3 MYB transcription factor DUO1 activates a male germline-specific regulon essential for sperm cell differentiation in Arabidopsis. Plant Cell. 2011;2:534–49.

    Article  CAS  Google Scholar 

  83. Reňák D, Dupl'áková N, Honys D. Wide-scale screening of T-DNA lines for transcription factor genes affecting male gametophyte development in Arabidopsis. Sex Plant Reprod. 2012;1:39–60.

    Article  CAS  Google Scholar 

  84. Lee JH, Yoon HJ, Terzaghi W, Martinez C, Dai M, Li J, et al. DWA1 and DWA2, two Arabidopsis DWD protein components of CUL4-based E3 ligases, act together as negative regulators in ABA signal transduction. Plant Cell. 2010;6:1716–32.

    Article  CAS  Google Scholar 

  85. Froehlich JE, Wilkerson CG, Ray WK, McAndrew RS, Osteryoung KW, Gage DA, et al. Proteomic study of the Arabidopsis thaliana chloroplastic envelope membrane utilizing alternatives to traditional two-dimensional electrophoresis. J Proteome Res. 2003;4:413–25.

    Article  CAS  Google Scholar 

  86. Carroll AJ, Heazlewood JL, Ito J, Millar AH. Analysis of the Arabidopsis cytosolic ribosome proteome provides detailed insights into its components and their post-translational modification. Mol Cell Proteomics. 2008;2:347–69.

    Google Scholar 

  87. Mähönen AP, Higuchi M, Törmäkangas K, Miyawaki K, Pischke MS, Sussman MR, et al. Cytokinins regulate a bidirectional phosphorelay network in Arabidopsis. Curr Biol. 2006;11:1116–22.

    Article  CAS  Google Scholar 

  88. Suzuki T, Miwa K, Ishikawa K, Yamada H, Aiba H, Mizuno T. The Arabidopsis sensor His-kinase, AHk4, can respond to cytokinins. Plant Cell Physiol. 2001;2:107–13.

    Article  Google Scholar 

  89. Gruis DF, Selinger DA, Curran JM, Jung R. Redundant proteolytic mechanisms process seed storage proteins in the absence of seed-type members of the vacuolar processing enzyme family of cysteine proteases. Plant Cell. 2002;11:2863–82.

    Article  CAS  Google Scholar 

  90. Yang H, Matsubayashi Y, Nakamura K, Sakagami Y. Diversity of Arabidopsis genes encoding precursors for phytosulfokine, a peptide growth factor. Plant Physiol. 2001;3:842–51.

    Article  Google Scholar 

  91. Etchells JP, Provost CM, Turner SR. Plant vascular cell division is maintained by an interaction between PXY and ethylene signalling. PLoS Genet. 2012;11, e1002997.

    Article  CAS  Google Scholar 

  92. Gidda SK, Miersch O, Levitin A, Schmidt J, Wasternack C, Varin L. Biochemical and molecular characterization of a hydroxyjasmonate sulfotransferase from Arabidopsis thaliana. J Biol Chem. 2003;20:17895–900.

    Article  CAS  Google Scholar 

  93. Libault M, Wan J, Czechowski T, Udvardi M, Stacey G. Identification of 118 Arabidopsis transcription factor and 30 ubiquitin-ligase genes responding to chitin, a plant-defense elicitor. Mol Plant Microbe Interact. 2007;8:900–11.

    Article  CAS  Google Scholar 

  94. Huibers RP, de Jong M, Dekter RW, Van den Ackerveken G. Disease-specific expression of host genes during downy mildew infection of Arabidopsis. Mol Plant Microbe Interact. 2009;9:1104–15.

    Article  CAS  Google Scholar 

  95. Ahlfors R, Brosché M, Kollist H, Kangasjärvi J. Nitric oxide modulates ozone-induced cell death, hormone biosynthesis and gene expression in Arabidopsis thaliana. Plant J. 2009;1:1–12.

    Article  Google Scholar 

  96. Liu J, Elmore JM, Lin ZJ, Coaker G. A receptor-like cytoplasmic kinase phosphorylates the host target RIN4, leading to the activation of a plant innate immune receptor. Cell Host Microbe. 2011;2:137–46.

    Article  CAS  Google Scholar 

  97. Ogawa T, Ueda Y, Yoshimura K, Shigeoka S. Comprehensive analysis of cytosolic Nudix hydrolases in Arabidopsis thaliana. J Biol Chem. 2005;26:25277–83.

    Article  CAS  Google Scholar 

  98. Brown DM, Zhang Z, Stephens E, Dupree P, Turner SR. Characterization of IRX10 and IRX10-like reveals an essential role in glucuronoxylan biosynthesis in Arabidopsis. Plant J. 2009;4:732–46.

    Article  CAS  Google Scholar 

  99. Thines B, Katsir L, Melotto M, Niu Y, Mandaokar A, Liu G, et al. JAZ repressor proteins are targets of the SCF(COI1) complex during jasmonate signalling. Nature. 2007;7154:661–5.

    Article  CAS  Google Scholar 

  100. Ghelis T, Bolbach G, Clodic G, Habricot Y, Miginiac E, Sotta B, et al. Protein tyrosine kinases and protein tyrosine phosphatases are involved in abscisic acid-dependent processes in Arabidopsis seeds and suspension cells. Plant Physiol. 2008;3:1668–80.

    Article  CAS  Google Scholar 

  101. Salzman RA, Fujita T, Zhu-Salzman K, Hasegawa PM, Bressan RA. An improved RNA isolation method for plant tissues containing high levels of phenolic compounds or carbohydrates. Plant Mol Biol Rep. 1999;17:11–7.

    Article  CAS  Google Scholar 

  102. Azevedo H, Lino-Neto T, Tavares RM. An improved method for high quality RNA isolation from needles of adult maritime pine trees. Plant Mol Biol Rep. 2003;21:333–8.

    Article  CAS  Google Scholar 

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

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  104. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P. Molecular biology of the cell. 4th ed. New York, NY: Garland Science; 2002.

    Google Scholar 

  105. Dahlberg JE, Lund E. Does protein synthesis occur in the nucleus? Curr Opin Cell Biol. 2004;16:335–8.

    Article  CAS  PubMed  Google Scholar 

  106. Go YM, Jones DP. Redox compartmentalization in eukaryotic cells. Biochim Biophys Acta. 2008;1780(11):1273–90.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  107. Lodish H, Berk A, Zipursky SL, Matsudaira P, Baltimore D, Darnell J. Molecular Cell Biology. 4th ed. New York, NY: W. H. Freeman; 2000.

    Google Scholar 

  108. Ochs RS. Biochemistry. Burlington, MA: Jones & Bartlett Learning; 2012.

    Google Scholar 

  109. Walter P, Blobel G. Translocation of proteins across the endoplasmic reticulum. II. Signal recognition protein (SRP) mediates the selective binding to microsomal membranes of in-vitro-assembled polysomes synthesizing secretory protein. J Cell Biol. 1981;91:551–6.

    Article  CAS  PubMed  Google Scholar 

  110. Walter P, Blobel G. Translocation of proteins across the endoplasmic reticulum III. Signal recognition protein (SRP) causes signal sequence-dependent and site-specific arrest of chain elongation that is released by microsomal membranes. J Cell Biol. 1981;91:557–61.

    Article  CAS  PubMed  Google Scholar 

  111. Stephens SB, Nicchitta CV. Divergent regulation of protein synthesis in the cytosol and endoplasmic reticulum compartments of mammalian cells. Mol Biol Cell. 2008;2:623–32.

    Google Scholar 

  112. Niu Y, Zhao X, Wu YS, Li MM, Wang XJ, Yang YG. N 6-methyl-adenosine (m6A) in RNA: an old modification with a novel epigenetic function. Genomics Proteomics Bioinformatics. 2013;1:8–17.

    Article  Google Scholar 

  113. Mohammad HP, Baylin SB. Linking cell signaling and the epigenetic machinery. Nat Biotechnol. 2010;10:1033–8.

    Article  CAS  Google Scholar 

  114. Perry RP, Kelley DE, Friderici K, Rottman F. The methylated constituents of L cell messenger RNA: evidence for an unusual cluster at the 5′ terminus. Cell. 1975;4:387–94.

    Article  CAS  PubMed  Google Scholar 

  115. Meyer KD, Jaffrey SR. The dynamic epitranscriptome: N 6-methyladenosine and gene expression control. Nat Rev Mol Cell Biol. 2014;5:313–26.

    Article  CAS  Google Scholar 

  116. Fu Y, Dominissini D, Rechavi G, He C. Gene expression regulation mediated through reversible m6A RNA methylation. Nat Rev Genet. 2014;5:293–306.

    Article  CAS  Google Scholar 

  117. Stoltzfus CM, Dane RW. Accumulation of spliced avian retrovirus mRNA is inhibited in S-adenosylmethionine-depleted chicken embryo fibroblasts. J Virol. 1982;42:918–31.

    CAS  PubMed Central  PubMed  Google Scholar 

  118. Camper SA, Albers RJ, Coward JK, Rottman FM. Effect of undermethylation on mRNA cytoplasmic appearance and half-life. Mol Cell Biol. 1984;4:538–43.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  119. Tuck MT, Wiehl PE, Pan T. Inhibition of 6-methyladenine formation decreases the translation efficiency of dihydrofolate reductase transcripts. Int J Biochem Cell Biol. 1999;31:837–51.

    Article  CAS  PubMed  Google Scholar 

  120. DeLeon SBT, Davidson EH. Gene regulation: Gene control network in development. Annu Rev Biophys Biomol Struct. 2007;36:191–212.

    Article  CAS  Google Scholar 

  121. Rapaport F, Khanin R, Liang Y, Pirun M, Krek A, Zumbo P, et al. Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data. Genome Biol. 2013;9:R95.

    Article  CAS  Google Scholar 

  122. Wang Y, Zhang WZ, Song LF, Zou JJ, Su Z, Wu WH. Transcriptome analyses show changes in gene expression to accompany pollen germination and tube growth in Arabidopsis. Plant Physiol. 2008;3:1201–11.

    Article  CAS  Google Scholar 

  123. McMahon MJ, Kofranek AM, Rubatzky VE. Plant Science. Upper Saddle River, NJ: Pearson Education, Inc.; 2011.

    Google Scholar 

  124. Bengough AG, Bransby MF, Hans J, McKenna SJ, Roberts TJ, Valentine TA. Root responses to soil physical conditions; growth dynamics from field to cell. J Exp Bot. 2006;2:437–47.

    Google Scholar 

  125. Feschotte C. Transposable elements and the evolution of regulatory networks. Nat Rev Genet. 2008;9:397–405.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  126. Schnable PS, Ware D, Fulton RS, Stein JC, Wei F, Pasternak S, et al. The B73 maize genome: complexity, diversity, and dynamics. Science. 2009;5956:1112–5.

    Article  CAS  Google Scholar 

  127. Bucher E, Reinders J, Mirouze M. Epigenetic control of transposon transcription and mobility in Arabidopsis. Curr Opin Plant Biol. 2012;5:503–10.

    Article  CAS  Google Scholar 

  128. Belancio VP, Hedges DJ, Deininger P. Mammalian non-LTR retrotransposons: for better or worse, in sickness and in health. Genome Res. 2008;18:343–58.

    Article  CAS  PubMed  Google Scholar 

  129. Slotkin RK, Martienssen R. Transposable elements and the epigenetic regulation of the genome. Nat Rev Genet. 2007;8:272–85.

    Article  CAS  PubMed  Google Scholar 

  130. van de Lagemaat LN, Landry JR, Mager DL, Medstrand P. Transposable elements in mammals promote regulatory variation and diversification of genes with specialized functions. Trends Genet. 2003;19:530–6.

    Article  PubMed  CAS  Google Scholar 

  131. Marino-Ramirez L, Lewis KC, Landsman D, Jordan IK. Transposable elements donate lineagespecific regulatory sequences to host genomes. Cytogenet Genome Res. 2005;110:333–41.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  132. Kang SG, Liu WH, Lu P, Jin HY, Lim HW, Shepherd J, et al. MicroRNAs of the miR-17 ∼ 92 family are critical regulators of T(FH) differentiation. Nat Immunol. 2013;8:849–57.

    Article  CAS  Google Scholar 

  133. Gambiono G, Perrone I, Gribaudo I. A rapid and effective method for RNA extraction from different tissues of grapevine and other woody plants. Phytochem Anal. 2008;19:520–5.

    Article  CAS  Google Scholar 

  134. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;4:R36.

    Article  CAS  Google Scholar 

  135. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;16:2078–9.

    Article  CAS  Google Scholar 

  136. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;7:923–30.

    Article  CAS  Google Scholar 

  137. Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative Genomics Viewer. Nat Biotechnol. 2011;1:24–6.

    Article  CAS  Google Scholar 

  138. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;5:511–5.

    Article  CAS  Google Scholar 

  139. Berardini TZ, Mundodi S, Reiser R, Huala E, Garcia-Hernandez M, Zhang P, et al. Functional annotation of the Arabidopsis genome using controlled vocabularies. Plant Physiol. 2004;2:1–11.

    Google Scholar 

Download references

Acknowledgments

This work was completed in Prof. Jian-kang Zhu’s lab. We highly appreciate Prof. Jian-kang-Zhu for his strong supports to this project including financial support for all experiments and experts in his group for analysis of a large number of datasets, and his valuable suggestions for improvement of the manuscript and correction of some scientific conceptions in the draft. We thank Xingang Wang and Huiming Zhang for providing valuable suggestions on improvement of RIP experiments in this study, and Rebecca Ann Stevenson for her valuable work on English polishing of the manuscript, to Prof. Renyi Liu and Yangqiang Li for their estimation of m6A peaks using the Moving-window program. This work was funded by the 2012 Shaanxi Province Fund for Returnees Scientists from Foreign Study (A289021201), Chinese Academy of Sciences, China Scholarship Council Project (22861057), USA NIH grants (R01GM070795 and R01GM059138) to Professor Jian-kang Zhu, and National Science Foundation of China (31101166).

Author information

Authors and Affiliations

Authors

Corresponding authors

Correspondence to Yizhen Wan or Zhaobo Lang.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

YW, KT, and XZ conceived and designed the experiment. YW grew Arabidopsis plants and collected the plant materials. YW and DZ performed the RNA extractions. YW carried out the m6A RNA IP experiment. KT and SX performed the mapping analysis. YW, KT, and SX designed the algorithm for analysis of m6A differential patterns among organs, and carried out analyses of m6A patterns, differential gene transcripts, and differential m6A methylation among organs. DZ and ZW participated in some statistical analyses. YW carried out qRT-PCR. YW and DZ analyzed the data. YW drafted the manuscript, SX, ZL, KT, ZW, and XZ provided critical comments. All authors read and approved the final manuscript.

Yizhen Wan, Kai Tang and Dayong Zhang contributed equally to this work.

Additional files

Additional file 1: Table S1.

The sequenced and mapped reads in the m6A-seq, mRNA-seq, and input RNA-seq samples. (DOC 43 kb)

Additional file 2: Figure S1.

The m6A peak and adenosine peak deduced from the HPLC-MS/MS analysis. a The relative m6A peak height (upper) and adenosine peak height (lower) in the standard sample. b The relative m6A peak height (upper) and adenosine peak height (lower) in the input sample. c The relative m6A peak height (upper) and adenosine peak height (lower) in the RIP sample. (DOC 50 kb)

Additional file 3: Table S2.

Number of transcripts in the mRNA-seq and the m6A-seq (RIP) samples and proportion of the m6A modified transcripts in the three organs of Arabidopsis. (DOC 33 kb)

Additional file 4: Table S3.

Number of m6A sites detected in the three organs of Arabidopsis. (DOC 37 kb)

Additional file 5: Table S4.

Ratio of m6A/A in the three organs of Arabidopsis. (DOC 38 kb)

Additional file 6: Table S5.

Category of the m6A modified transcripts based on the number of m6A sites per transcript. (DOC 37 kb)

Additional file 7: Table S6.

Proportion of two types of m6A distributing feature in mRNA. (DOC 30 kb)

Additional file 8: Table S7.

The transcriptome-wide normalized read depth in the 60 bins of the gene representing the overall m6A patterns in the different regions of the genes. Non-significant differences were found between three organs (P = 0.761). (DOC 108 kb)

Additional file 9: Table S8.

The primers used for qRT-PCR. (DOC 35 kb)

Additional file 10: Figure S2.

The relative abundance (RA) of m6A RNA deduced from qRT-PCR and the expected abundance (EA) of m6A RNA deduced from the m6A-seq dataset. a RA for ‘AT1G35710’, b EA for ‘AT1G35710’, c RA for ‘AT3G07610’, d EA for ‘AT3G07610’, e RA for ‘AT4G14410’, f EA for ‘AT4G14410’, g RA for ‘AT2G28490’, h EA for ‘AT2G28490’, i RA for ‘AT1G03880’, j RA for ‘AT1G03880’, k RA for ‘AT1G33700’, l EA for ‘AT1G33700’, m RA for ‘AT2G07836’, n EA for ‘AT2G07836’, o RA for ‘AT5G20960’, p EA for ‘AT5G20960’, q RA for ‘AT5G22700’, r EA for ‘AT5G22700’, s RA for ‘AT3G13400’, t EA for ‘AT3G13400’, u EA for ‘AT4G38120’, v EA for ‘AT4G38120’. (DOC 1342 kb)

Additional file 11: Figure S3.

RNA QC results of the total RNA and the RIP RNA for m6A-seq samples. a RNA quality for the total RNA sample was high with RIN over 8.5. b RNA fragmentation for the m6A-seq samples was consistent in the experiments, with an average length of 106 nt. (DOC 313 kb)

Additional file 12:

Spreadsheet files show the m 6 A peaks detected in three organs using ‘Moving-window’ program and their positions in the genome. (ZIP 547 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wan, Y., Tang, K., Zhang, D. et al. Transcriptome-wide high-throughput deep m6A-seq reveals unique differential m6A methylation patterns between three organs in Arabidopsis thaliana . Genome Biol 16, 272 (2015). https://0-doi-org.brum.beds.ac.uk/10.1186/s13059-015-0839-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s13059-015-0839-2

Keywords