Detailed Sequencing Methods to Measure mRNA Chemical Changes Across the Whole Transcriptome

Among all internal mRNA modifications in mammals, m⁶A is the most abundant, exhibiting ∼0.4–0.6% m⁶A/A abundance in mRNA.^14,15 m⁶A affects almost every aspect of mRNA processing and metabolism, impacting a variety of biological processes in a living cell.^5,7 m⁶A methylation on mammalian mRNA is mediated by the METTL3/METTL14 methyltransferase complex as the main m⁶A “writer” protein.¹⁶ FTO and ALKBH5 were identified by our laboratory as the m⁶A “eraser” proteins in 2010 and 2013,^17,18 respectively, which reverse m⁶A methylation. Proteins that preferentially bind m⁶A, or “reader” proteins, were also discovered, with YTHDF1–3^{19 −21} and YTHDC1–2^{22 −24} possessing an evolutionally conserved YTH domain that binds selectively to m⁶A. These proteins bind m⁶A-methylated RNA and regulate mRNA stability, translation, pre-mRNA splicing, nuclear export, and other biological processes. The discovery of m⁶A “writer,” “eraser,” and “reader” proteins largely promote the recent explosion of m⁶A epitranscriptome research; the application of next-generation sequencing technology to map mRNA m⁶A distribution has been an indispensable part in nearly every study related to m⁶A biology.

For over a decade, the antibody-based m⁶A-seq or m⁶A-MeRIP-seq strategy^14,15 has been a prevalent method to study the distribution of m⁶A methylations at 100–200 nucleotide (nt) resolution. Relying on the anti-m⁶A antibody, miCLIP²⁵ and m⁶A-LAIC-seq²⁶ provided updated versions. However, these antibody-based methods were unable to provide single-base resolution, m⁶A stoichiometry, nor the sensitivity to compare m⁶A methylation dynamics. MAZTER-seq²⁷ and m⁶A-REF-seq,²⁸ exploring MazF RNase, selectively cleave RNA at ACA motif without m⁶A methylation. Although this approach could detect m⁶A methylation within the ACA motif, it only accounts for around 15% of the overall m⁶A sites in mRNA and lacks sensitivity. DART-seq,²⁹ m⁶A-SEAL,³⁰ and m⁶A-label-seq³¹ have also been reported, but they either lack stoichiometric information at m⁶A-modified sites or cannot be applied transcriptome-wide.

We have recently reported m⁶A-selective allyl chemical labeling and sequencing (m⁶A-SAC-seq) that reads out m⁶A as mutated signals.¹ Our approach is based on the unique activity of the Methanocaldococcus jannaschii homologue MjDim1, which is known to convert A to m⁶A and then m⁶A to m⁶₂A.³² We found that this enzyme can employ a chemically modified allylic-SAM to install an allyl group at the N6 position of m⁶A. With the allylic-SAM as the cofactor, MjDim1 exhibited a ∼10-fold preference for m⁶A over A in the allyl group transfer reaction, converting m⁶A into allyl-modified m⁶A (N⁶-allyl,N⁶-methyladenosine, a⁶m⁶A; Figurea).¹ To further induce misincorporation signatures at the generated a⁶m⁶A sites, the subsequent I₂ treatment converts a⁶m⁶A and a⁶A into the corresponding N¹,N⁶-ethanoadenine and N¹,N⁶-propanoadenine derivatives, respectively, which can be read out as misincorporation signatures by human immunodeficiency virus 1 (HIV-1) reverse transcriptase (HIV RT; Figurea). HIV RT generated ∼10-fold higher mutation rates at the cyclized a⁶m⁶A sites (m⁶A sites) than the cyclized a⁶A sites (unmodified A sites) in almost all sequence contexts in NNm⁶ ANN versus NNANN. Meanwhile, HIV RT could also read through cyclized NNa⁶m⁶ ANN without noticeable RT stops. Therefore, the MjDim1-catalyzed allyl transfer using allylic-SAM shows a ∼10-fold preference for m⁶A over A, and the cyclized a⁶m⁶A adduct (generated from m⁶A) induces another ∼10-fold higher misincorporation rate than the cyclized a⁶A adduct formed from unmodified A. Overall, m⁶A-SAC-seq offers a ∼100-fold higher selectivity at m⁶A over unmodified A sites¹. Before MjDim1 labeling, RNA fragments could also be split into two sections, “untreated” and “FTO treated”, in which the FTO treatment erases a large portion of mRNA m⁶A as a background control to further eliminate false positives.

m⁶A-SAC-seq maps internal m⁶A methylomes at base resolution, and the misincorporation rates readout at each modified site could be used for assessing the m⁶A modification stoichiometry. By applying m⁶A-SAC-seq to the mixtures of synthetic NNm⁶ANN and NNANN oligo probes in different ratios, we established a spike-in calibration system for each sequencing library and monitored m⁶A modification fractions transcriptome-wide through computation. The initial m⁶A-SAC-seq protocol starts with ∼30 ng of poly(A) RNA or rRNA-depleted RNA. After further optimization of library preparation and bioinformatic pipelines, the latest m⁶A-SAC-seq starts with ∼2 ng of poly(A) RNA with high reproducibility,³³ revealing ∼30 000–130 000 m⁶A sites that overlap well with m⁶A profiles obtained by antibody-based approaches. For the first time, m⁶A-SAC-seq demonstrated the quantitative base-resolution maps of internal m⁶A modifications at diverse motif contexts unbiasedly, and set the technological basis for sensitively monitoring m⁶A dynamics in diverse biological processes.

As a complement sequencing tool to m⁶A-SAC-seq, which involves enzymic reactions targeting m⁶A bases over unmodified adenosines, through collaboration we have also helped develop evolved TadA-assisted N⁶-methyladenosine sequencing (eTAM-seq).⁹ eTAM-seq employs highly efficient global adenosine deamination of unmethylated A sites mediated by TadA8.20 and reads out these unmodified sites as guanosines (G) in next-generation sequencing (NGS) data, while m⁶A is resistant to this enzymic deamination and thus is still read as A (Figureb). In this way, eTAM-seq achieves transcriptome-wide, base-resolution detection and quantification of m⁶A. It uncovered ∼35 000 m⁶A-modified sites in HeLa mRNA,⁹ which overlaps very well with m⁶A-SAC-seq results. Note that eTAM-seq enables site-specific and quantitative sequencing of m⁶A with as few as 10 cells,⁹ which requires much lower input RNA than other existing quantitative sequencing methods. Similar to the deamination principle in eTAM-seq but employing chemical reactions instead, Liu et al. developed glyoxal and nitrite-mediated deamination of unmethylated adenosines (GLORI) to quantitatively map m⁶A at base precision.³⁴ GLORI maps m⁶A methylomes transcriptome-wide in mouse and human cells, revealing clustered m⁶A dynamics with information on site distribution and modification stoichiometry. The current version requires several hundred nanograms of RNA as the starting material, but further improvements should be able to lower the input amount substantially.

m⁶A-SAC-seq has two advantages over deamination-based methods such as eTAM-seq⁹ and GLORI:³⁴ (i) only the m⁶A sites show mutation signatures, which preserves the sequence complexity and allows more accurate sequence mapping; (ii) the positive readout of m⁶A significantly reduces sequencing costs. A main limitation is the sequence context preference of the MjDim1 enzyme. Calibration probes are highly recommended for each sequencing library to reveal modification stoichiometry at different sequence contexts, since MjDim1 does exhibit sequence preference toward the GA motif. Note that while m⁶A-SAC-seq requires much less sequencing depth compared with eTAM-seq and GLORI, sequence bias needs to be calibrated. Taken together, m⁶A-SAC-seq, eTAM-seq, and GLORI currently serve as the three quantitative mapping tools targeting m⁶A (Table S1↗).

Following m⁶A, pseudouridine (Ψ) is the second most abundant mRNA modification in mammalian mRNA, with ∼0.2% Ψ/U levels measured by mass spectrometry. Ψ modifications are known to distribute broadly in abundant noncoding RNAs, such as rRNA, tRNA, and small nuclear RNA (snRNA). Ψ is also known to exist in mammalian mRNA; however, sequencing Ψ within low-abundance mRNA was challenging using traditional experimental methods. Based on the chemical reaction with N-cyclohexyl-N′-(2-morpholinoethyl) carbodiimide methyl-p-toluenesulfonate (CMC) to yield CMC-modified Ψ and subsequent RT truncation signatures induced during reverse transcription,^{35 −37} several next-generation sequencing methods were developed to map the transcriptome-wide distribution of cellular Ψ modifications, including Pseudo-seq,³⁵ Ψ-seq,³⁶ and PSI-seq.³⁷ The RT truncation signatures are difficult to detect, especially at low- and medium-modified sites. Only modest numbers of Ψ sites were identified previously in mammalian mRNA with a low overlap among different data sets. An azide-modified CMC was applied to enrich Ψ-containing mRNA fragments in CeU-seq,³⁸ with many more Ψ sites detected, but this method could not reveal Ψ stoichiometry transcriptome-wide. HydraPsi-seq,³⁹ a method that relies on the resistance of pseudouridine toward hydrazine/aniline cleavage, was also reported to map Ψ modifications in yeast mRNA, but it cannot reveal Ψ stoichiometry either. The lack of a reliable method to map Ψ transcriptome-wide has been a bottleneck in functional investigations of Ψ in mRNA and other low-abundance RNA species.

We reported bisulfite-induced deletion sequencing (BID-seq)² to quantitatively map Ψ at single-base resolution in 2022. Inspired by Khoddami et al., who reported RBS-seq (a modification of RNA bisulfite sequencing),⁴⁰ a modified version of RNA bisulfite sequencing that enables the simultaneous detection of m⁵C, Ψ, and m¹A at single-base resolution transcriptome-wide. A key discovery in RBS-seq was the observation of Ψ-dependent deletion signatures generated by a Ψ-bisulfite adduct during RT.^40,41 Fleming et al. further investigated the chemical products⁴¹ generated at Ψ-modified sites after the bisulfite reaction in RBS-seq and identified two main adducts that were shown to induce the opening of the ribose ring and cause deletion signatures during reverse transcription. N¹-methylpseudouridine (m¹Ψ) in mRNA vaccines was shown to react similarly with bisulfite to yield ribose-opening products as well.⁴²

RBS-seq still only uncovered very limited numbers of Ψ sites with weak signatures, mainly because of the conversion ratio of Ψ into Ψ-bisulfite adducts under conventional bisulfite reaction conditions. In the conventional bisulfite reaction, the protonation at the N3 position of cytosine typically requires an acidic pH, which facilitates the attack of BS to the C6 position to generate the C-BS adduct and subsequent deamination. We reasoned that the acidic condition is not optimal for the formation of Ψ-bisulfite adducts critical to the deletion signature during RT. We hypothesized that a neutral pH could not only enhance the production of Ψ-BS adducts but also inhibit C-to-U conversion to avoid reduced sequence complexity, achieving high Ψ detection sensitivity and accuracy² (Figurea).

Indeed, when testing the synthetic RNA oligo probes of AGΨGA versus AGUGA and AGCGA, the matrix-assisted laser desorption/ionization-time-of-flight (MALDI-TOF) MS measurement clearly showed that bisulfite offered an almost quantitative conversion of Ψ to Ψ-BS adducts under neutral conditions,² but the unmodified uridines were not affected after bisulfite treatment and desulphonation. No detectable C-to-U conversion was observed. Furthermore, we systematically screened all commercially available RT enzymes and self-made evolved RTs and observed that SuperScript IV reads out Ψ-BS adducts as deletion signatures with the highest deletion ratios among all known RTs.² We then constructed BID-seq libraries using the NGS platform (Figureb). With the fully modified Ψ site within a synthetic oligo containing NNΨNN, 232 out of 256 motifs in NNΨNN gave deletion ratios over 50% at the Ψ sites after bisulfite treatment, with nearly all motifs showing >25% deletion ratios. We further validated 42, 53, and two known Ψ sites in HeLa 18S, 28S, and 5.8 rRNAs, respectively, without any false positives observed. BID-seq was applied to mRNA samples from human cell lines and mouse tissues, uncovering thousands of Ψ sites with stoichiometric information at each modified site.²

Thirteen pseudouridine synthase (PUS) enzymes are encoded in the human genome, with several specific PUS enzymes reported to install Ψ modifications in human mRNA.⁴³ The quantitative feature of BID-seq enabled us to monitor Ψ stoichiometry change in PUS-depleted cells versus the control, and we found that mRNA Ψ sites could be either installed by a single specific PUS enzyme or by multiple PUS proteins.² BID-seq revealed more than 100 Ψ-modified stop codons in 12 mouse tissues and confirmed the role of Ψ in promoting stop codon readthrough in vivo.² Collectively, BID-seq set the stage for functional and mechanistic investigation of Ψ in low-abundance mRNA. In 2023, a similar approach, PRAISE, was also reported.⁴⁴ In a further optimized BID-seq protocol, we achieved even lower background deletions at unmodified uridines after bisulfite treatment and uncovered around 8500 Ψ sites starting from as little as 10 ng poly(A) RNA from mouse embryonic stem cells (mESC;⁴⁵Table S2↗).

m⁵C exists in diverse RNA species including rRNA, tRNA, mRNA, and various noncoding RNAs. The antibody-based mapping method for RNA m⁵C, such as m⁵C-RIP-seq⁴⁶ and 5-azacytidine-mediated RNA immunoprecipitation (Aza-IP),⁴⁷ could provide neither single-base resolution nor m⁵C stoichiometry information, while miCLIP⁴⁸ requires overexpression of the mutant enzyme. In recent years, BS-seq has been increasingly utilized to examine m⁵C modifications, with several commercial RNA BS conversion kits available, including the EZ RNA Methylation Kit from Zymo Research and the Methylamp RNA BS Conversion Kit from Epigentek. Using BS-seq, recent studies^{49 −52} have revealed that m⁵C modification in mRNA and its regulator proteins impacts diverse cellular functions and plays crucial roles in development and cancer. However, the exact level and stoichiometry of m⁵C on mRNA have been a subject of debate due to the absence of a sensitive, robust, and quantitative sequencing method. Discrepancies were observed when conventional BS-seq was applied to low-abundance mRNA, with some studies detecting thousands of m⁵C sites in mRNAs⁵³ while other studies discovered only a few sites.⁵⁴ More recent studies have reported only a few hundred m⁵C sites in human and mouse transcriptomes using an improved bisulfite sequencing method and a more stringent computational approach.^50,51 These inconsistent findings have raised the need to develop more sensitive and robust methods for identifying and quantifying real m⁵C sites in mRNA.

A major challenge for RNA m⁵C BS sequencing has been the high false positive rates or high background caused by incomplete C-to-U conversion due to reduced reaction temperature and reduced reaction time to avoid severe RNA degradation. The reduced temperature is also ineffective in denaturing local secondary structures of highly structured RNAs, leading to further reduced C-to-U conversion at structured regions. Mechanistically, two competing pathways exist in bisulfite conversion of RNA, with one giving the desired C-to-U conversion and the other leading to the undesired RNA degradation. The protonated N3 nitrogen under acid conditions facilitates cytosine’s reaction with BS to give the C-BS adduct, which is converted to U-BS adduct by deamination. Subsequent desulphonation of the U-BS adduct under basic conditions generates U, completing C-to-U conversion. Alternatively, the U-BS adduct may undergo spontaneous depyrimidination to cause DNA degradation. To improve C-to-U conversion efficiency and reduce RNA damage, we developed ultrafast BS sequencing (UBS-seq)³ using a new recipe with a high BS concentration (∼10 M) and a high reaction temperature (98 °C; Figure). Applying UBS-seq to highly structured rRNA as a model, we showed that the average detected fraction for the two known m⁵C sites was >95%, while no false positive site was detected when a 5% cutoff for the unconverted rate of C was used. This new approach allowed detection of m⁵C stoichiometry in highly structured RNA species and outperformed all the reported BS conditions in terms of lower background and higher sensitivity in detecting real m⁵C sites in structured rRNA.³

When UBS-seq was applied to polyA⁺-enriched RNA from HeLa and HEK293T cell lines, input mRNA as low as 10–20 ng could yield transcriptome-wide m⁵C quantification, identifying 2723 and 2404 m⁵C sites with a modification fraction ≥5%, respectively.³ The quantitative nature of UBS-seq allowed us to reveal sequence motifs of m⁵C sites in mRNA and assign NSUN2 as the main m⁵C methyltransferase that installs ∼90% m⁵C sites to HeLa mRNA. To further validate the detected m⁵C sites with low modification fractions, with NSUN2 and NSUN6 as potential mRNA m⁵C “writer” proteins, we conducted rescue experiments by transfecting the corresponding methyltransferase plasmids back to the HeLa cells with NSUN2 or NSUN6 depletion and sequenced the isolated poly(A)-tailed RNA. Indeed, we observed that the decreased m⁵C fractions in the depleted strains were mostly rescued, further confirming that these lowly modified m⁵C sites are real. In addition, our results showed that m⁵C sites deposited by NSUN2 but not by NSUN6 are enriched in 5′-UTR regions in both HeLa and HEK293T mRNA, suggesting that m⁵C modification or its binding proteins may be involved in regulating mRNA translation.³ This new method and the data sets will aid future functional investigations on RNA m⁵C (Table S3↗).

We and others previously detected N¹-methyladenosine (m¹A) at ∼0.02% m¹A/A abundance in mammalian poly(A)-tailed RNA.^55,56 Different from m⁶A and Ψ, which cannot naturally induce misincorporation signatures during RT, the methyl group at the N1 position of the m¹A base can disrupt base pairing and induce misincorporation in the presence of many commercially available RT enzymes, such as HIV RT, AMV RT, SuperScript II RT, and SuperScript IV RT. The mutation signals tend to be weak with notable RT stops observed. The TGIRT-based m¹A-MAP displays excellent performance in m¹A site detection in tRNA, but the low turnover of the TGIRT enzyme impedes the actual application of m¹A-MAP for longer RNAs.⁵⁷ We and our collaborators set up an evolution platform to evolve engineered RT enzymes for efficient readthrough of m¹A with high rates of mutation signatures. We developed a fluorescence-based RT evolution platform for the direct evolution of RT enzymes, to select enzymes that give high RT misincorporation ratios of any given type of RNA modification.¹⁰ We started with HIV RT and identified RT1306, an evolved HIV RT to give robust readthrough and high misincorporation rates at m¹A sites. Taking advantage of the evolved RT1306 and the AlkB-mediated demethylation at m¹A sites as controls, we developed m¹A-quant-seq¹⁰ (Figurea, which uncovered several hundred m¹A sites in human polyA-tailed RNA, with stoichiometric information (Table S3↗). m¹A-quant-seq also uncovered an array of mRNAs and lncRNAs containing highly modified m¹A sites (above 30% modification stoichiometry),¹⁰ such as MALAT1, mt-ND5, PRUNE, etc., which is consistent with the previous publications and suggested functional roles.⁵⁸

During our studies of m¹A in mammalian cytosolic tRNA and mitochondrial tRNAs, we noticed m¹A methylation on nascent mitochondrial RNA (mt-RNA). Mitochondrial transcription is unique in human cells, with bidirectional transcription producing a long precursor mt-RNA that contains two mt-rRNA, 13 mt-mRNA, and 22 mt-tRNA serving as junctions within the mitochondrial polycistronic RNA. In addition to m¹A methylations, other base methylations appear to also play roles in nascent mt-RNA processing. We therefore developed demethylation-assisted multiple methylation sequencing (DAMM-seq),⁴ which starts with ∼10 ng of input RNA and quantitatively maps not only m¹A but also N³-methylcytidine (m³C), N¹-methylguanosine (m¹G), and N²,N²-dimethylguanosine (m²₂G) in mitochondrial polycistronic RNA. These modifications all block the Watson–Crick base pairing interface and exist in limited sequence motifs in tRNAs; HIV RT already reads through these known motifs well and induces high misincorporation rates for quantitative stoichiometry determination.⁴

DAMM-seq utilizes the misincorporation ratios obtained at each methylated site for estimating the methylation stoichiometry, enabling the quantitative characterization of m¹A, m¹G, m³C, and m²₂G in nascent mt-RNAs.⁴ As an application example of quantitative DAMM-seq, we applied DAMM-seq to sequence mitochondrial polycistronic RNA under different cellular treatments, particularly ALKBH7 depletion and overexpression. DAMM-seq sensitively monitored the methylation level changes at methylated m¹A, m¹G, m³C, and m²₂G sites within mitochondrial polycistronic RNA and uncovered an ALKBH7 demethylation effect at m¹A and m²₂G sites within pre-tRNA regions of mt-Leu1 and mt-Ile, respectively. This study identified ALKBH7 as an RNA demethylase regulating mitochondrial RNA processing and mitochondrial activity.⁴ DAMM-seq serves as an effective approach for the quantitative investigation of multiple RNA methylations within nascent RNA species of an input RNA amount as low as 10 ng (Table S3↗). A similar method named PANDORA-seq was also reported around the same time.⁵⁹

Among the expanding list of functionally relevant post-transcriptional RNA modifications, methylation of the 2′-OH of an RNA ribose (2′-O-methylation, N_m) is both abundant and unique. 2′-O-methylation is one of the most common RNA modifications present in many cellular RNAs, such as tRNA, rRNA, and small nuclear RNA (snRNA). It also exists in mammalian mRNA and is about 10-fold less abundant compared to m⁶A. Unlike base modifications, which are specific for one of the four RNA bases, N_m methylation occurs at the 2′-OH position of each ribonucleoside. N_m modifications can modulate the secondary structure of rRNA and tRNA and affect mRNA stabilization and translation.^60,61 While functional roles of N_m modifications on abundant tRNAs and rRNAs have been well documented, studies of their impacts on lower-abundance RNAs, such as mRNAs and long noncoding RNAs (lncRNAs), have been hampered by the lack of effective antibodies and robust sequencing methods that map N_m at base resolution.

Previous base-resolution approaches using NGS have been successful at identifying 2′-O-methylation in highly abundant RNA.^{62 −66} However, most of these methods cannot provide stoichiometric information at the modified sites. RiboMethSeq⁶⁵ is the only quantitative method developed thus far, but very deep sequencing depth is required for confident N_m detection in low-abundance RNAs and at low-stoichiometry N_m sites. In addition, most current N_m mapping methods mainly rely on RT truncation signatures induced by chemical-assisted cleavage, which may include substantial false positives arising from RT stops at structured regions in RNA or unmodified RNA ends from fragmentation biases or RNA processing. Compared with RNA truncation signatures, modification-dependent misincorporations as readouts are generally more reliable in modification detection and quantification. However, since N_m modification occurs on the ribose and does not directly modulate the Watson–Crick base pairing interface during reverse transcription, all known reverse transcriptases (RT) are unable to induce misincorporations, which has precluded the use of mutation-based analysis methods for N_m mapping. Leveraging the fluorescence-based RT evolution platform to generate RTs that are mutagenic at specific RNA modification of interest,⁸ we speculated that evolved RTs might be selected to generate misincorporations opposite of N_m-modified bases.

Indeed, a new HIV RT variant, RT-41B4, was selected to yield notable misincorporation signatures at N_m-modified sites. In the presence of adjusted dATP/dNTP ratios during RT, we developed N_m-Mut-seq (an RT-mutation-based N_m mapping method)¹¹ that maps C_m, G_m, and A_m methylations at single-base resolution with stoichiometry information (Figureb). To validate N_m-Mut-seq, we mapped almost all known C_m, G_m, and A_m sites on human rRNAs with high misincorporation rates and without observing any false positives at unmodified sites or other rRNA modifications.¹¹ Applying N_m-Mut-seq to HepG2 cellular mRNA, we uncovered around 1000 C_m, G_m, and A_m sites, a subset of which were validated by multiple orthogonal methods,¹¹ further confirming the fidelity of the method. Fibrillarin (FBL) was shown to be a main N_m “writer” protein for HepG2 mRNA, and the methylation is guided by snoRNAs that interact with the target mRNAs. Note that the current N_m-Mut-seq still requires ∼200–800 ng of input polyA⁺ RNA (Table S3↗); otherwise, the high PCR cycle numbers lead to dramatic PCR duplicates in the obtained NGS data. N_m-Mut-seq, therefore, provides an effective method to not only detect N_m at base resolution in low abundant RNAs such as mRNA and lncRNA but also measure the stoichiometry of the modified sites transcriptome-wide.

m⁷G methylation is a well-known modification in the cap structure of the 5′ end of mammalian mRNA, which impacts diverse biological processes such as mRNA stability, splicing, nuclear export, and translation. Highly modified internal m⁷G sites were also found at cytoplasmic tRNA and 18S rRNA.² However, whether m⁷G methylation could be deposited internally to mammalian mRNA remained unknown for a long time. In 2019, we and another group demonstrated the existence of internal m⁷G modification in mammalian mRNA and miRNA.^12,67 We reported the antibody-based m⁷G-MeRIP-seq and a single-base resolution m⁷G sequencing method m⁷G-seq¹² (Table S3↗).

Using m⁷G-MeRIP-seq, we identified METTL1 as the main “writer” protein that installs mRNA internal m⁷G. m⁷G is positively charged and could be converted into a reduced state in the presence of sodium borohydride; the subsequent heating under acidic conditions leads to depurination at the reduced m⁷G site, yielding an abasic site (AP site).¹² After screening different commercially available RT enzymes, we found that the abasic site generated at the m⁷G site could be read as RT misincorporation signals using HIV RT; the abasic sites could be further captured and enriched by biotin-tagged hydrazine to yield a higher misincorporation rate. Overall, m⁷G-seq includes three libraries for each sequencing library: zero mutations at internal m⁷G sites in the “Input” library, a moderate mutation rate induced at abasic sites in “Before pulldown” libraries, and a high mutation rate induced at biotin-enriched abasic sites in “Pulldown” libraries.¹² The base-resolution m⁷G-seq confirmed the presence of mRNA internal m⁷G methylations in cancer cell lines. Note that a recent study revealed that quaking proteins (QKIs) bind to internal m⁷G-modified mRNAs with GA-rich motifs, as the first reported internal m⁷G “reader” proteins,⁶⁸ suggesting functional roles of mRNA internal m⁷G methylations in mammals.

Besides the role of METTL1 as the 'writer' protein for mRNA internal m⁷G installation, METTL1 is well-known to mediate m⁷G methylation at position 46 of cytoplasmic tRNA, and this methylation has been shown to link with cancer progression and tumorigenesis,^{69 −71} requiring a quantitative method for monitoring tRNA m⁷G dynamics in pathological processes. Based on the chemical principle in m⁷G-seq, we further optimized it to be more quantitative (Figure), or m⁷G-quant-seq,¹³ without the need for the biotin pulldown enrichment. In m⁷G-quant-seq, we performed the depurination at a more acidic pH (∼2.9) under mild heating. The internal m⁷G sites could be almost completely converted to abasic sites when extending the incubation time of both KBH₄ reduction and depurination.¹³ Note that the current m⁷G-quant-seq requires ∼200 ng of cellular small RNA as input, to avoid high PCR cycle number and the consequent PCR duplicates in NGS¹³ (Table S3↗).

Li-Sheng Zhang earned his B.S. in Materials Chemistry from Peking University and received his Ph.D. in Chemistry at The University of Chicago where he worked with Professor Chuan He. He worked as a postdoctoral scholar at Professor He’s laboratory before starting his independent research in Hong Kong. Currently, he is a tenure-track Assistant Professor of chemistry and life science at The Hong Kong University of Science and Technology (HKUST). His research interests focus on quantitative sequencing technology and biological functions of RNA modifications.

Qing Dai earned his B.S. degree in chemistry from Central China Normal University and completed his Ph.D. in organic chemistry at Nankai University. He then carried out postdoctoral work in nucleic acid chemistry at the University of Chicago. Currently, he is a Research Professor at the University of Chicago. His research interests include developing new sequencing methods for DNA and RNA modifications and their applications.

Chuan He received his Ph.D. from Massachusetts Institute of Technology (MIT), where he worked with Professor Stephen J. Lippard. He did his postdoctoral work with Professor Gregory L. Verdine at Harvard University. Currently, He is a John T. Wilson Distinguished Service Professor at The University of Chicago and an investigator at Howard Hughes Medical Institute.