Deoxyribozyme-based Method for Site-specific Absolute Quantification of N6-methyladenosine Modification Fraction

N6-methyladenosine (m6A) is the most prevalent modified base in eukaryotic messenger RNA (mRNA) and long noncoding RNA (lncRNA). Although candidate sites for m6A modification are identified at the transcriptomic level, site-specific quantification methods for m6A modifications are still limited. Herein, we present a facile method implementing deoxyribozyme that preferentially cleaves the unmodified RNA. We leverage reverse transcription and real-time quantitative PCR along with key control experiments to quantify the absolute methylation fraction of specific m6A sites. We validate the accuracy of the method using synthetic RNA with controlled methylation fraction and apply our method on 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, expanding the current toolkit for studying RNA modifications.


INTRODUCTION
Over 100 types of RNA modifications have been identified to date. Among them, N 6methyladenosine (m 6 A) is most prevalent in messenger RNA (mRNA) and various long noncoding RNA (lncRNA) in higher eukaryotes (Boccaletto et al. 2018). m 6 A modifications are widely involved in post-transcriptional gene regulation. The complex and dynamic nature of m 6 Amediated 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 (Nachtergaele and He 2018;Maity and Das 2016). The single methyl group is commonly deposited by either a methyltransferase writer complex composed of METTL3, METTL14, and WTAP (Liu et al. 2014) or by METTL16 methyltransferase (Warda et al. 2017) and is removed by either FTO (Jia et al. 2011) or ALKBH5 demethylase .
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 (Kasowitz et al. 2018;Shi et al. 2019). Despite m 6 A modification having a consensus DRACH motif (D=A, G or U; R=G or A; H=A, C or U) (Meyer et al. 2012;Dominissini et al. 2012), the sub-stoichiometric 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 . Being able to quantify m 6 A modification fraction 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.
Due to 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 highthroughput sequencing-based methods utilizing antibodies and chemical crosslinking (Meyer et al. 2012;Dominissini et al. 2012;Chen et al. 2015;Linder et al. 2015). Although these sequencingbased methods can map m 6 A candidate sites at the transcriptomic level, they cannot provide the fraction of modification at each site, due to factors such as antibody binding efficiency, specificity and cross-linking reactivity (Helm and Motorin 2017). Real-time quantitative PCR (qPCR) was previously applied for locus specific detection of pseudouridine () modification though chemical labelling of  residue, causing a shift in the melting peak of the resulting qPCR amplicons (Lei and Yi 2017). 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 (Dai et al. 2007;Liu et al. 2018), differential reverse transcription activity of Tth and BstI reverse transcriptases (Harcourt et al. 2013;Wang et al. 2016), and a combination of selective elongation of DNA polymerase and ligation (Xiao et al. 2018). 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 sequencedependence of these enzymatic activities requires caution for general applications to these methods (Potapov et al. 2018;Harada and Orgel 1993). 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 qPCRindependent method that can provide absolute quantification of m 6 A fraction site-specifically is SCARLET (site-specific cleavage and radioactive labeling followed by ligation-assisted extraction and TLC) (Liu et al. 2013). However, the sophistication of the method and its requirement for radioactive labeling prevents its broad application. Very recently, endoribonuclease digestionbased sequencing methods have been developed, which rely on selective cleavage of unmethylated A at the ACA motif (Zhang et al. 2019;Garcia-Campos et al. 2019). These approaches provide single-base 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 (Garcia-Campos et al. 2019). To address these challenges, we present an easyto-implement method for quantifying m 6 A fraction at specific loci from the extracted total RNA using a deoxyribozyme (DR) with strong preference for unmethylated RNA.

Specificity and sequence-dependence of DR cleavage efficiency
The method utilizes recently reported DR to discriminate between A and m 6 A containing RNA (Sednev et al. 2018). However, to accurately determine the modification fraction, it is important to have high cleavage efficiency for one form, and minimal for the other form to reduce the potential false positive and false negative. In addition, in order to achieve reproducible quantification, it is preferred to have discrimination efficacy less sensitive to the reaction condition.
We therefore chose VMC10 DR (Fig. 1A), as it can cleave unmethylated A with reasonably high efficiency while the cleavage of m 6 A remains low even after long incubation time (Sednev et al. 2018).
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 as "GB RNA" hereafter), and 35 to 41-mer synthetic RNA fragments with sequences around MALAT1 2515 site, MALAT1 2577 site, and ACTB 1216 site (Supplemental Table S1). Each of these targets has either m 6 A or A at the respective m 6 A sites. The RNAs were treated with corresponding 40-mer DR and subsequently analyzed by denaturing polyacrylamide gel electrophoresis (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. 2A-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 reverse transcription (RT) and 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 deoxyribozyme treatment, the amount of cleaved RNA should be inversely proportional to the methylation fraction (Fm) of RNA at the target site. The remaining RNA can be quantified using RT-qPCR. In order to control for the initial RNA input, we use RT-qPCR to also detect levels of adjacent uncleaved regions on the target RNA as an internal reference.
Theoretically, the modified fraction calculated from qPCR can be written as

2
(1), in which However, the measured only reflects the digested fraction of the RNA substrate, i.e., Eq.
(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 false positive, and cleavage of the m 6 A will lead to 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 (Supplemental 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 false negative 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 due to 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 due to 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 non-target RNAs from the total RNA extract. We define FDR as a correction factor to count for the incomplete DR digestion efficiency, which has to be determined for each m 6 A target (Fig. 2).
We can determine FDR 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 Eq.
(2). We define FN 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 due to the presence of non-target RNAs. We can determine FN by performing DR digestion using the same in vitro transcribed unmodified RNA mixed with total RNA, and compare with FDR from Eq. (3): in which ΔΔCt is determined as in Eq.
(2). With the quantification of FDR and FN, the corrected modification fraction follows: We therefore can calculate Fm as: (6).

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. 3A,B). Next, we tested whether we can use RT-qPCR to quantify the absolute methylation fraction. As this quantification is performed on in vitro purified RNA, only FDR is needed to correct for Fm. Based on Eq.
(3), using the 100% unmodified GB RNA, we measured and at the m 6 A site and a nearby reference site in the DR treated sample, and and 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 compared to , there was a consistent difference between and . 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 non-functional version of DR ("dead" DR or dDR) (

DR cleavage efficiency in presence of nonspecific RNAs
Next, we evaluated how the presence of total RNA affects the cleavage efficiency of the DR. The presence of the large amount of non-specific 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 FN 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)  (3) (Fig. 4A,B). Then the in vitro transcribed GB RNA and the PLAC2 RNA fragment were spiked into the total RNA respectively to determine FN based on Eq. (4). For the unmethylated A site in the ACTB mRNA, FN was determined by measuring the total RNA directly. In order to increase the binding specificity of DR, we also compared a 60-mer DR and a 40-mer DR. We found that FDR values of 60-mer DR were higher than those of 40-mer DR (Supplemental Fig. S2), likely due to a higher hybridization efficiency by 60-mer DR. FN values were consistently high for all tested RNAs, with the lowest FN values being 0.78 ± 0.02 for 40-mer DR and 0.93 ± 0.02 for 60-mer DR, demonstrating that 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 FN values were determined to be 0.94 ± 0.1 for 40-mer DR and 0.98 ± 0.05 for 60-mer DR. Therefore, we can simplify Eq. (6) to be (7).

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 A sites by RNA sequencing from more than one study: MALAT1 2515 (chr11 65500276 The methylation fractions of the endogenous sites were determined to range from 0.20 to 1.0 (Fig. 5B). Notably, four of the targets (MALAT1 2515, MALAT1 2577, and MALAT1 2611 and ACTB 1216) were previously measured using the SCARLET assay (Liu et al. 2013), and our results show comparable methylation fractions. While the generally consistent results between our methods and SCARLET assay help validate our assay, we did notice that 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 oligo needs to be ligated to the RNase H cleaved RNA carrying either unmodified A or m 6 A at the 5' end (Liu et al. 2013). 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.

Potential limitations and other considerations of the method
In summary, we have established a method for quantifying the absolute methylation fraction of potential m 6 A sites of specific transcripts using deoxyribozyme digestion, expanding the toolkit for site-specific quantification of m 6 A. As the VMC10 DR selectively cleaves the unmodified A, it can potentially be used to discriminate other modifications, such as m 1 A (Tserovski et al. 2016).
We, therefore, expect the deoxyribozyme-based quantification method can be easily applied to site-specific absolute quantifications of other RNA modifications.
While this method is easy to implement, there are several limitations that need to be considered. Firstly, the assay utilizes VMC10 DR, which has high cleavage efficiencies only on DGACH sequences, limiting its application on a subset of m 6 A sites with the DAACH sequences (Sednev et al. 2018). Secondly, DR digestion efficiency varies among different sequences.
Although low DR cleavage efficiency can be corrected by determining FDR for each modification site of interest using in vitro transcribed RNA, low DR efficiency can lead to less accurate quantification due to two reasons. (1) 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.
(2) 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 FDR (Supplemental Discussion and Supplemental Fig. S1). A lower FDR will result in a larger underestimation. When FDR 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.
To test the effect of the nearby modifications, 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 ( Figure 6, Supplemental Figure S6). 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 m 6 A fraction in RNA that contain m 6 A modifications in clusters.
Furthermore,  modification had 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.
To improve the accuracy of the measurement, there are also a few factors to note. Firstly, for the synthetic RNA, we observed equal quality of Fm estimation using samples treated with dDR or samples lacking any DR as a negative control ( Fig. 2C; Supplemental 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. Secondly, we recommend using 60-mer DR for quantification, as 60-mer DR overall has higher digestion efficiencies of unmethylated RNAs potentially due to a higher hybridization efficiency. Thirdly, 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 variability in measurements originates from the RT step (comparing error bars in Figure 3B and C). We, thereby, recommend performing multiple RT reactions for each DR treated sample to reduce measurement error.

Cell culture and RNA extraction
HeLa cells were cultured in DMEM 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% CO2. The total RNA was extracted using 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  Supplemental Table S2.

In vitro transcription of 0% and 100% methylated GB RNA
A gene block containing 460 nt random sequence with 51% GC content and one adenosine was purchased from Genewiz. The gene block sequence is listed in Supplemental 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 oligo
Unmodified phosphoramidites were purchased from Glen Research. Phosphoramidite of N 6methyladenosine was synthesized by following previously published procedure (Dai et al. 2007).
RNA oligos were synthesized using Expedite DNA synthesizer at 1 umol scale. After deprotection, RNA oligos were purified by PAGE.

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

Reverse transcription
For    Determined m 6 A modification fractions of the 12 endogenous sites. All error bars report mean ± s.d. for at least 3 biological replicates. for 3 independent DR cleavage reactions.