Deoxyribozyme-based method for absolute quantification of N6-methyladenosine fractions at specific sites of RNA

N6-Methyladenosine (m6A) is the most prevalent modified base in eukaryotic mRNA and long noncoding RNA. Although candidate sites for the m6A modification are identified at the transcriptomic level, methods for site-specific quantification of absolute m6A modification levels are still limited. Herein, we present a facile method implementing a deoxyribozyme, VMC10, which preferentially cleaves the unmodified RNA. We leveraged reverse transcription and real-time quantitative PCR along with key control experiments to quantify the methylation fraction of specific m6A sites. We validated the accuracy of this method with synthetic RNA in which methylation fractions ranged from 0 to 100% and applied our method to several endogenous sites that were previously identified in sequencing-based studies. This method provides a time- and cost-effective approach for absolute quantification of the m6A fraction at specific loci, with the potential for multiplexed quantifications, expanding the current toolkit for studying RNA modifications.

Over 100 types of RNA modifications have been identified to date. Among them, N 6 -methyladenosine (m 6 A) 2 is most prevalent in mRNA and various long noncoding RNA (lncRNA) in higher eukaryotes (1). m 6 A modifications are widely involved in post-transcriptional gene regulation. The complex and dynamic nature of m 6 A-mediated regulation enables timely responses to signaling cues and large-scale modulation of gene expression. Therefore, m 6 A has been shown to be essential for development and associated with many human diseases (2,3). The single methyl group is commonly deposited by either a methyltransferase writer complex composed of METTL3, METTL14, and WTAP (4) or by METTL16 methyltransferase (5) and is removed by either FTO (6) or ALKBH5 demethylase (7). Through its effects on RNA secondary structure and its interactions with m 6 A-binding proteins, m 6 A modifications affect essentially all known steps during an RNA's lifetime, including alternative splicing, polyadenylation, RNA export, translation, and degradation (8,9). Despite m 6 A modification having a consensus DRACH motif (D ϭ A, G, or U; R ϭ G or A; and H ϭ A, C, or U) (10,11), the substoichiometric nature of m 6 A modification potentially creates large compositional heterogeneity in a single RNA species, i.e. each RNA of the same species may selectively carry m 6 A modification at one or a few DRACH motifs among all (12). Being able to quantify the extent of m 6 A modification at precise sites can greatly advance our current understanding of how changes in the m 6 A modification pattern (the site and fraction) are modulated by signaling cues and are then linked to various functional consequences.
Because of its important roles, techniques have been developed and applied to detect and quantify m 6 A modification. Detection of m 6 A modification is primarily facilitated by various high-throughput sequencing-based methods utilizing antibodies and chemical cross-linking (10,11,13,14). Although these sequencing-based methods can map m 6 A candidate sites at the transcriptomic level, they cannot provide the fraction of modification at each site, because of factors such as antibody binding efficiency, specificity, and cross-linking reactivity (15). Real-time quantitative PCR (qPCR) was previously applied for locus-specific detection of pseudouridine () modification through chemical labeling of residue, causing a shift in the melting peak of the resulting qPCR amplicons (16). Similar quantitative methods were recently developed for detection of m 6 A. These methods utilize enzymatic activities followed by qPCR, including differential ligation efficiency of T3 and T4 DNA ligases (17,18), differential reverse transcription activity of Tth and BstI reverse transcriptases (19,20), and a combination of selective elongation of DNA polymerase and ligation (21). Although these polymer elongation and ligation-based methods are successful at modification discrimination and can report the relative m 6 A abundance change, absolute quantification using these methods were only applied on MALAT1, an abundant lncRNA. In addition, the potential sequence dependence of these enzymatic activities requires caution for general applications to these methods (22,23). Considering these potential pitfalls, absolute quantification using these methods would require calibration curves using fully modified and fully unmodified RNA for each target m 6 A site, which is expensive. The only available qPCR-independent method that can provide This work was supported by funds from the Searle Scholars Program, National Institutes of Health Director's New Innovator Award 1DP2GM128185-01, and a pilot grant from the National Institutes of Health Center for the Dynamic RNA Epitranscriptomes at the University of Chicago (to J. F.). The authors declare that they have no conflicts of interest with the contents of this article. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This article contains Tables S1-S4 and Figs. S1-S7. 1 To whom correspondence should be addressed. E-mail: jingyifei@ uchicago.edu. 2 The abbreviations used are: m 6 A, N 6 -methyladenosine; lncRNA, long noncoding RNA; qPCR, quantitative PCR; SCARLET, site-specific cleavage and radioactive labeling followed by ligation-assisted extraction and TLC; DR, deoxyribozyme; GB RNA, 460-nt in vitro transcribed RNA from a gene block sequence with only one adenine in the sequence; dDR, dead DR.
cro ARTICLE absolute quantification of m 6 A fraction site-specifically is sitespecific cleavage and radioactive labeling followed by ligationassisted extraction and TLC (SCARLET) (24). However, the sophistication of the method and its requirement for radioactive labeling prevents its broad application. Very recently, endoribonuclease digestion-based sequencing methods have been developed that rely on selective cleavage of unmethylated A at the ACA motif (25,26). These approaches provide singlebase resolution for identification of modifications site with relative quantitative information but are limited to m 6 A sites carrying the ACA motif, as well as regions that contain relatively sparse ACA motifs (26).
To address these challenges, we present an easy-to-implement method for quantifying m 6 A fraction at specific loci from extracted total RNA. The method utilizes a deoxyribozyme (DR), VMC10, recently reported by Sednev et al. (27) to discriminate between A-and m 6 A-containing RNA (Fig. 1A). We chose VMC10 DR because it cleaves unmethylated A with reasonably high and robust efficiency, whereas its cleavage of m 6 A remains low even after a long incubation time (27). Sednev et al. (27) demonstrated that the remainder of intact RNA after VMC10 DR treatment is correlated with the known methylation level of specific sites of abundant endogenous RNAs. However, quantification of the m 6 A fraction requires additional characterization and correction of potential false positives and false negatives caused by various factors. Without these additional characterizations, modification levels of different m 6 A sites cannot be compared because of the sequence dependence of DR. In our method, we quantify and correct for the potential errors caused by sequence-dependent incomplete cleavage of the DR and the presence of nontarget RNAs from the total RNA. We show that our method can be used to robustly quantify the m 6 A fraction at specific loci on endogenous RNAs with a broad range of cellular abundance. We further show that the method can be adjusted for multiplexed quantification of m 6 A sites. Finally, we extensively discuss the limitations of the method and factors that need to be considered when applying this method.

Specificity and sequence dependence of DR cleavage efficiency
We first verified the cleavage efficiency of DR on a variety of fully modified or unmodified sites. For this purpose, we employed a 460-nt in vitro transcribed RNA from a gene block sequence with only one adenine in the sequence (referred to as "GB RNA" hereafter) and 35-41-mer synthetic RNA fragments with sequences around the MALAT1 2515, MALAT1 2577, and ACTB 1216 sites (Table S1). Each of these targets has either  ). B, unmodified RNA is selectively cleaved by DR upstream of the target site, whereas m 6 A-modified RNA remains uncleaved. The remaining uncleaved RNAs are then quantified using RT with gene-specific reverse primer and qPCR. To control for variations in RNA input, an adjacent region on target RNA is also quantified with RT and qPCR as an internal reference. C, in the negative control sample, RNA is treated with a nonfunctional version of DR (dDR). Both m 6 A modified and unmodified RNA targets remain uncleaved and are subsequently quantified with RT and qPCR.
Absolute quantification of m 6 A fraction m 6 A or A at the respective m 6 A sites. The RNAs were treated with corresponding 40-mer VMC10 DR (referred to as DR for simplification hereafter) and subsequently analyzed by denaturing PAGE. For all the targets, the RNA fragments with unmethylated A were cleaved with high efficiency, and the cleavage efficiencies of modified RNAs were consistently below 5% (Fig. 2, A-D). In addition, the cleavage efficiencies on the unmodified RNAs were sequence-dependent (Fig. 2E), ranging from 50 to 82% for our tested cases.

Method for absolute quantification of m 6 A fraction
We designed a quantification assay using RT-qPCR. As shown in Fig. 1, DR is designed for each modification site based on a VMC10 construct. The total RNA is subjected to DR treatment, during which only unmethylated RNAs upstream of the target site are cleaved. Thus, after the DR treatment, the amount of cleaved RNA should be inversely proportional to the methylation fraction (F m ) of RNA at the target site. The remaining RNA can be quantified using RT-qPCR. To control for any variations in the initial RNA input, we use RT-qPCR to also detect levels of adjacent uncleaved regions on the same target RNA as an internal reference.
Theoretically, the modified fraction calculated from qPCR can be written as, where Ct ϩDRϪm 6 A and Ct ϩDRϪref are the qPCR Ct values at the m 6 A site and a nearby reference site in the DR treated sample, whereas Ct ϪDRϪm 6 A and Ct ϪDRϪref are the Ct values at the m 6 A and the reference site without DR digestion.
However, the measured ⌬⌬Ct only reflects the digested fraction of the RNA substrate; i.e. Equation 1 only holds when digestion efficiency of the unmodified template is 100% and digestion efficiency of the modified template is 0%. Incomplete cleavage of unmodified A will lead to a false positive, and cleavage of the m 6 A will lead to a false negative. Based on the previous study and our tested cases (Fig. 2), the cleavage of VMC10 DR on m 6 A sequence is minimal, leading to insignificant error caused by false negative (Fig. S1). In addition, it is practically difficult and expensive to generate in vitro purified template containing 100% modified m 6 A to account for the exactly falsenegative error at each m 6 A site of interest; we therefore left out the correction factor for false-negative error in our final calculation. On the other hand, the false-positive error can be significant because of the sequence-dependent incomplete cleavage of unmodified RNA by the DR and needs to be corrected for each m 6 A site of interest.
We therefore consider two major factors that may contribute to the false-positive error caused by incomplete digestion and quantify the effect of the two factors to extract the true modification fraction: the intrinsic sequence-dependent digestion efficiency and the presence of a large amount of nontarget RNAs from the total RNA extract. We define F DR as a correction factor to account for the incomplete DR digestion efficiency, which has to be determined for each m 6 A target (Fig. 2). We can determine F DR at each m 6 A site of interest by performing the DR digestion followed by RT-qPCR using the in vitro transcribed unmodified RNA, in which ⌬⌬Ct is determined as in Equation 2. We define F N as the ratio of DR digestion efficiency of an RNA target in total RNA over digestion efficiency of a pure RNA target, to account for the potential drop of DR efficiency caused by the presence of nontarget RNAs. We can determine F N by performing DR

Absolute quantification of m 6 A fraction
digestion using the same in vitro transcribed unmodified RNA mixed with total RNA and compare with F DR from Equation 3, in which ⌬⌬Ct is determined as in Equation 2. With the quantification of F DR and F N , the corrected modification fraction follows.
We therefore can calculate F m as follows.

Validation of the absolute quantification of m 6 A fraction using pure RNA
To test the feasibility of the method to quantify the m 6 A methylation fraction, we used GB RNA with methylation fractions ranging from 0 to 100%. We performed DR treatment on GB RNA and estimated the cleaved fractions by denaturing PAGE. The cleaved fraction was linearly dependent on the methylation fraction of the input RNA (Fig. 3, A and B). Next, we tested whether we can use RT-qPCR to quantify the absolute methylation fraction. Because this quantification is performed on in vitro purified RNA, only F DR is needed to correct for F m . Based on Equation 3, using the 100% unmodified GB RNA, we measured Ct ϩDRϪm 6 A and Ct ϩDRϪref at the m 6 A site and a nearby reference site in the DR-treated sample and Ct ϪDRϪm 6 A and Ct ϪDRϪref at the m 6 A and the reference site using a negative control containing identical amount of RNA but without the DR. We found that in addition to the expected larger Ct ϩDRϪm 6 A compared with Ct ϪDRϪm 6 A , there was a consistent difference between Ct ϩDRϪref and Ct ϪDRϪref . We speculated that this difference in Ct values at the reference site might be due to changes in the RNA secondary structure upon DR binding that can affect RT efficiency. To create a more accurate negative control, we designed a nonfunctional version of DR ("dead" DR or dDR) (Fig. 1, A and C), which has mutations in the AGC triplet, CG dinucleotide, and position 19 important for the catalytic activity of 8 -17 family of enzymes (28). We tested the activity of dDR on multiple targets, for all of which digestion of RNA was undetectable (Figs. S2 and S3). Indeed, using the dDR-treated RNA as a negative control, the difference between Ct ϩDRϪref and Ct ϪDRϪref was eliminated (Fig. S4). Based on Equation 3, we determined F DR of the synthetic RNA to be 0.49 Ϯ 0.08 (mean Ϯ S.D.). With the F DR correction, we showed that the estimated F m correlated well with the input m 6 A methylation fractions ( Fig. 3C and Fig. S5).

DR cleavage efficiency in presence of nontarget RNAs
Next, we evaluated how the presence of total RNA affects the cleavage efficiency of the DR. The presence of the large amount of nontarget RNAs may compete for DR binding, consequently decreasing its cleavage efficiency at the target site in total RNA as opposed to purified RNA. We accounted for this potential decrease in efficiency with the F N correction factor, which we measured using three RNA transcripts: 1) the GB RNA used above, which is naturally missing in total RNA; 2) a PLAC2 RNA fragment containing two target sites, which is of low abundance in HeLa cell line; and 3) an unmethylated A site in the endogenous ACTB mRNA. The A1165 site on ACTB mRNA was chosen as the unmethylated A site because it was not detected in the sequencing-based studies (4,29), nor does it contain the DRACH consensus motif. F DR of these three RNAs were measured with in vitro transcribed RNAs based on Equation 3 (Fig.  4, A and B). Then the in vitro transcribed GB RNA and the PLAC2 RNA fragment were spiked into the total RNA respectively to determine F N based on Equation 4. For the unmethylated A site in the ACTB mRNA, F N was determined by measuring the total RNA directly.
To increase the binding specificity of DR, we compared a 60-mer DR and a 40-mer DR. We found that F DR values of 60-mer DR were higher than those of 40-mer DR (Fig. S2 and Fig. 4, A and B), likely because of a higher hybridization efficiency by 60-mer DR. F N values were consistently high for all tested RNAs, with the lowest F N values being 0.78 Ϯ 0.02 for 40-mer DR and 0.93 Ϯ 0.02 for 60-mer DR, demonstrating that the ability of our method to quantify m 6 A status should not be compromised by the presence of total RNA and that 60-mer can slightly outperform 40-mer DR (Fig. 4C). Overall, the average F N values were determined to be 0.94 Ϯ 0.1 for 40-mer DR and 0.98 Ϯ 0.05 for 60-mer DR. In addition, we compared the effect of DR concentration on F N values using PLAC2 m 6 A 2 DR as an example. We found that the cleavage efficiency is consis-

Absolute quantification of m 6 A fraction
tent at a wide range of DR concentrations as long as the molar concentration of DR is at least 10-fold higher than the molar concentration of RNA target (Fig. S6A). Likewise, the F N values stayed consistently at ϳ1 for all tested concentrations of DR that were saturated compared with target RNA (Fig. S6B). Given these results, we can simplify Equation 6 to be as follows.

Quantification of m 6 A fraction of endogenous sites
Having developed and validated our method, we applied it to determine the methylation fraction of several endogenous sites that were identified as potential m 6 (4,14,29). The selected RNAs vary from low to high abundance in HeLa cells, and some of them contain more than one modification site. Specifically, four of the targets (MALAT1 2515, MALAT1 2577, and MALAT1 2611, and ACTB 1216) were previously measured using the SCARLET assay (24), therefore serving as an additional validation of our method. To apply our method, a DR and a dDR were designed for each site. Because of the higher F DR and F N values with 60-mer DR, we chose to use the 60-mer DR for all endogenous RNAs. For each target site, we first generated in vitro transcribed RNAs containing the m 6 A sites of interest and performed DR digestion on these in vitro transcribed unmethylated RNAs to get F DR for each site. The F DR values were all greater than 0.49 and again varied among different RNAs (Fig. 5A and Fig. S3).
The methylation fractions of the endogenous sites were determined to range from 0.13 to 0.92 (Fig. 5B). Notably, our results show comparable methylation fractions for MALAT1 2515, MALAT1 2577, MALAT1 2611, and ACTB 1216 as in the SCARLET assay (24). Although the generally consistent results between our methods and SCARLET assay help validate our assay, we did notice that the values measured in our assay are slightly higher than those from SCARLET. One possible explanation for this slight variation can be the splint ligation step used in the SCARLET assay, in which the DNA oligonucleotide needs to be ligated to the RNase H cleaved RNA carrying either unmodified A or m 6 A at the 5Ј end (24). It is possible that the splint ligation is less efficient for the m 6 A-containing RNA and therefore underestimates the m 6 A fraction in SCARLET assay.
Having validated that the method is able to quantify m 6 A fractions one site at a time, we investigated whether the method can be utilized in a multiplexed way. As a proof of concept, we chose to remeasure three m 6 A sites with varying levels of methylation: MALAT1 2611, ACTB 1216, and INCENP 912. For each biological replicate, the three corresponding active DRs were combined in one reaction, and the three dDRs were combined in another reaction. Similarly, RT reactions were also performed with combined RT primers for all three targets at either m 6 A site or internal control site, followed by separate qPCR for each target RNA. The multiplexed measurements of the methylation fractions were comparable with the ones we previously obtained in individual measurements (Fig. 5C). These results support that the method is accurate at measuring m 6 A fractions in a multiplexed fashion and can be further adapted for highthroughput measurements.

Effect of the nearby modifications on the DR cleavage efficiency
m 6 A modifications often exist in clusters (10). In addition, other types of RNA modifications are identified in mRNAs and lncRNAs (2). Therefore, the possible effect of the nearby modifications needs to be considered when applying this method. We designed synthetic RNA containing a nearby m 6 A, m 1 A, or and measured their effects on cleavage efficiency of DR by PAGE analysis (Fig. 6 and Fig. S7). The cleavage efficiency was unaffected by the presence of m 6 A modification 2 and 4 nt upstream and downstream of the target site, suggesting that the method can be used to quantify the m 6 A fractions in RNA that

Absolute quantification of m 6 A fraction
contain m 6 A modifications in clusters. Furthermore, modification caused a very minimal decrease in the cleavage efficiency of DR 2 nt away from target site and had no effect on the cleavage when present 4 nt away from the target site. Finally, m 1 A modification, which can affect the Watson-Crick base pairing, significantly decreased the cleavage efficiency at 2 nt away from the target site and moderately decreased the cleavage efficiency at 4 nt away from the target site. Overall, the results indicate that other nearby RNA modifications that do not affect the base pairing with the DR are not likely to affect the DR cleavage efficiency even when placed as close as only 2 nt away from the target site. However, nearby RNA modifications that weaken the base pairing with the DR will have a larger effect on the DR activity, but the effect decreases when the modification is more distal from the target site.

Discussion
In summary, here we present a method for quantifying the absolute methylation fraction of potential m 6 A sites using a previously developed VMC10 DR (27), expanding the toolkit for site-specific quantification of m 6 A. In addition, the method can be adjusted for high-throughput quantification of m 6 A sites. Furthermore, as the VMC10 DR selectively cleaves the unmodified A, it can potentially be used to discriminate other modifications, such as m 1 A (30). We therefore expect that the DR-based quantification method can be easily applied to sitespecific absolute quantifications of other RNA modifications.
Although this method is easy to implement, there are several limitations that need to be considered. First, the assay utilizes VMC10 DR, which has high cleavage efficiencies only on DGACH sequences, limiting its application on a subset of m 6 A

Absolute quantification of m 6 A fraction
sites with the DAACH sequences (27). Second, DR digestion efficiency varies among different sequences. Although low DR cleavage efficiency can be corrected by determining F DR for each modification site of interest using in vitro transcribed RNA, low DR efficiency can lead to less accurate quantification because of two reasons. First, a higher digestion efficiency leads to a larger ⌬⌬Ct that reduces the measurement variation by qPCR. Conversely, low digestion efficiency will make the ⌬⌬Ct too small to be accurately detected by qPCR. Second, the Ͻ5% cleavage efficiency on the modified RNA can lead to underestimation of the m 6 A fraction, and the percentage of underestimation depends on F DR (Fig. S1). A lower F DR will result in a larger underestimation. When F DR is 50%, a 5% cleavage of the modified RNA will result in a 10% underestimation of the m 6 A. Finally, the presence of a nearby modified nucleotide may affect the DR cleavage efficiency, depending on the type of the modification (Fig. 6).
To improve the accuracy of the measurement, there are also a few factors to note. First, for the synthetic RNA, we observed equal quality of F m estimation using samples treated with dDR or samples lacking any DR as a negative control ( Fig. 3C and Fig.  S5). Nevertheless, we still recommend using dDR-treated sample as a negative control, because it corrects for potential changes in the RNA secondary structure caused by DR binding that can affect RT efficiency. Second, we recommend using 60-mer DR for quantification, because 60-mer DR overall has higher digestion efficiencies of unmethylated RNAs potentially because of a higher hybridization efficiency. Third, the quality of the primers used for RT and qPCR should be verified by performing calibration curves. Finally, we noticed that the largest source of technical variability in measurements originates from the RT step (compare error bars in Fig. 3, B and C). We therefore recommend performing multiple RT reactions for each DR-treated sample for a more accurate characterization.

Cell culture and RNA extraction
HeLa cells were cultured in Dulbecco's modified Eagle's medium (Gibco) supplemented with 10% (v/v) fetal bovine serum (Thermo Fisher Scientific). The cells were grown at 37°C under humidified conditions with 5% CO 2 . The total RNA was extracted using the RNeasy mini kit (Qiagen) according to the manufacturer's instructions.

In vitro transcription of endogenous RNA target fragments
The dsDNA templates for in vitro transcription were prepared by PCR with primers (Integrated DNA Technologies) that contain T7 promoter sequence and cDNA generated from total HeLa RNA. 1 g of dsDNA templates were added to 100-l reactions containing final concentrations of 2.5 mM each rNTP (New England Biolabs), 1ϫ T7 polymerase reaction buffer (New England Biolabs), 4 mM MgCl 2 , 0.5 unit/l SUPERase-In RNase inhibitor (Thermo Fisher Scientific), and 14 units/l T7 polymerase (a kind gift from Dr. D. Bishop's Group). The reactions were incubated at 37°C for 1 h. The transcript products were treated with DNase I recombinant (Roche) at 37°C for 30 min and purified by phenol-chloroform extraction and ethanol precipitation. The primers are listed in Table S2.

In vitro transcription of 0 and 100% methylated GB RNA
A gene block containing a 460-nt random sequence with 51% GC content and one adenosine was purchased from Genewiz. The gene block sequence is listed in Table S1. The dsDNA template was amplified with primers containing an upstream T7 promoter sequence (Table S2). The in vitro reactions were carried out in the same conditions as for endogenous RNA targets, except that N 6 -methyladenosine-5Ј-triphosphate (Trilink Biotechnologies) was used instead of rATP for the generation of 100% methylated RNA. The transcript products were treated with DNase I recombinant (Roche) at 37°C for 30 min and purified by phenol-chloroform extraction, 7% denaturing polyacrylamide gel, and ethanol precipitation.

Synthesis of RNA oligonucleotide
Unmodified phosphoramidites were purchased from Glen Research. Phosphoramidite of N 6methyladenosine was synthesized by following previously published procedure (17). RNA oligonucleotides were synthesized using Expedite DNA synthesizer at 1 mol scale. After deprotection, RNA oligonucleotides were purified by PAGE.

Deoxyribozyme digestion
Total RNA (500 ng to 2 g) or in vitro transcribed RNA fragments (50 nM) were mixed with 55.6 M of either DR or dDR (Integrated DNA Technologies), 55.6 mM Tris-HCl (pH 7.5), and 166.7 mM NaCl in a final volume of 9 l. The annealing of DR to the target site was facilitated by 5 min of incubation at 95°C, followed by slow cooling to room temperature. After annealing, 1 l of 200 mM MgCl 2 was added to each reaction and incubated at 37°C for 12 h. The final concentrations of the reagents in the incubation buffer are 50 M of DR, 50 mM Tris-HCl (pH 7.5), 150 mM NaCl, and 20 mM MgCl 2 in a 10-l reaction. The DR treatment of 35-40-nt MALAT1 2515, MALAT1 2577, and ACTB 1216 was carried out following the same protocol, except with stepwise cooling (95°C for 5 min and 25°C for 10 min) instead of slow cooling. To remove the DR after the digestion, 1.33 l of 10ϫ TURBO DNase buffer and 2 l of TURBO DNase (Thermo Fisher Scientific) were added to the 10-l DR reactions. The samples were incubated at 37°C for 2 h. Subsequently, the DNase enzyme was inactivated by addition of EDTA (pH 7.5) to 15 mM final concentration and incubated at 75°C for 10 min. The DR sequences are listed in the Table S3.

RNA digestion analyzed by PAGE
DR digestion reactions containing 50 -100 ng of in vitro transcribed RNA or 35-40-nt synthetic RNAs were run on either 7% or 15% denaturing (7 M urea) PAGE, respectively. The gels were stained with SYBR green II RNA gel stain (Thermo Fisher Scientific) for 10 min and imaged with ChemiDoc TM Imaging System (Bio-Rad). The cleavage efficiencies were analyzed with ImageJ using intensities of bands corresponding to the fulllength RNA and the longer cleaved product.

Reverse transcription
For each DR-treated sample, separate RT reactions were performed with gene-specific reverse primers for the m 6 A region and internal reference site. Because of the presence of excess EDTA after DNase inactivation, either the reactions were significantly diluted or extra MgCl 2 was added for maximum reverse transcriptase activity. The RNA was denatured at 70°C for 5 min and then added to freshly prepared RT buffer with final concentration of 1 mM dNTPs (Thermo Fisher Scientific), 10% DMSO (Thermo Fisher Scientific), 10 mM DTT (Sigma-Aldrich), 250 nM of gene-specific reverse primer (Integrated DNA Technologies), and 20-fold dilution of reverse transcriptase from an iScript cDNA synthesis kit (Bio-Rad). The reactions were incubated at 25°C for 5 min and at 46°C for 20 min and heat-inactivated at 95°C for 1 min. All primers are listed in Table S4. qPCR 1 l of cDNA was added into reaction mixture, containing 250 nM of each forward and reverse primers and 1ϫ SsoAdvanced TM Universal SYBR Green Supermix (Bio-Rad) in a final volume of 20 l. The qPCRs were performed with CFX real-time PCR system (Bio-Rad), using preincubation of 95°C for 30 s, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. The reactions were then subjected to melting curve analysis: 95°C for 10 s and 65°C for 5-s increments by 0.5°C to 95°C for 5 s. The data were analyzed with the supporting Bio-Rad CFX Maestro software. All primers are listed in Table S4. All error bars in the figures are means Ϯ S.D. of multiple biological replicates. For in vitro prepared GB RNA with different input m6A fraction, biological replicates are defined as independently mixed GB RNA samples. For the cases of endogenous mRNAs, biological replicates are defined as independently extracted total RNA samples. The m 6 A fraction calculated for each biological replicate is from the average values of multiple technical replicates defined by independently performed RT reactions for each RNA sample.

Data availability
All data are contained within the manuscript and the supporting information file.