Transcriptome-wide high-throughput deep m6A-seq reveals unique differential m6A methylation patterns between three organs in Arabidopsis thaliana
Genome Biology volume 16, Article number: 272 (2015)
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.
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.
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.
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 , mammals [4, 6], insects , and plants . 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 , 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 . 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.
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)  and in the rice (23 to 47 million reads for each sample) , 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  and (approximately 50 %) plant . 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).
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) , but higher than the estimated in the plant (0.7 to 1.0 m6A site per transcript) . 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 .
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 .
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.
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.
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.
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 . Both of our two groups used the similar methods to analyze the relationship between the m6A methylation extent and the gene expression level . 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.
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).
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).
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) , 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.
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).
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.
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 . 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. . However, both our study from six Arabidopsis m6A-seq datasets (Fig. 5) and a recent report from rice m6A-seq of two different tissues  did not observer this dominant m6A pattern as the authors claimed a unique m6A pattern in the plant in their paper . 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.  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 . 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 . 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 . 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 . 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 . 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 . 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 . 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).
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 , 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 . 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 . 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 . 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 . 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 . Because the fragmented transcripts were more likely to be degraded than the intact transcripts , 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 .
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 , 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.
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 . 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.
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 . 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.
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’ . The potential PCR duplicates were removed by the parameter ‘rmdup’ rooted in SAMtools . The fragment numbers for each transcript were estimated using the featureCounts with a parameter of ‘-p’ .
The peaks and distributing patterns of m6A in the Arabidopsis transcriptome were visualized using free software, Integrative Genomics Viewer (IGV2.3, Boston, MA, USA ).
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 . The consensus m6A motif sequences were figured out by Luo et al.’s protocol  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 .
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 . 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 .
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/) . 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/).
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).
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.
The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE72706.
bovine serum albumin
cetyl trimethylammonium bromide
false discovery rate
fragments per kilobase of transcript per million fragments mapped of the transcript in mRNA-seq
High-performance liquid chromatography
the Kyoto Encyclopedia of Genes and Genomes
RNA sequencing based on m6A RNA immunoprecipitation
modified fragments per kilobase of transcript per million fragments mapped of the transcript in m6A-seq
RNA sequencing from oligo-dT RNA
nicotinamide adenine dinucleotide
MFPKM in m6A-seq divided by LOG (FPKM in mRNA-seq, 2)
quality control test
quantitative real-time polymerase chain reaction
RNA integrity number
small nuclear RNA and small nucleolar RNA
the Arabidopsis Information Resource (https://www.arabidopsis.org/)
transposable element gene
transcription start site
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.
Pan T. N 6-methyl-adenosine modification in messenger and long non-coding RNA. Trends Biochem Sci. 2013;4:204–9.
Grosjean H. Fine-tuning of RNA functions by modification and editing. Berlin: Springer-Verlag; 2005.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Csepany T, Lin A, Baldick Jr CJ, Beemon K. Sequence specificity of mRNA N6-adenosine methyltransferase. J Biol Chem. 1990;265:20117–22.
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.
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.
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.
Heyndrickx KS, Vandepoele K. Systematic identification of functional plant modules through the integration of complementary data sources. Plant Physiol. 2012;159:884–901.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Hare PD, Moller SG, Huang LF, Chua NH. LAF3, a novel factor required for normal phytochrome A signaling. Plant Physiol. 2003;4:1592–604.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Mikkelsen MD, Thomashow MF. A role for circadian evening elements in cold-regulated gene expression in Arabidopsis. Plant J. 2009;2:328–39.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Ogawa T, Ueda Y, Yoshimura K, Shigeoka S. Comprehensive analysis of cytosolic Nudix hydrolases in Arabidopsis thaliana. J Biol Chem. 2005;26:25277–83.
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.
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.
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.
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.
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.
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.
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.
Dahlberg JE, Lund E. Does protein synthesis occur in the nucleus? Curr Opin Cell Biol. 2004;16:335–8.
Go YM, Jones DP. Redox compartmentalization in eukaryotic cells. Biochim Biophys Acta. 2008;1780(11):1273–90.
Lodish H, Berk A, Zipursky SL, Matsudaira P, Baltimore D, Darnell J. Molecular Cell Biology. 4th ed. New York, NY: W. H. Freeman; 2000.
Ochs RS. Biochemistry. Burlington, MA: Jones & Bartlett Learning; 2012.
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.
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.
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.
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.
Mohammad HP, Baylin SB. Linking cell signaling and the epigenetic machinery. Nat Biotechnol. 2010;10:1033–8.
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.
Meyer KD, Jaffrey SR. The dynamic epitranscriptome: N 6-methyladenosine and gene expression control. Nat Rev Mol Cell Biol. 2014;5:313–26.
Fu Y, Dominissini D, Rechavi G, He C. Gene expression regulation mediated through reversible m6A RNA methylation. Nat Rev Genet. 2014;5:293–306.
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.
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.
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.
DeLeon SBT, Davidson EH. Gene regulation: Gene control network in development. Annu Rev Biophys Biomol Struct. 2007;36:191–212.
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.
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.
McMahon MJ, Kofranek AM, Rubatzky VE. Plant Science. Upper Saddle River, NJ: Pearson Education, Inc.; 2011.
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.
Feschotte C. Transposable elements and the evolution of regulatory networks. Nat Rev Genet. 2008;9:397–405.
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.
Bucher E, Reinders J, Mirouze M. Epigenetic control of transposon transcription and mobility in Arabidopsis. Curr Opin Plant Biol. 2012;5:503–10.
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.
Slotkin RK, Martienssen R. Transposable elements and the epigenetic regulation of the genome. Nat Rev Genet. 2007;8:272–85.
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.
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.
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.
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.
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.
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.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;7:923–30.
Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative Genomics Viewer. Nat Biotechnol. 2011;1:24–6.
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.
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.
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).
The authors declare that they have no competing interests.
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.
The sequenced and mapped reads in the m6A-seq, mRNA-seq, and input RNA-seq samples. (DOC 43 kb)
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)
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)
Number of m6A sites detected in the three organs of Arabidopsis. (DOC 37 kb)
Ratio of m6A/A in the three organs of Arabidopsis. (DOC 38 kb)
Category of the m6A modified transcripts based on the number of m6A sites per transcript. (DOC 37 kb)
Proportion of two types of m6A distributing feature in mRNA. (DOC 30 kb)
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)
The primers used for qRT-PCR. (DOC 35 kb)
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)
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)
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)
About this article
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://doi.org/10.1186/s13059-015-0839-2