Shaping the landscape of N⁶-methyladenosine RNA methylation in Arabidopsis

Abstract

N⁶-methyladenosine (m⁶A) modification on messenger RNAs (mRNAs) is deposited by evolutionarily conserved methyltransferases (writers). How individual m⁶A writers sculpt the overall landscape of the m⁶A methylome and the resulting biological impact in multicellular organisms remains unknown. Here, we systematically surveyed the quantitative m⁶A methylomes at single-nucleotide resolution and their corresponding transcriptomes in Arabidopsis (Arabidopsis thaliana) bearing respective impaired m⁶A writers. The m⁶A sites associated with the five Arabidopsis writers were located mostly within 3′ untranslated regions with peaks at around 100 bp downstream of stop codons. m⁶A predominantly promoted the usage of distal poly(A) sites but had little effect on RNA splicing. Notably, impaired m⁶A writers resulted in hypomethylation and downregulation of transcripts encoding ribosomal proteins, indicating a possible correlation between m⁶A and protein translation. Besides the common effects on mRNA metabolism and biological functions uniquely exerted by different Arabidopsis m⁶A writers compared with their counterparts in human cell lines, our analyses also revealed the functional specificity of individual Arabidopsis m⁶A writers in plant development and response to stresses. Our findings thus reveal insights into the biological roles of various Arabidopsis m⁶A writers and their cognate counterparts in other multicellular m⁶A methyltransferase complexes.

Introduction

Chemical modifications of RNA are emerging as important regulatory components of gene expression. Among all these modifications, N⁶-methyladenosine (m⁶A), formed by the methylation of the N6 position of internal adenosine residues, represents the most widespread internal modification in eukaryotic messenger RNAs (mRNAs). m⁶A influences nearly every aspect of mRNA metabolism, including pre-mRNA processing (Ke et al., 2015; Haussmann et al., 2016; Zhou et al., 2019), nuclear export (Fustin et al., 2013; Roundtree et al., 2017), translation (Wang et al., 2015; Zhou et al., 2015), and decay (Wang et al., 2014; Shen et al., 2016), thereby modulating a wide range of fundamental biological processes in eukaryotes.

Deposition of m⁶A to target mRNAs requires evolutionarily conserved methyltransferases, which are collectively termed m⁶A writers. In the model plant Arabidopsis (Arabidopsis thaliana), so far two types of m⁶A writers have been suggested to install m⁶A on mRNA. The first multicomponent writer complex possibly consists of the mRNA adenosine methylase (MTA) (Zhong et al., 2008), its homolog MTB (Růžička et al., 2017), FKBP12 INTERACTING PROTEIN 37KD (FIP37) (Shen et al., 2016), VIRILIZER (VIR), and HAKAI (Růžička et al., 2017). Like their mammalian counterparts, these components install m⁶A preferentially near the stop codon and 3′ untranslated regions (UTRs), which is predominantly responsible for m⁶A levels detectable in Arabidopsis. By contrast, FIONA1 (FIO1), the only Arabidopsis ortholog of the human methyltransferase METTL16, acts separately to deposit m⁶A modifications mainly in U6 small nuclear RNA and a subset of mRNAs contributing modestly to the overall m⁶A levels (Wang et al., 2022; Xu et al., 2022).

The effects of m⁶A on mRNA stability are the best-characterized roles of m⁶A in plants to date and ultimately influence various aspects of plant developmental processes and stress responses (Anderson et al., 2018; Wei et al., 2018; Parker et al., 2020; Arribas-Hernandez et al., 2021; Shao et al., 2021; Cheng et al., 2022). Notably, knocking out any of the components except for HAKAI in the Arabidopsis writer complex results in embryonic lethality (Tzafrir et al., 2004; Vespa et al., 2004; Zhong et al., 2008; Růžička et al., 2017), which implicates an essential role of m⁶A in plants, and also reveals a technical challenge of elucidating the function of various writer components in post-embryonic development. In multicellular organisms, including Arabidopsis, some studies have been performed to understand specific m⁶A methylomes and the relevant biological responses associated with the disruption of individual m⁶A writers in a writer complex (Bodi et al., 2012; Shen et al., 2016; Parker et al., 2020; Qin et al., 2022; Zhang et al., 2022). However, as these individual m⁶A methylomes were surveyed using different methodologies and experimental materials collected under different growth conditions, it is practically unreliable to compare these methylomes with distinct detection resolutions and biological backgrounds. Thus, by far, it remains largely unknown how individual m⁶A writers in a methyltransferase complex contribute to the overall landscape of m⁶A methylomes, the global transcriptome profiles, and the associated biological processes in a given developmental context.

Among the methodologies commonly used for mapping m⁶A methylomes, both the methylated RNA immunoprecipitation sequencing (MeRIP-seq/m⁶A-seq) and m⁶A individual-nucleotide-resolution cross-linking and immunoprecipitation (miCLIP) combine immunoprecipitation using m⁶A-specific antibody with short-read cDNA-based sequencing (Dominissini et al., 2013; Linder et al., 2015). Although miCLIP can detect m⁶A at single-base resolution, it is unable to quantify differential modifications among different samples (Koh et al., 2019). In contrast, nanopore long-read direct RNA sequencing (DRS) offers a promising approach to quantitatively locate m⁶A methylation at single-nucleotide resolution and compare m⁶A methylomes across different conditions (Garalde et al., 2018; Parker et al., 2020; Pratanwanich et al., 2021). In nanopore DRS, RNA modifications cause changes in intensity levels of the electric current when a native, unamplified RNA strand traverses the protein nanopore, thus permitting precise localization of modified bases in a quantitative manner by computational methods (Liu et al., 2019; Lorenz et al., 2020).

Here we report a systematic atlas of quantitative m⁶A methylomes at single-nucleotide resolution and their corresponding transcriptomes mediated by the known individual components, including MTA, MTB, FIP37, VIR, and HAKAI, in the Arabidopsis m⁶A writer complex using nanopore DRS. Comparison of these methylomes and transcriptomes obtained from respective genetic materials grown under the same conditions permits us to uncover important characteristics of common or distinct m⁶A modification effects conferred by the different Arabidopsis m⁶A writers in the same developmental context, which provides insights into the biological roles of various Arabidopsis m⁶A writers and their cognate counterparts in other eukaryotic m⁶A methyltransferase complexes.

Results

m⁶A writers affect plant development

To systematically compare the effects of different m⁶A writers on plant development and m⁶A deposition, we created genetic materials with impaired expression or activity of MTB, VIR, and HAKAI in addition to the previously reported fip37-4 LEC1:FIP37 and AmiR-mta in which FIP37 and MTA were knocked down (Shen et al., 2016). First, we isolated two heterozygous T-DNA insertional mutants of mtb from ABRC, designated mtb-1/+ and mtb-2/+, which produced aborted mtb-1 and mtb-2 homozygous seeds (Supplemental Figure 1, A and B). We then partially complemented mtb-2 with the ABI3:MTB transgene, in which the coding sequence of MTB was driven by the embryo-specific promoter of ABA INSENSITIVE 3 (ABI3) (Parcy et al., 1994) to study the post-embryonic function of MTB. One mtb-2/+ ABI3:MTB line containing a single copy of the transgene based on the segregation ratio was selected for further characterization (Supplemental Figure 1B). Among its progenies, the homozygous mtb-2 ABI3:MTB seedlings were isolated for subsequent analyses according to their distinctive growth phenotype (Figure 1A). Second, as loss of function of VIR also led to embryonic lethality (Růžička et al., 2017), we generated its knockdown lines (AmiR-vir) by artificial microRNA (AmiR) interference (Supplemental Figure 1A) (Schwab et al., 2006). Third, we generated a mutant of HAKAI by CRISPR/Cas9-mediated gene editing using a single guide RNA targeting a region within the first exon of HAKAI. From the T3 population of the transformants, we selected a homozygous null mutant without the CRISPR/Cas9 transgene, hakai-3, which contains a 1-bp deletion of guanine (G) resulting in a premature stop codon immediately in the first exon (Supplemental Figure 1, A and C).

Next, we examined the post-embryonic phenotypes of AmiR-mta, mtb-2 ABI3:MTB, fip37-4 LEC1:FIP37, AmiR-vir, and hakai-3. All these plants, excluding hakai-3 exhibited a clear delay in leaf initiation and growth and produced leaves with irregular margins and abnormal outgrowth (Figure 1A). Strikingly, both mtb-2 ABI3:MTB and fip37-4 LEC1:FIP37 generated fasciated shoot apical meristems at the early seedling stage (Figure 1A). They subsequently developed necrotic lesions on the cotyledons and the first pair of leaves and eventually died without flowering (Figure 1A). AmiR-mta and AmiR-vir also produced large shoot apical meristems even with meristem outgrowth on lobed leaves visible at the late seedling development stage (Supplemental Figure 2A) but were eventually able to flower and set seeds. In contrast, like two previously reported hakai mutants (Růžička et al., 2017), hakai-3 did not show obvious growth defects during seedling development (Figure 1A and Supplemental Figure 2B).

In developing seedlings, MTB and FIP37 expression remained at very low levels in mtb-2 ABI3:MTB and fip37-4 LEC1:FIP37 compared to that in wild-type plants, respectively, while considerable reduction of MTA and VIR was observed in their respective AmiR lines (Figure 1B), confirming the impaired expression of m⁶A writer genes in these materials. Interestingly, we also observed reduced HAKAI expression in hakai-3 containing a point mutation (Figure 1B and Supplemental Figure 1C), implying possible nonsense-mediated decay of the mutated HAKAI transcripts harboring a premature stop codon or self-regulation of HAKAI. We then examined m⁶A levels in these genetic materials by liquid chromatography-tandem mass spectrometry (LC-MS/MS) and revealed that m⁶A was reduced by 44%–88% compared to wild-type (Figure 1C). Levels of m⁶A reduction were associated with the phenotypic severity of these plants as exemplified by a comparison between fip37-4 LEC1:FIP37 displaying the lowest m⁶A levels and the most severe growth defects and hakai-3 bearing the lowest decrease of m⁶A levels without obvious growth defects (Figure 1, A and C). Taken together, these observations demonstrate that these m⁶A writers govern m⁶A deposition associated with plant development.

Nanopore DRS of m⁶A writer mutants

To quantitatively examine the occurrence and regulation of m⁶A at single-nucleotide resolution mediated by different writers, we performed nanopore DRS on three biological replicates of 6-day-old wild-type, AmiR-mta, mtb-2 ABI3:MTB, fip37-4 LEC1:FIP37, AmiR-vir, and hakai-3 seedlings. The total number of high-quality reads (quality score ≥ 7) for each sample ranged from 1.6 to 3.4 million, and at least 96.2% of these reads were successfully mapped to the TAIR10 reference transcriptome using Minimap2 (Li, 2018) (Supplemental Figure 3A). The majority of nanopore reads had a high-quality score at around 20 with an average read length of about 840–971 nt (Supplemental Figure 3B). The longest read length detected was 12,155 nucleotides corresponding to At2g17930 which spans 35 exons (Supplemental Figure 3C). The expression of 19,035 genes was identified with good reproducibility across replicates (Supplemental Figure 4). These observations indicate the high quality of our nanopore data set used for subsequent analyses. Moreover, nanopore reads confirmed reduced transcript levels of MTA and VIR, truncated transcripts of MTB and FIP37, and the 1-bp deletion of G in HAKAI in their respective genetic materials (Figure 1, D–H), further corroborating the reliability of our nanopore DRS reads.

Precise and quantitative detection of hypomethylated m⁶A sites

Information on m⁶A modification can be extracted from the nanopore reads through several computational algorithms (Huang et al., 2019; Lorenz et al., 2020; Parker et al., 2020; Pratanwanich et al., 2021), among which xPore not only identifies m⁶A sites at the single base resolution, but also quantifies the fraction of modified reads known as the modification rate, thus providing an unprecedented detail of epitranscriptome (Pratanwanich et al., 2021). We thus analyzed our nanopore reads using xPore to profile m⁶A sites at single-nucleotide resolution and to compute the differential modification rates (DMRs) representing the differences in modification rates between the respective writer genetic materials and wild-type seedlings. Based on stringent criteria (P < 0.001; one-directional signal shifts; detectable in all biological replicates), we identified 6,899, 6,022, 3,650, 6,594, and 3,015 high-confidence hypomethylated m⁶A sites in AmiR-mta, mtb-2 ABI3:MTB, fip37-4 LEC1:FIP37, AmiR-vir, and hakai-3, respectively (Figure 2A and Supplemental Figure 5, Supplemental Data Set 1). While most of the DMRs in hakai-3 were distributed between −0.50 and 0, the patterns in other genotypes were similar to the majority of DMRs being distributed between −1.00 and −0.25. The lowest number of hypomethylated sites and overall lower DMRs observed in hakai-3 were consistent with the mild reduction of m⁶A levels as measured by LC-MS/MS (Figure 1C). In total, 15,551 unique hypomethylated m⁶A sites were identified, corresponding to 4,315 protein-encoding genes (Supplemental Data Set 1). We randomly selected a few candidates and validated the sites using m⁶A–IP–qPCR (Supplemental Figure 6), indicating the authenticity of hypomethylated sites revealed in our study. Furthermore, a comparison of our data with the miCLIP data set (Parker et al., 2020) revealed that 44.3%–56.5% of our hypomethylated sites were within 5 nt of miCLIP peaks (Figure 2, B and C), while 57.5%–77.4% of our hypomethylated m⁶A sites were located within m⁶A peaks identified by m⁶A-seq data (Shen et al., 2016; Duan et al., 2017) (Supplemental Figure 7A). Non-overlapping sites revealed in our nanopore DRS versus the previous studies could be due to distinct approaches used for identifying m⁶A sites and/or plant materials grown at different conditions.

We next examined the distribution pattern of the hypomethylated sites along the transcripts and found that m⁶A was markedly enriched near the stop codon and 3′ UTR for all genetic materials examined (Figure 2D), which is in accordance with previous observations (Shen et al., 2016; Parker et al., 2020). A closer inspection revealed the enrichment of the hypomethylated sites at ∼100 bp downstream of the stop codon in all data sets (Figure 2E). Notably, the highest percentage of the hypomethylated sites were in the region 100 bp upstream or downstream of the stop codon in all genetic materials except for hakai-3, in which more than one-third of the hypomethylated sites were in the coding region (Figure 2F). Among hypomethylated transcripts found in this study, more than half of the genes contained only one hypomethylated m⁶A site in hakai-3, and this percentage was higher than those in other genetic backgrounds (Supplemental Figure 7B). In contrast, AmiR-mta had the highest percentage of genes (63%) containing more than one hypomethylated site (Figure 2G and Supplemental Figure 7B), reflecting the role of MTA as the major catalytically active subunit of the m⁶A methyltransferase complex.

Comparative m⁶A methylomes mediated by different m⁶A writers

To evaluate the sequence features of hypomethylated m⁶A sites associated with individual m⁶A writers, we performed motif search of a ± 20 bp region surrounding the hypomethylated site using MEME (Bailey and Elkan, 1994). The most frequently detected motif was RRACH (R = G/A, H = A/U/C) for all genotypes except for hakai where a more specific motif, ARACH, was identified (Figure 3A), implying that HAKAI is likely involved in m⁶A deposition on a smaller subset of target transcripts. Meanwhile, UGGAY (Y = U/C) and GGAUU were other significant motifs for a smaller number of consensus sequences (Supplemental Figure 7C).

Notably, we found that only 8.5%–19.6% (relative ratios of individual writers) of the hypomethylated sites were common to all five writers (Figure 3B and Supplemental Figure 8A), while on a transcript basis, the ratios of commonly modified transcripts were higher, ranging from 28.0% to 51.3% (Figure 3C). The ratios of common sites or transcripts were increased when we compared the profiles among a smaller number of writers. These observations indicate that individual m⁶A writers in the methyltransferase complex may modify different targets by acting in various combinations (Figure 3, B and C).

In contrast, each impaired m⁶A writer also contained some unique hypomethylated sites and transcripts, indicating their respective indispensable roles in the methyltransferase complex. For example, we found that 3,024 and 2,151 hypomethylated m⁶A sites were unique to AmiR-mta and mtb-2 ABI3:MTB, respectively (Figure 3B and Supplemental Figure 8, B and C). Since both MTA and MTB contain the core catalytic motif “DPPW” (Supplemental Figure 9A) within the MT-A70 methyltransferase domain that is essential for m⁶A methylation (Shih et al., 2019), MTA and MTB may function both cooperatively and independently to deposit m⁶A to their target transcripts. To this end, we examined their in vitro m⁶A methyltransferase activity and found that MTB exhibited an even stronger m⁶A methyltransferase activity than MTA (Supplemental Figure 9B). Gene ontology (GO) enrichment analysis of hypomethylated transcripts unique to impaired MTA and MTB also revealed some distinct GO terms between these two writers (Supplemental Figure 9C), indicating that they independently contribute to m⁶A modification associated with different aspects of plant growth and stress responses.

Functional enrichment analysis of hypomethylated transcripts

To understand the biological context of m⁶A modification, functional enrichment analysis was carried out on the hypomethylated genes associated with each impaired writer using clusterProfiler R package (Wu et al., 2021). Top GO terms that were common to all five writers included ribosome biogenesis, response to cytokinin, cytoplasmic translation, and responses to cadmium ion and salt stress (Figure 3D), demonstrating that these writers methylate target transcripts with similar functions in some key cellular processes and responses to stresses. As the overall reduction in the methylation rate was lower in hakai-3 than in the other genetic backgrounds (Figure 2A and Supplemental Data Set 1), the numbers of genes in these GO terms in hakai-3 were usually the lowest among all genotypes examined.

Notably, when we interrogated the list of hypomethylated genes for the enrichment of cellular compartments, those genes associated with chloroplast/plastid appeared as the most significantly enriched GO terms (Supplemental Figure 10A). Loss of methylation in these genes could contribute to the severe developmental defects observed in plants with impaired m⁶A writers except hakai-3 (Figure 1A). These genes included those encoding essential proteins that constitute a functional chloroplast. For example, CELL GROWTH DEFECT FACTOR 1 (AT5G23040; Supplemental Data Set 1) is required for chlorophyll synthesis (Reinbothe et al., 2015), which is a functional GO term appearing in all genotypes except hakai-3 (Figure 3D). As depletion of this gene results in embryo arrest at the globular stage (Kawai-Yamada et al., 2014), its hypomethylated status might be associated with developmental arrest displayed by the seedlings with impaired m⁶A writers except hakai-3. Another example is CRUMPLED LEAF (CRL, AT5G51020; Supplemental Data Set 1) which encodes a chloroplast protein and its mutation leads to severe developmental defects with the constitutive expression of stress-related genes (Hudik et al., 2014), which resembles the transcriptomes of plants with impaired m⁶A writers (Figure 4). Indeed, a comparison of their transcriptomes showed a moderate positive correlation [Spearman correlation coefficient (ρ) = 0.62–0.75] between crl and the plants with impaired m⁶A writers except for hakai-3 (Supplemental Figure 10B). While these observations provide circumstantial evidence linking chloroplast function with m⁶A modification, further experiments are required to test whether loss of m⁶A at chloroplast-related genes could be associated with chloroplast malfunction and the resulting severe developmental defects exhibited by the seedlings with impaired m⁶A writers except hakai-3.

Besides the common GO terms, we also found uniquely enriched GO terms for hypomethylated transcripts in plants with individual impaired writers (Supplemental Figure 11), indicating that individual m⁶A writers also specifically deposit m⁶A on different target transcripts involved in distinct biological processes in Arabidopsis.

Comparative transcriptomes affected by individual m⁶A writers

To understand the effects of individual m⁶A writers on the transcriptome, we performed differential gene expression analysis of the nanopore reads between impaired writers and wild-type plants using DESeq2 (Love et al., 2014). Hierarchical clustering of differentially expressed genes revealed that the transcriptomes of wild-type and hakai-3 plants were grouped together, while those of the other four genotypes contained two separate clusters with the mtb-2 ABI3:MTB transcriptome as an independent cluster (Figure 4A and Supplemental Data Set 2). Consistently, the mutation in HAKAI had the least disturbance on the transcriptome with the lowest number of differentially expressed genes among the genotypes with impaired m⁶A writers (Figure 4B). There were 211 and 264 commonly up- and downregulated genes, respectively, in all genotypes except hakai-3, while these numbers considerably dropped to only 83 and 87, respectively, if hakai-3 was included in the comparison (Figure 4C), substantiating that HAKAI influences a smaller group of genes than other writers.

GO analysis on differentially upregulated or downregulated transcripts revealed that responses to various abiotic and biotic stimuli, including oxidative stress, fungus, salt stress, karrikin and cadmium ion, were among the top five functional categories common to all genotypes (Figure 4D), implying that m⁶A greatly affects plant responses to external stimuli. Interestingly, in contrast to a smaller number of genes associated with two out of the five top GO terms in the downregulated list, all transcripts associated with the top five GO terms in the upregulated list were relevant to stress responses (Figure 4D), indicating that m⁶A modification is more prone to suppress the upregulation of stress-related transcripts, which likely occurs both directly and indirectly as some of these transcripts are hypomethylated (Supplemental Data Set 3).

Since response to salt stress was the top GO term associated with both the hypomethylated and differentially upregulated transcripts (Figures 3D and 4D), we tested the effect of defective m⁶A on plant response to salt stress. Seedlings with impaired m⁶A writers exhibited greater susceptibility to salt treatment than wild-type seedlings (Supplemental Figure 12), as shown in a previous study (Hu et al., 2021), implicating m⁶A in regulating salt stress responses.

We subsequently examined the top five enriched GO terms that were unique in each of the differentially expressed data sets (Supplemental Figure 13). For upregulated transcripts, genes associated with GO terms related to primary metabolism and photosynthesis were uniquely prevalent in AmiR-mta and fip37-4 LEC1:FIP37 transcriptomes, respectively. Growth-related GO terms were associated distinctly with upregulated transcripts in hakai-3, while some GO terms associated with stress and defense were only found in AmiR-vir. For downregulated transcripts, GO terms related to carbohydrate derivative biosynthesis, and several metabolic processes were prevalent in mtb-2 ABI3:MTB, while responses to hypoxia and oxygen were notable GO terms unique in hakai-3. These observations substantiate that individual m⁶A writers affect distinct aspects of plant growth and stress response despite some of their common functions.

m⁶A hypomethylation mainly decreases the stability of transcripts

To examine the potential association between m⁶A modification and mRNA abundance, we compared differentially expressed genes with hypomethylated genes. The transcripts containing hypomethylated m⁶A sites were predominantly downregulated in all plants with impaired m⁶A writers (Figure 5A), indicating that m⁶A modification mainly stabilizes mRNA in Arabidopsis. This is also consistent with a previous report showing that m⁶A primarily protects Arabidopsis mRNA from endonucleolytic cleavage (Anderson et al., 2018). Subsequent examination of the enriched GO terms for hypomethylated transcripts (DMR < −0.3) that were differentially regulated disclosed a greater overlap of functional categories among the different genotypes for hypomethylated transcripts that were downregulated than those upregulated (Supplemental Figure 14).

One stark example of the top functional category associated with hypomethylated transcripts that were downregulated includes those encoding the highly conserved ribosomal proteins (RPs) that make up ribosomal subunits involved in the cellular process of translation (Figure 5B). There are 79 RPs in Arabidopsis, with each encoded by several paralogs, thereby giving rise to ribosomes that are heterogeneous in their composition (Martinez-Seidel et al., 2020). Among the 242 genes encoding different RPs (Martinez-Seidel et al., 2020), transcripts corresponding to 207 genes (85.5%) that cover all RP families in Arabidopsis were hypomethylated in at least one type of the plants with impaired m⁶A writers (Supplemental Data Set 1). In particular, most of these hypomethylated RP transcripts (188 out of 207) were downregulated in at least one genotype (Supplemental Data Set 1), indicating that m⁶A methylation mainly stabilizes RP paralogs.

To verify the effect of m⁶A modification on mRNA decay of RP genes, we measured the levels of two randomly selected downregulated RP transcripts in 6-day-old seedlings with impaired m⁶A writers treated with the transcription inhibitor, actinomycin D. Both transcripts displayed obvious decreases in their cumulative RNA stability in these seedlings versus wild-type seedlings after actinomycin D treatment, whereas the stability of non-m⁶A targeted transcripts, such as β-tubulin (TUB8), remained almost unchanged (Figure 5C). This result confirms that m⁶A methylation decreases the decay of these two RP transcripts in seedlings. It will be interesting to further investigate whether the reduction in mRNA abundance for many RP paralogs mediated by impaired m⁶A methylation could influence their protein levels, thus affecting their participation in functional ribosomes and the subsequent translation processes.

m⁶A hypomethylation disrupts mRNA 3′-end processing

As m⁶A was evidently enriched near the stop codon and 3′ UTR (Figure 2D), we subsequently examined the impact of defective m⁶A deposition on two major steps of pre-mRNA processing, splicing and 3′-end polyadenylation. Only a small number of differential alternative splicing events, such as intron retention and exon skipping, were detected in plants with impaired m⁶A writers (Supplemental Figure 15A; Supplemental Data Set 4). In addition, the numbers of hypomethylated and non-hypomethylated transcripts that undergo different types of alternative splicing were mostly similar in the plants with defective m⁶A deposition (Supplemental Figure 15B). These observations suggest that m⁶A has a little direct effect on RNA splicing.

Polyadenylation of pre-mRNA is a key processing step in mRNA maturation and more than 70% of Arabidopsis genes have at least two poly(A) sites, leading to the expression of alternative polyadenylation (APA) isoforms (Wu et al., 2011). Since most APA sites were found within the 3′ UTR (Wu et al., 2011) where the majority of m⁶A sites are located (Figure 2D), we examined the effect of defective m⁶A deposition on APA and revealed clear defects in RNA 3′-end formation in 2159, 2168, 1414, 2545, and 595 genes, respectively, in the plants with different impaired m⁶A writers (Figure 6A and Supplemental Data Set 5). Strikingly, most of these transcripts (91.3%–94.6%) displayed a shift to the usage of proximal poly(A) sites (pPAS) (Figure 6A), which is also observed in the hypomorphic vir mutants (Parker et al., 2020). Nearly half of these genes (40.9%–51.6%; hypergeometric P value ≤ 1.49E−67) also contained hypomethylated m⁶A sites (Figure 6, B and C), implying that m⁶A could affect the usage of alternative poly(A) sites in plants. One example of hypomethylated genes with an increased pPAS usage in plants with impaired m⁶A writers except hakai-3 is AT1G77490 (Figure 6D), which encodes a thylakoid ascorbate peroxidase that plays a major role in regulating H₂O₂ level and thereby regulates H₂O₂-triggered retrograde signaling from chloroplasts to nucleus in response to stress (Maruta et al., 2012).

As m⁶A has been reported to protect the transcriptome integrity in plants (Pontier et al., 2019), we examined possible readthrough at poly(A) sites but found that very few transcripts were significantly affected in plants with impaired m⁶A writers (P value < 0.05; Supplemental Data Set 6). Overall, these results indicate a predominant effect of m⁶A modification on promoting the usage of distal poly(A) sites.

Comparison of m⁶A modification in plants and human cell lines

To gain a better understanding on the regulatory roles of m⁶A, we compared m⁶A features uncovered from this study with published data on the orthologs of plant writers in human cell lines, including METTL3 (ortholog of MTA), METTL14 (ortholog of MTB), Wilm's tumor 1-associating protein (WTAP) (ortholog of FIP37), VIRMA (ortholog of VIR), and HAKAI. As only one nanopore DRS data set on human embryonic kidney cell line with METTL3 knockout (HEK293T METTL3 KO versus HEK293T) was available (Pratanwanich et al., 2021), we also included published meRIP data derived from defective HeLa cell lines created either by CRISPR/CAS9 knockout of METTL3, METTL14 and VIRMA, or by siRNA knockdown of all five writers (Wang et al., 2014; Yue et al., 2018).

Although the m⁶A motif, RRACH, is conserved in animals and plants (Shao et al., 2021), the distribution of m⁶A enrichment along mRNA transcripts differed between Arabidopsis and human cell lines (Figure 2, E and 7A). m⁶A was highly enriched near stop codons with almost comparable distribution in surrounding coding regions and 3′ UTR in human HeLa cell lines (Wang et al., 2014; Yue et al., 2018) (Figure 7A), while m⁶A was mostly enriched in 3′ UTR with peaks downstream of stop codons in Arabidopsis (Figure 2E). A small number of m⁶A sites in coding regions of Arabidopsis mRNAs could be associated with weak effects of m⁶A modification on alternative splicing in Arabidopsis (Supplemental Figure 15A). In contrast, m⁶A effects on alternative splicing in animals are still under debate because of contradictory findings so far reported (Dominissini et al., 2012; Geula et al., 2015; Ke et al., 2017).

Besides the known m⁶A effect on accelerating mRNA decay in HeLa cells (Wang et al., 2014; Yue et al., 2018), we further assessed the influence of defective m⁶A deposition on the APA of mRNAs. Although there was a pronounced 3′ UTR lengthening of transcripts associated with hypomethylated peaks resulting from knockout of METTL3 or VIRMA (Yue et al., 2018) (Figure 7B), the numbers of affected transcripts were much lower than those found in their counterparts in Arabidopsis (Figure 6A). In addition, such an effect was either undetectable in siWTAP (Wang et al., 2014) or less obvious in the knockout of METTL14 and depletion of HAKAI (Yue et al., 2018) (Figure 7B). Surprisingly, analysis of APA in HEK293T cell line defective in METTL3 using nanopore DRS did not reveal m⁶A effect on poly(A) site choice (Pratanwanich et al., 2021) (Figure 7C). While these paradoxical findings indicate that different human cell lines and sequencing approaches may affect the detection of the effects of various m⁶A writers on APA, the overall results so far demonstrate a less prevalent effect of m⁶A on APA in human cell lines than in plants.

We also performed functional enrichment analysis on the hypomethylated genes associated with each impaired writer in human HeLa cell lines (Wang et al., 2014; Yue et al., 2018). Top GO terms common to all five writers included those associated with histone modification, chromatin modification, mitochondrial protein and gene expression, and protein modification (Figure 7D and Supplemental Figure 16). These categories are different from the top enriched functional GO terms in Arabidopsis (Figure 3D) implying different biological roles of m⁶A modification in plants and human cell lines.

Discussion

The m⁶A writer complex comprising multiple components catalyzes m⁶A deposition to target mRNAs in eukaryotes. Although some studies have been performed to understand the roles of individual writer components in determining m⁶A deposition on their targets, variations in biological materials and sequencing approaches have so far prevented a systematic survey on the common or distinct functions of various writer components in sculpting the landscape of m⁶A methylome in a given biological context. In this study, we have carried out nanopore DRS to compare the methylomes and transcriptomes generated from respective Arabidopsis genetic materials where the five writer components are impaired. As these plants were grown and harvested under the same conditions, the quantitative m⁶A methylomes at single-nucleotide resolution and their corresponding transcriptomes permit an unambiguous elucidation of different and complementary roles of individual writer components in determining m⁶A methylomes and the associated biological effects in Arabidopsis (Figure 8).

Among the known components of the m⁶A writer complex in Arabidopsis, HAKAI has the lowest effects on m⁶A modification and plant development. First, unlike other writers, knocking out HAKAI does not cause obvious phenotypes during seed development and post-embryonic development (Tzafrir et al., 2004; Vespa et al., 2004; Zhong et al., 2008; Růžička et al., 2017) (Figure 1A and Supplemental Figure 2B). Second, among the plants bearing the impaired writers, disruption of HAKAI in hakai-3 results in the lowest reduction of m⁶A levels (Figure 1C) associated with the lowest number of hypomethylated sites (Figure 2A). Third, the transcriptome of wild-type plants more closely resembles that of hakai-3 than other impaired writers (Figure 4A), demonstrating that HAKAI influences a smaller number of transcripts than other writers. These pieces of evidence suggest that compared to the other four writers, HAKAI mediates m⁶A deposition to the least number of transcripts, conferring a minimal effect on the landscape of m⁶A methylome and transcriptome as well as the resulting plant phenotypes.

Despite various strengths exhibited by the five writers on m⁶A deposition, their effects on mRNA exhibit some common features uniquely in Arabidopsis compared with those in human cell lines (Figure 8). Most of the m⁶A sites associated with the five Arabidopsis writers are located within 3′ UTR with peaks at 100 bp downstream of stop codons, which contrasts with a comparable distribution of m⁶A sites in coding regions and 3′ UTR surrounding the peaks near stop codons in human cell lines. The differential m⁶A localization pattern is coupled with distinct effects on mRNA metabolism between Arabidopsis and human cell lines. m⁶A deposited by the Arabidopsis writers mainly stabilizes mRNA and promotes the usage of distal poly(A) sites, whereas m⁶A accelerates mRNA decay but has a less prominent effect on determining APA sites in human cell lines, and if any, m⁶A promotes the usage of pPAS (Wang et al., 2014; Yue et al., 2018; Pratanwanich et al., 2021). In addition, Arabidopsis writers deposit m⁶A on the transcripts with similar functions in some key biological processes, such as ribosome biogenesis and response to stresses. These functions are different from the top enriched ones associated with m⁶A writers in human cell lines. Together, these common features displayed by Arabidopsis m⁶A writers imply that different components of the m⁶A methyltransferase complex may have evolved to act in concert to install m⁶A on some key transcripts, which affects RNA metabolism and the resulting biological functions distinctively in multicellular plants.

Besides the common features, our analyses of nanopore DRS have also revealed the functional specificity of individual Arabidopsis m⁶A writers. The plants bearing individual impaired m⁶A writers all contain a group of unique hypomethylated and differentially expressed transcripts that are absent in other genetic backgrounds. These transcripts are enriched for different biological processes (Figure 8), indicating that the individual m⁶A writers specifically affect certain aspects of plant development and response to stresses by directly catalyzing m⁶A modification of the relevant transcripts and/or indirectly influencing gene expression as a result of m⁶A modification.

In addition to enhancing mRNA stability, it is possible that m⁶A modification influences other RNA metabolic processes, such as nuclear transport and translation, especially for many transcripts that are hypomethylated without any changes in their steady-state levels (Figure 3, D and 5A). This possibility is supported by two observations. First, there is a predominant effect of m⁶A modification on promoting the usage of distal poly(A) sites in Arabidopsis. This alters 3′ UTR length that could influence post-transcriptional gene regulation, including nuclear export and protein translation (Tian and Manley, 2017; Turner et al., 2018). Second, a high proportion of hypomethylated transcripts that are differentially expressed in the plants with impaired writers are those encoding RPs (Figure 5B), which are essential subunits in various functional ribosomes for translation. This observation implies that alteration in protein translation could be associated with m⁶A effects on many transcripts because of defective ribosomes. This potential role of m⁶A in ribosomal code is likely specific to plants as there is a lack of mRNAs encoding RPs among the m⁶A-containing mRNAs in yeast and mammalian cell cultures (Schwartz et al., 2014).

In conclusion, our study has provided unprecedented insights into the common and unique features of the methylomes and the corresponding transcriptomes associated with the five m⁶A writers in the Arabidopsis m⁶A methyltransferase complex. In addition to some common effects on mRNA metabolism uniquely exerted by these Arabidopsis writers versus their human counterparts, our findings also reveal the functional specificity of individual Arabidopsis writers in different biological processes. Further elucidation of the mechanisms underlying the target specificity of m⁶A deposition for various m⁶A writers in multicellular organisms, including Arabidopsis, requires in-depth investigation of the function and dynamics of RNA modification-related proteins that interact with an individual or grouped m⁶A writers in given biological contexts.

Materials and methods

Plant materials and growth conditions

Arabidopsis (Arabidopsis thaliana) plants were grown on soil or Murashige and Skoog (MS) plates under a 16 h light/8 h dark-light cycle at 23°C ± 2°C. Plants used were in the Columbia (Col) background. Seeds of mtb-1 (CS16197) and mtb-2 (CS850592) were obtained from the Arabidopsis Biological Resource Center (http://www.arabidopsis.org/).

Vector construction and plant transformation

To construct ABI3:MTB, the coding sequence of MTB was inserted downstream of the 5.0 kb ABI3 promoter in the modified pENTR vector (Shen et al., 2016). To construct AmiR-vir, the AmiR was designed by WMD3 Web MicroRNA Designer (http://wmd3.weigelworld.org/). Based on the AmiR sequence, a set of four primers were generated and used for PCR amplification according to the published protocol (Schwab et al., 2006). The resulting PCR fragment was digested with EcoRI and BamHI and cloned into a modified pENTR vector containing the 2 × 35S promoter. To construct CAS9-HAKAI, the CAS9-sgRNA cassette was amplified and cloned into the modified pENTR vector. The resulting vector was digested with BbsI and ligated with the synthesized sgRNA oligos. Agrobacterium tumefaciens-mediated transformation of Arabidopsis was performed with the floral dipping method. To generate mtb-2/+ ABI3:MTB, the ABI3:MTB construct was transformed into mtb-2/+ plants and selected by kanamycin on MS plates. T1 transformants harboring both ABI3:MTB and mtb-2/+ were maintained as independent lines, from which one representative line that contained only a single ABI3:MTB transgene based on the segregation ratio was selected for further investigation. AmiR-vir and CAS9-HAKAI were transformed into wild-type Col plants and selected by Basta on soil.

Expression analysis

Total RNA was isolated with the RNeasy Plus Mini Kit (Qiagen) and reversed transcribed with the M-MLV Reverse Transcriptase (Promega). Reverse transcription-quantitative PCR (RT–qPCR) was performed with three biological replicates on the 7900HT Fast Real-Time PCR systems (Applied Biosystems) using PowerUp SYBR Green Master Mix (Applied Biosystems). The differences between the cycle threshold (Ct) of target genes and the Ct of TUBULIN 2 (TUB2) (ΔCt = Ct_{target gene—}Ct_TUB2) served as an internal control for determining the normalized expression of target genes. Primers used for gene expression analysis are listed in Supplemental Data Set 7.

Measurement of m⁶A level via LC-MS/MS

Total RNA of 6-day-old seedlings of wild-type, AmiR-mta, mtb-2 ABI3:MTB, fip37-4 LEC1:FIP37, AmiR-vir, and hakai-3 was extracted using Trizol reagent (Invitrogen). mRNA was subsequently isolated using Dynabeads Oligo(dT)₂₅ (Invitrogen)and was digested to single ribonucleosides following a published protocol (Thüring et al., 2017). Multiple reaction monitoring quantification was performed on a NanoLC 425 (Eksigent) coupled to a QTrap 6500 mass spectrometer (SCIEX) to detect A and m⁶A with mass transitions at 268.0–136.0 and 282.0–150.1, respectively. Solvent A was 2% (v/v) acetonitrile with 0.1% (v/v) formic acid, while solvent B consisted of 98% (v/v) acetonitrile with 0.1% (v/v) formic acid. The samples were resolved on a Kinetex F5 analytical column (2.6 μm, 0.3 × 100 mm, Phenomenex) using an elution gradient of 0%–30% solvent B over 5 min at a flow rate of 8 µl/min.

Nanopore DRS

Total RNA was extracted from three biological replicates of 6-day-old seedlings of wild-type, AmiR-mta, mtb-2 ABI3:MTB, fip37-4 LEC1:FIP37, AmiR-vir, and hakai-3 using Trizol reagent (Invitrogen). mRNA was then isolated using the Dynabeads mRNA purification kit (Invitrogen). The quality and integrity of mRNA were assessed with an Agilent Bioanalyzer system. Each library was prepared with 750 ng of mRNA using the nanopore direct RNA sequencing kit (SQK-RNA002, Oxford Nanopore Technologies). The prepared libraries were loaded onto FLO-MIN106 flow cells and sequenced with the GridION sequencer. The raw fast5 data were basecalled by Guppy 4.2.3 with the high accuracy mode to generate FASTQ files.

Identification of m⁶A sites

The FASTQ reads were mapped to the reference transcriptome of TAIR10 by Minimap2 2.20 (Li, 2018). Alignment was subsequently converted to a BAM file by Samtools (Li et al., 2009) and Nanopolish Eventalign v0.13 (Loman et al., 2015) was subsequently used for signal segmentation. Both aligned reads and events with the reference annotation (Ensemble plant release 50) were processed with xPore (Pratanwanich et al., 2021) to determine differential m⁶A modification sites.

Differential expression analysis with nanopore reads

Read counts for wild-type plants and plants with impaired m⁶A writers were performed with Subread v2.0.3 and featureCounts v2.0.2 in the long-read mode (Liao et al., 2014, 2019). Differentially expressed genes were then identified using DESeq2 v1.32 (Love et al., 2014). Heatmaps for transcriptomics data were created via the pheatmap package using log-transformed expression value (http://rpackages.ianhowson.com/cran/pheatmap/). Hierarchical clustering of genes and samples in heatmaps was done using the complete-linkage clustering method and Euclidian distances.

Identification of alternative 3′ end positions

The poly (A) and mRNA regions were segmented using the poly (A) script of Nanopolish v0.13 (Loman et al., 2015) that identifies authentic 3′end for each read. To identify the differential usage of each position between wild-type plants and plants with impaired m⁶A writers, we performed the Kolmogorov–Smirnov test followed by multiple testing correction with the Benjamini–Hochberg method. The differential 3′ end usage results were filtered for FDR < 0.05 and normalized to the total number of reads. The difference between normalized numbers of plants with impaired m⁶A writers and wild-type plants for each site was used to obtain the minimum and maximum change in 3′UTR length, which represents the most reduced and increased 3′ end usage sites, respectively. The difference between the most reduced and increased 3′ end usage sites was then obtained to estimate the direction and distance of each change.

Alternative splicing analysis

Full-length alternative isoform analysis of RNA (FLAIR v1.5) (Tang et al., 2020) was performed on the nanopore reads to detect the differential alterative splicing events between wild-type plants and plants with impaired m⁶A writers.

m⁶A–IP–qPCR

m⁶A–IP–qPCR was performed using 6-day-old seedlings in different genetic backgrounds as previously described (Dominissini et al., 2013; Shen et al., 2016). m⁶A-containing candidate genes and an internal control (TUB2) mRNA were detected by qPCR on the 7900HT Fast Real-Time PCR systems (Applied Biosystems) with PowerUp SYBR Green Master Mix (Applied Biosystems). The primers used for m⁶A–IP–qPCR analysis are listed in Supplemental Data Set S7.

Analysis of human data

MeRIP data (GSE102493 and GSE46705) (Wang et al., 2014; Yue et al., 2018) were downloaded from GEO (Supplemental Data Set 8, A–E). The raw SRA files are converted to FASTQ files using SRA Toolkit. The data was processed and quality-controlled by nf-core/rnaseq pipeline (https://github.com/nf-core/rnaseq). The reads were aligned to human hg19 genome (version GRCh37.87) by HISAT2 (Kim et al., 2019). The differential methylation peaks between wild-type and the respective cell lines with defective writers were detected using exomePeak (Meng et al., 2013) and the peak distribution was generated by Guitar (Cui et al., 2016). Hypomethylated peaks with log₂(fold change) < 0 and FDR < 0.05 were selected for further analysis. DaPars was used for the de novo identification of APA (Xia et al., 2014) with FDR < 0.05 and at least two-fold increase or decrease of PDUIs between wild-type and defective cell lines (Supplemental Data Set 8, F–J). For human nanopore DRS data (PRJEB40872) (Pratanwanich et al., 2021), the FASTQ files were downloaded from ENA. Nanopore reads were aligned to the hg19 genome (version GRCh37.87) using minimap2 (Li, 2018). The APA sites were identified as described for our Arabidopsis nanopore data, while the m⁶A peaks were retrieved from published data sets (Pratanwanich et al., 2021).

Accession numbers

Sequence data from this article can be found in the Arabidopsis Information Resource (https://www.arabidopsis.org) under the following accession numbers: MTA (AT4G10760), MTB (AT4G09980), FIP37 (AT3G54170), VIR (AT3G05680), and HAKAI (AT5G01160).

Data availability

All data from this study are available in the main text and Supplemental Information. The nanopore direct RNA sequencing data described in this study have been deposited in NCBI Sequence Read Archive (SRA) database with the accession numbers, PRJNA815899 (plants with impaired m⁶A writers) and PRJNA749003 (wild-type plants).

Supplemental data

The following materials are available in the online version of this article.

Supplemental Figure S1. Generation of various genetic materials where the expression or activity of MTB, VIR, and HAKAI is impaired.

Supplemental Figure S2. Phenotypes of 1-month-old AmiR-mta, AmiR-vir, and hakai-3.

Supplemental Figure S3. Analysis of nanopore DRS reads.

Supplemental Figure S4. Heatmap showing the correlation of gene expression across three biological replicates of all six samples in different genetic backgrounds.

Supplemental Figure S5. Distribution of differential methylation rate (DMR) and P value for hypomethylated sites identified from the comparison between the respective genotypes with impaired m⁶A writers and wild-type plants.

Supplemental Figure S6. Verification of m⁶A sites identified through nanopore DRS analysis for five randomly selected candidates.

Supplemental Figure S7. Characteristics of hypomethylated sites identified in plants with impaired m⁶A writers.

Supplemental Figure S8. Examples of hypomethylated sites with corresponding P values in genes identified in plants with impaired m⁶A writers.

Supplemental Figure S9. MTA and MTB possess m⁶A methyltransferase activity in vitro.

Supplemental Figure S10. Enrichment analysis of hypomethylated genes for cellular components.

Supplemental Figure S11. Functional enrichment analysis of hypomethylated transcripts that are unique in individual genetic backgrounds.

Supplemental Figure S12. Susceptibility of plants with impaired m⁶A writers to salt stress.

Supplemental Figure S13. Top 5 GO terms that are distinct in the upregulated or downregulated gene list (log₂ [fold change] > 0.5 or < −0.5; P < 0.05) from different genetic backgrounds.

Supplemental Figure S14. Functional enrichment analysis of differentially expressed transcripts (log₂ [fold change] > 0.5 or < −0.5; P < 0.05) containing hypomethylated m⁶A site(s) (DMR < −0.3).

Supplemental Figure S15. m⁶A has a minimal effect on alternative splicing of transcripts.

Supplemental Figure S16. Functional enrichment analysis of human transcripts containing hypomethylated m⁶A site(s).

Supplemental Data Set S1. List of hypomethylated sites identified from genetic materials with impaired m⁶A writers.

Supplemental Data Set S2. Differentially expressed genes in genetic materials with impaired m⁶A writers compared to wild-type plants [log2(Fold change) ≤ −0.5 or log2(Fold change) ≥ 0.5; P < 0.05].

Supplemental Data Set S3. Stress-related transcripts that are hypomethylated and upregulated in genetic materials with impaired m⁶A writers compared to wild-type plants.

Supplemental Data Set S4. List of differential alternative splicing events identified from genetic materials with impaired m⁶A writers (P < 0.05).

Supplemental Data Set S5. List of genes with significant alteration in 3′ UTR length identified from genetic materials with impaired m⁶A writers (FDR < 0.05).

Supplemental Data Set S6. List of differential readthrough transcripts identified from genetic materials with impaired m⁶A writers (Fisher's exact test: P < 0.05).

Supplemental Data Set S7. List of primers used in this study.

Supplemental Data Set S8. Human genes with differential methylation peaks [A–E; log2(fold change) < 0 and FDR < 0.05] or distal polyA site usage index (PDUI; F–J; FDR < 0.05 and at least two-fold increase or decrease of PDUIs).

Acknowledgments

The authors thank the Arabidopsis Biological Resource Centre for providing seeds and Genome Institute of Singapore, A*STAR for the nanopore direct RNA sequencing service. The authors thank the Protein and Proteomics Centre (PPC) in the Department of Biological Sciences, National University of Singapore for the LC-MS/MS service.

Funding

This work was supported by the National Research Foundation Competitive Research Programme (NRF-CRP22-2019-0001), Singapore Food Story R&D Programme (SFS_RND_SUFP_001_04), and the intramural research support from Temasek Life Sciences Laboratory and National University of Singapore.

References

Author notes