Transcriptome-wide Mapping of Internal N⁷-methylguanosine Methylome in Mammalian Messenger RNA

Quantitative Detection of Internal N⁷-methylguanosine Sites within Mammalian mRNA

To assess the possible presence of m⁷G in mRNA internally, we employed a liquid chromatography-tandem mass (LC-MS/MS) method to quantify the modification level of thoroughly purified poly(A)-tailed RNA (with two rounds of poly(A) enrichment followed by rRNA depletion). Aiming to monitor m⁷G signals from both internal and 5’-cap regions, we applied a 3-step method for RNA digestion (Chu et al., 2018), which includes the endonuclease cleavage (nuclease S1) followed by a treatment with phosphodiesterase I to completely expose m⁷G signal from the cap structure (Figure S1A). We observed an m⁷G/G level of approximately 0.45% in HEK293T cells (“polyA⁺ RNA” on the left in Figure 1B), including m⁷G sites at both the 5’-cap and internal regions. We then removed the cap using tobacco decapping enzyme to specifically quantify the internal m⁷G level (see STAR Methods). The m⁷G level decreased dramatically after the first round of decapping treatment but remained at a stable level around 0.03% consistently, which did not diminish after two more rounds of cap removal (Figure 1B). We also performed RNA digestion in a two-step manner without phosphodiesterase I digestion, and found the m⁷G level detected was comparable to that after two to three rounds of decapping treatments (“polyA⁺ RNA” on the right in Figure 1B), which suggests the digestion of mRNA with only endonuclease cleavage (nuclease S1) would not affect cap structure and can ensure reliable exposure of internal m⁷G level in mRNA. We measured m⁷G/G levels ranging from approximately 0.02% – 0.05% across several human and mouse cell lines (Figure 1C), which represent roughly 5–10% of the m⁶A/A ratios in the same cells. Thus, m⁷G is present internally within mammalian mRNA.

Transcriptome-wide Mapping of Internal N⁷-methylguanosine with m⁷G-MeRIP-seq

To study the distribution of internal m⁷G methylome, we performed m⁷G-MeRIP-seq to map the transcriptome-wide distributions of internal m⁷G sites. We first conducted dot blot to characterize that the selected antibody (MBL) binds to m⁷G specifically over all known methylated G (including m¹G, m²G, m²₂G, and m⁶G) or unmethylated G (Figure S1B). Then we briefly immunoprecipitated methylated RNA after cap removal (with two rounds of decapping treatments) using the anti-m⁷G specific antibody (MBL), followed by high-throughput sequencing (Figure S1C). All the sequencing samples for mRNA were prepared from total RNA with two rounds of polyA⁺ enrichment where tRNAs were depleted and rRNAs were eliminated to below 0.1% of the original abundance (Figure S1D). Based on m⁷G-MeRIP-seq, we found that m⁷G was enriched around 12 folds in the anti-m⁷G antibody-bound fraction, starting with cap-depleted and fragmented poly(A)-tailed RNA from HEK293T cells. In contrast, other known mRNA modifications including m⁶A, Nm (Gm and Am) or m¹A, were not enriched (Figure 1D). By applying m⁷G-MeRIP-seq to poly(A)-tailed RNA (cap- and rRNA-depleted) from human and mouse cells in two replicates, we identified 7,469, 7,911, 9,243, 6,184 and 6,743 potential internal m⁷G peaks in HeLa, HEK293T, HepG2, MEF and mES cells, respectively, in adequately expressed transcripts shared in replicates. We detected 3,823 m⁷G peaks that overlapped in the three human cell lines, indicating a relatively high confidence of antibody enrichment, and 2,217 m⁷G peaks overlapped in the two mouse cell lines (Figure 1E). About 50% of the methylated transcripts are methylated once, and each of them carries around 1.9 peaks on average (Figure 1F). A statistically significant GA-rich motif was found in the MeRIP-seq peak regions (Figure 1G).

We investigated the distribution of methylated genes (1≤FPKM≤50) in ten averagely divided expression bins. When calculating the percentage ratio of methylated genes versus all genes in each expression decile, the methylated transcripts exhibit a progressively larger fraction as gene expression level increases in both human and mouse cell lines (Figure 2A); we thus excluded the possibility that the increased fractions of the methylated genes are likely due to increased detection of highly expressed genes (Figure S2A). We detected internal m⁷G peaks in all three transcript segments of 5′ untranslated regions (UTR), coding sequences (CDS), and 3′ UTRs. Similar distribution patterns of internal m⁷G peaks were observed among all three human cell lines and two mouse cell lines (Figure 2B and Figure S2B), with the same case for the distributions of the genes bearing m⁷G (Figure S2C). We also plotted the fractions of genes with m⁷G peaks in each of the segments relative to expression level and found that the transcripts of relatively high expressed genes were more likely to be methylated at CDS and 3’ UTR segments than 5’ UTR (Figure S2D). Meta-gene profiles of internal m⁷G peaks along mRNA demonstrated a major accumulation area located at 3’ UTR and a minor one within 5’ UTR (Figure 2C). Although m⁷G is enriched at the 3’ UTR, there is nearly no overlap between m⁶A and m⁷G methylated peaks in HeLa cell line (Figure S2E). Furthermore, we investigated the peak conservation between human and mouse, and found that 2,557 identified m⁷G peaks in mESC could be mapped faithfully to adequately expressed orthologous in human HepG2 cells; 1,016 of them had an m⁷G peak at the orthologous positions, representing 39.7% overall conservation. The CDS regions exhibited the highest degree of conservation with around 41–43% conserved peaks, compared to other gene segments with nearly 15–20% conservation (Figure 2D, and Figure S2F). Gene ontology (GO) analysis revealed that the methylated genes are enriched under GO terms related to RNA splicing, mRNA processing, cell-cell adhesion and cell division (Figure 2E).

Design and Development of m⁷G-seq at Base Resolution

To achieve mapping of m⁷G at base-resolution with an orthogonal approach, as has been done for other modifications (Li et al., 2017), we developed m⁷G-seq by taking advantage of the unique chemical reactivity of m⁷G in a reduction-induced depurination reaction (Figure 3A). The positive charge on the five-membered ring makes m⁷G (3a–1 in Figure 3A) especially susceptible to NaBH₄-mediated reduction, which eliminates the aromaticity of the five-membered ring attached to the ribose without affecting unmodified G (Wintermeyer and Zachau, 1970). Subsequent heating (55 °C) at pH 4.5 induces depurination of the reduced m⁷G alone (3a–2 in Figure 3A) and generates an RNA abasic site (3a–3), which can be further captured by biotin-ligated hydrazide to yield a biotinylated RNA in a one-pot reaction (3a–4). After pull-down and reverse transcription using HIV reverse transcriptase, the biotinylated sites would be mutated to predominantly T as well as other bases (Küpfer and Leumann, 2005; Küpfer and Leumann, 2007), and thus we can identify the internal m⁷G sites based on these mutations at single-base resolution.

We validated the m⁷G-seq assay including the chemical reactions as well as the pull-down enrichment step using model RNAs (Figure S3A). 125-mer RNA probes with a single G (unmodified or m⁷G modified) incorporated site-specifically were synthesized and subjected to dephosphorylation, reduction and biotinylation consecutively. Dephosphorylation was necessary to avoid potential side reactions of free phosphate groups (at the 5’ end of the RNA) with biotin-ligated reagents. After biotin pull-down and reverse transcription, the oligos with m⁷G modification presented a highly enriched signal compared with the unmodified ones at full length (Figure S3B), which indicated that the reverse transcriptase could read through the biotinylated abasic sites with misincorporation at those sites as demonstrated by Sanger sequencing of the further PCR-amplified products (Figure S3C).

We then optimized the chemical assay to be compatible with library construction and sequencing (Figure 3B). mRNA fragments (around 60–80 mer) were subject to cap removal and PNK end repair to leave the hydroxyl group at the 3’ end. Then dephosphorylation further ensured the depletion of phosphate groups at both ends. After the 3’ adapter installation to fragmented poly(A)-tailed RNA, the ligated mRNA fragments underwent a two-step chemical treatment before further enrichment by the biotin tag. We then performed the reverse transcription followed by cDNA 3’-adaptor ligation and PCR amplification (see STAR Methods). In this way, as m⁷G sites were converted to biotinylated abasic sites, we can observe misincorporation at the specific sites that were further enriched in the pull-down samples compared with input samples and those before biotin enrichment. To evaluate this method before applying it to mRNA, we designed a 60-mer RNA oligo sequence which is similar to mRNA fragment size and bears one single fully methylated m⁷G in the AAGAU sequence. After applying the aforementioned protocol, we observed substantial mutations at the modified site after generation of the abasic site (~28% overall mutation rate), and higher mutation rates after further pull-down enrichment (~51% overall mutation rate), whereas the adjacent bases did not exhibit any misincorporation above background signals (Figure S3D). Note that the chemical reduction and depurination reactions are not stoichiometric with fractions of m⁷G escaping reduction and depurination. The pull-down step further enriches modified abasic sites, leading to increased mutation rates.

While the biotin-based enrichment allows confident assignments of the modified m⁷G sites, the modification fraction could in principle be measured using the misincorporation rate before enrichment and correlate with a calibration curve generated using RNA probes that bear m⁷G at different modification fractions and undergo the same chemical treatment and high-throughput sequencing. A series of probes with the same AAGAU sequence but different m⁷G/G ratios was prepared from in vitro transcription and then subjected to the base-resolution m⁷G-seq protocol (Figure 3B). A calibration curve was generated that can normalize the misincorporation rate (before enrichment) to the methylation status, which can then be used to estimate internal m⁷G methylation fraction in biological RNA samples (Figure S3E).

Single-Nucleotide Resolution m⁷G Maps in Human rRNA and tRNAs

We first applied m⁷G-seq to rRNA and tRNAs isolated from human cells and examined the known m⁷G₁₆₃₉ site in human rRNA (Haag et al., 2015). We observed highly increased mutation only at the known m⁷G site. With reduction and biotinylation, we could obtain a misincorporation pattern of 18.0% G→T, 9.9% G→C, and 2.8% G→A (~31% overall mutation rate) at m⁷G₁₆₃₉ site in HeLa rRNA; after pull-down enrichment of biotinylated fragments, we detected a distribution of misincorporation with 38.2% G→T, 13.0% G→C and 3.7% G→A (~55% overall mutation rate) (Figure S4A), implying a predominant G to T or C mutation pattern induced by the chemical treatment followed by RT using HIV reverse transcriptase. The overall mutation rate is around 30% at the m⁷G₁₆₃₉ site before enrichment across all three human cell lines (Figure S4A), suggesting nearly 100% modification at this site (Figure S3E).

Our method also identified m⁷G sites in 22 human tRNAs at the position 46 (Figure 4A, Figure S4B, and S4C). These tRNAs, with different sequence contexts around m⁷G₄₆, varied in the misincorporation rate in a range of 5–31% before biotin enrichment, suggesting an estimated high methylation status ranging from 60%–100% (Figure 4A, 4B, and Figure S3E). Our observation of m⁷G modification in tRNA is also consistent with a recent report based on a cleavage assay (Lin et al., 2018).

As tRNAs tend to be heavily modified, we also characterized various mutation patterns of other G modifications to avoid potential mis-assignments (Figure 4C). The dimethylated m²₂G₂₆, which blocks Watson-Crick base pairing, showed high misincorporation rates throughout all three samples with a majority of G→T mutation up to around 90% as expected; however, the monomethylated m²G₁₀ showed relatively low mutation rates below 5%. When we examined m¹G, the misincorporation patterns differed based on the modification locations. While both m¹G₉ and m¹G₃₇ displayed the overall mutation rate around 40–60%, m¹G₉ tended to give higher G→T mutation but m¹G₃₇ caused more G→C mutation in our assay. Importantly, none of these G methylations showed an increased misincorporation after reduction and depurination treatment compared with m⁷G; we observed nearly no mutation in input samples, but high mutation in samples before enrichment and after pull-down only for m⁷G. Therefore, we compare mutation rates of Gs before and after our chemical manipulation as a criterion to identify m⁷G sites in our later analysis to exclude other guanosine modifications. As O⁶-methylguanosine (m⁶G) was also observed in both DNA and RNA as a highly mutagenic base (Hudson et al., 2015), we constructed a 60-mer probe with one O⁶-methylguanosine (m⁶G) in the middle by in vitro transcription. We subjected this probe to m⁷G-seq treatment and high-throughput sequencing. The m⁶G sited exhibited a mutation pattern with mainly G→A mutation (~7.5%) in both input and before pulldown samples but could not be enriched with biotin pulldown (Figure S4D). Therefore, we can further exclude m⁶G in m⁷G-seq.

Because mammalian m⁷G₄₆ sites in tRNAs are known to be installed by the METTL1/WDR4 complex (Lin et al., 2018; Alexandrov et al., 2002; Cartlidge et al., 2005). We applied m⁷G-seq to METTL1 knockdown cells with both transient knockdown (siRNA) in HeLa cells and stable knockdown (shRNA) in HepG2 cells to further validate the applicability of our method. Both siMETTL1 and shMETTL1 cells showed obvious decreases in misincorporation rate at all m⁷G₄₆ sites (Figure 4D and 4E), showing reduced m⁷G modification fraction at these sites with METTL1 knockdown. Interestingly, the extent of reduction of m⁷G is not universal with METTL1 knockdown: while tRNA^Ala(UGC) and tRNA^Thr(UGU) showed nearly 30% decreases in the methylation level, tRNA^Pro(HGG) and tRNA^Val(AAC) only showed a few percent m⁷G decreases. Either METTL1/WDR4 has preferential tRNA targets or there could be additional m⁷G methyltransferase(s) that could compensate m⁷G methylation by METTL1, which should be investigated in the future.

Single-Nucleotide Resolution mRNA m⁷G Methylomes

We next applied m⁷G-seq without pull-down enrichment to mRNAs purified from HeLa and HepG2 cells, as an attempt to identify frequently methylated sites. We used 3% mutation rate above background input mutation as a cutoff and computationally identified over 90 sites in mRNA without the enrichment step; these sites potentially possess estimated methylation levels around or above 20% based on estimations using our calibration curve (Figure S3E) and comparison to tRNAs. We further identified ~50 sites manually with modifications estimated to be around or above 50% (Table S2). These sites showed misincorporation rates and possess m⁷G methylation levels approaching those in tRNAs (Figure 4 and Figure 5A). Representative mRNA internal m⁷G sites conserved in HeLa and HepG2 cells are shown in Figure 5A. The AG(m⁷G)A motif shown in CNOT2 (5’ UTR) and EIF2S2 (CDS) resembles the same motif for rRNA m⁷G₁₆₃₉. Our m⁷G-MeRIP-seq found dominant GA-enriched motif, which was also discovered in a number of purine-enriched region, with GG(m⁷G)AA in 3’ UTR of GPR107 and UA(m⁷G)AA in 3’ UTR of RTN4 as examples. We also found motifs such as AUCG(m⁷G)A and UU(m⁷G)AU with high modification fractions conserved in both cell lines, suggesting presence of multiple methyltransferases that specifically install these sites.

Considering that m⁷G may play functions at specific cellular locations (i.e. chromatin associated RNA or nascent RNA), we also wanted to identify lowly to moderately modified m⁷G sites from the whole-cell mRNA and compare with m⁷G-MeRIP-seq results. We then performed enrichment of the biotin-tagged mRNA fragments obtained from m⁷G reduction and depurination, followed by reverse transcription using HIV RT enzyme. We used 8.0% mutation rate increase over the background input control as our cutoff because the modified m⁷G sites show higher mutation rate after enrichment. We identified around 3500 internal mRNA m⁷G sites in HeLa cells from overlap of two replicates. We then overlaid HeLa m⁷G sites to those shared in HepG2 replicates and identified 801 overlapped ones (Table S3). Nearly 70% of these 801 single sites overlap with targets revealed from m⁷G-MeRIP-seq results, showing consistency between two orthogonal approaches. It should be noted that m⁷G-MeRIP-seq can enrich m⁷G sites ranging from low modification fraction to high modification fraction, as well as potential non-specific binding. The base resolution m⁷G-seq utilizes mutation over background input to detect m⁷G with much higher accuracy, but it may miss most lowly modified sites due to inadequate sequencing depth and our stringent mutation cutoff criterion. HIV reverse transcriptase itself is known to generate baseline mutation rate, thus applying a stringent mutation rate cutoff is important to ensure reliable assignments of m⁷G sites despite potential loss of lowly modified sites.

We reasoned that some of the m⁷G sites identified by m⁷G-seq without overlap with m⁷G-MeRIP-seq peaks may still be genuine m⁷G sites. They may exist in structured RNA which hindered the antibody-based pulldown. Nevertheless, to obtain more accurate assignments we decided to assign m⁷G sites not only shared in two human cell lines but also overlapped with m⁷G-MeRIP-seq peaks as confident sites (Figure S5A). Gene ontology (GO) analysis uncovered the similar enriched GO terms of cell-cell adhesion and RNA splicing from these 801 sites, which correlates well with the m⁷G-MeRIP-seq data (Figure 5B).

Around half of the more confident m⁷G sites or m⁷G methylated genes exist in CDS, while 5’ UTR showed a little increase in distribution rate, but still lower than that of 3’UTR (Figure 5C and Figure S5B). On average we observed mostly one methylated peak per gene and this is true to all three transcript segments (Figure S5C). The general misincorporation level ranged from 8%−35% after pull-down enrichment (Figure S5D). We then analyzed the sequence motifs of the base-resolution mRNA internal m⁷G sites. The main motifs were found to be G(m⁷G)A, A(m⁷G)A, C(m⁷G)C, G(m⁷G)G, C(m⁷G)U and G(m⁷G)C, which showed wider variations than those in tRNAs (Figure 5D). We picked up seven representative motifs and plotted them along 5’ UTR, CDS, and 3’ UTR. The G(m⁷G)G motif appears to be more likely at 3’ UTR, while the C(m⁷G)C motif accumulates more at 5’ UTRs and barely appears at 3’ UTR (Figure 5E).

A Subset of Internal m⁷G sites are Installed by METTL1

The similarity of some of the m⁷G motifs observed within mRNA to the known AG(m⁷G)H motifs (Okamoto et al., 2004; Lin et al., 2018) identified in tRNAs suggested that the METTL1-WDR4 heterodimer complex, the known tRNA m⁷G methyltransferase (Alexandrov et al., 2002), might introduce some of these internal m⁷G modifications in mammalian mRNAs. To test this hypothesis, we knocked down METTL1 in HeLa cells and quantified the corresponding internal m⁷G level within cap- and rRNA-depleted mRNA. We observed a ~25% reduction of the m⁷G/G ratio in siMETTL1 cells compared to controls (Figure S6A). To map changes in m⁷G peaks, we applied m⁷G-MeRIP-seq to mRNAs isolated from siControl vs. siMETTL1 HeLa cells as well as shControl vs. shMETTL1 HepG2 cells. The knockdown efficiency in both mRNA and protein levels was validated with RT-qPCR and western blotting, respectively (Figure S6B). We uncovered 1,563 altered peaks in transient knockdown cells and 3,939 in stable knockdown ones; in each case, around 74–77% of peaks became hypo-methylated upon METTL1 knockdown (fold change (FC) ≥ 1) (Figure S6C). In these 1,227 hypo-methylated peaks in HeLa cells, a majority were enriched in CDS and 3’ UTR, while METTL1 stable knockdown HepG2 cells showed a high enrichment of these potential METTL1 target m⁷G peaks in CDS (Figure 6A and Figure S6D). The metagene profile of these hypo-regulated targets (affected by METTL1 knockdown) also demonstrated a broad cumulation near stop codon both in CDS and 3’ UTR (Figure S6E). METTL1-related internal m⁷G peaks identified in two different cell lines upon knockdown shared similar GA-enrich motifs (Figure 6B) as the top motifs that were also found by m⁷G-MeRIP-seq (Figure 1G). The METTL1 knockdown resulted in ~54% global reduction of these 1,227 m⁷G peaks in HeLa cells, and ~61% reduction in HepG2 cells (Figure 6C), with the reduction occurring consistently in all three transcript segments (Figure S6F and S6G).

To further investigate whether the METTL1-WDR4 complex installs m⁷G onto mRNA, we performed photoactivatable ribonucleoside-enhanced crosslinking and immunoprecipitation (PAR-CLIP) of both proteins. We found that both proteins bind mRNA, with over 35% and around 4% of METTL1-dependent m⁷G-MeRIP-seq peaks overlapping with WDR4 or METTL1 PAR-CLIP peaks, respectively, while the whole set of HeLa m⁷G-MeRIP-seq peaks (Figure 1E) exhibits much lower overlap with these CLIP peaks (Figure S6H), implying the existence of other potential mRNA m⁷G methyltransferases. The observed hypo-methylation upon METTL1 knockdown and the correlation between METTL1 and WDR4 targets suggest that METTL1, together with WDR4, act as a “writer” complex of m⁷G within mRNA. WDR4 might facilitate RNA binding, with a stronger or direct interaction with RNA as suggested by the larger number of overlapped targets between its PAR-CLIP sites with the METTL1-dependent m⁷G-MeRIP-seq peaks, whereas METTL1 could serve as the catalytic component for methylation.

Gene Ontology analysis of the 1,227 hypo-methylated transcripts (upon siMETTL1 knockdown) uncovered a potential correlation of m⁷G methylation with RNA splicing as well as cell adhesion (Figure 6D). We then further conducted biochemistry to treat small RNA fraction (<200 nt; mostly tRNAs) isolated from shMETTL1 and mRNA isolated from siMETTL1 cells, respectively, using the recombinant protein complex of METTL1/WDR4 prepared with the transfection of both FLAG-METTL1 and FLAG-WDR4 plasmids in human HEK293T cells. Substantially elevated m⁷G/G levels were observed after incubation in the presence of SAM (Figure 6E). The m⁷G/G level was recovered to 0.8 % when small RNAs from shMETTL1 HepG2 cell were treated with the complex (Figure 6E), confirming tRNA methylation activity of METTL1/WDR4 complex. The mRNA m⁷G/G ratio went up to ~0.22% compared with the original 0.03% in the input mRNA, indicating in vitro mRNA methylation activity of the same complex. When we treated the same METTL1/WDR4 complex with 12–17 mer RNA probes bearing the consensus motif for m⁷G we failed to observe any activity (see STAR Methods). We then constructed a 144-mer RNA oligo probe that mimic the sequence of a m⁷G site in mRNA UTRs, and subjected it to the same biochemistry assay. A m⁷G/G level of 0.40% was obtained, which afforded 0.26 equivalent of m⁷G in each probe on average (Figure S6I). These results confirmed the m⁷G methylation activity of METTL1/WDR4 complex and suggested that METTL1/WDR4 does not work on short RNA probes and may recognize tRNA-like structures in order to mediate m⁷G methylation on mRNA.

METTL1 Knockdown Decreases Translation Efficiency of m⁷G-Modified Transcripts

To investigate potential functions of METTL1-mediated internal m⁷G mRNA methylation, we first assessed mRNA-ribosome association by performing ribosome profiling in HeLa cells transfected by either siMETTL1 or control siRNAs. We sequenced both the input and the ribosome-bound mRNA fragments, and analyzed translation efficiency (ribosome-bound mRNA normalized by the whole cell mRNA level) of five groups of transcripts: 1) 1,437 non-target mRNAs containing m⁷G-MeRIP-seq peaks that were unaffected by METTL1 knockdown; 2) 972 METTL1 potential target (PT) mRNAs including m⁷G-MeRIP-seq peaks that were hypomethylated upon METTL1 knockdown; 3) 397 PT mRNAs that overlap with either WDR4 or METTL1 PAR-CLIP peaks; and 4) 101 high confidence PT mRNAs that overlap with m⁷G base-resolution sites and METTL1/WDR4 PAR-CLIP peaks. We found that METTL1 depletion decreased translation efficiency of all four PT groups, suggesting that loss of METTL1 and thus reduction of internal m⁷G may suppress translation of some of METTL1 target transcripts relative to non-targets (Figure 7A and 7B, Figure S7A).

We next performed polysome profiling to examine the distribution of translating METTL1 target mRNAs in METTL1-knockdown cells versus control cells. Three groups of mRNA-protein particles (mRNPs) were separated: non-ribosome, 40S-80S, and polysome. Western blotting of each fraction showed that the cytoplasmic METTL1 and WDR4 were distributed mainly in the non-ribosome and 40S mRNPs portions (Figure S7B). We quantified changes in the distribution of specific mRNAs among non-ribosome fractions and polysome fractions using RT-qPCR. Upon METTL1 knockdown, the mRNA levels within polysome fractions decreased for the METTL1 targets (DDAH1, SLC20A1, SNAP23, and VCL), but not for a non-target control (ACTB) (Figure 7C and Figure S7C). These profiling results provided support that the METTL1-dependent internal m⁷G methylation could promote translation of the corresponding transcripts.

To investigate whether METTL1 knockdown influences the translation of target mRNAs via altered tRNA methylation, we examined the enrichment folds of specific tRNA codons related to m⁷G-methylated tRNAs within the METTL1 targets (group 2) vs. randomly selected non-targets (Figure S7D). We observed no enrichment of these tRNA codons in METTL1 targets vs. non-targets. These data suggest that METTL1 knockdown could specifically reduce the translation of METTL1 target mRNAs independently of effects on global translation and tRNA methylation.

Frequently Modified m⁷G Sites in Human tRNA and mRNA

Applying the m⁷G-seq approach (without enrichment) to tRNAs, we successfully characterized the base-resolution human cytosolic tRNA m⁷G₄₆ methylation, which clearly distinguished m⁷G from other guanosine modifications. These tRNA m⁷G sites all exist with high modification fractions. We also performed m⁷G-seq in METTL1 knockdown cells, and confirmed the tRNA m⁷G methylation function of METTL1. The ability to monitor the m⁷G methylation status by m⁷G-seq may aid future functional studies of tRNA and rRNA m⁷G methylation and monitor potential dynamics of these modified sites during cellular processes.

Our method also uncovered over 90 internal m⁷G sites within human mRNA that may exhibit around to above 20% modification fraction. Approximately 50 of them can be considered as highly methylated internal m⁷G sites in human mRNA. The identification of these highly methylated sites, conserved in two different human cell lines, will stimulate further functional characterizations of these modifications and the target mRNAs. We found internal m⁷G with high modification levels displays GA-enriched motifs near start codon and stop codon. Potential effects on translation regulation could be studied in the future with proper assays (Karijolich et al., 2011).

We further conducted biotin pull-down enrichment to unveil additional internal m⁷G sites and to verify the sites discovered before enrichment. After pull-down enrichment, the misincorporation rate difference (above the input baseline mutation rate potentially introduced by HIV RT) was elevated to 8–35%. After applying stringent cutoff, we were able to identify more than 800 sites shared by two human cell lines, from which we performed statistical analysis for features of these confident internal m⁷G sites. These sites exhibit additional sequence motifs beyond the known GA- or GG-containing motifs found in tRNA and uncovered for mRNA using m⁷G-MeRIP-seq, suggesting potential presence of multiple m⁷G methyltransferases that may install mRNA m⁷G.

METTL1-WDR4 as a m⁷G Methyltransferase Complex for mRNA

We found that the known tRNA m⁷G modification complex, METTL1-WDR4 (Leulliot et al., 2008; Alexandrov et al., 2002; Lin et al., 2018), while participating in tRNA methylation as also verified in our knockdown results, also installs a subset of internal m⁷G sites in human mRNA. The simultaneous methylation of the position 46 at the extra loop in human tRNAs and human mRNA by the METTL1-WDR4 complex may not be surprising. Other mRNA modification machineries such as eukaryotic ψ synthases PUS1, PUS7 and TRUB1 have been shown to work on both tRNA and mRNA (Carlile et al., 2014; Li et al., 2015; Lovejoy et al., 2014; Safra et al., 2017; Schwartz et al., 2014). NSun family m⁵C methyltransferases bind quite extensively to mRNA as well (Hussain et al., 2013; Yang et al., 2017; Khoddami et al., 2013). Here, we show that the heterodimer METTL1-WDR4 is responsible for a portion of internal mRNA m⁷G methylation. Besides our observation that the knockdown of METTL1 led to reduced m⁷G methylation within mRNA, the METTL1-WDR4 PAR-CLIP results also revealed the binding of this complex to mRNA and a substantial overlap of the binding sites with m⁷G-MeRIP-seq sites in the same cell line (Figure S6H). WDR4 appears to possess more binding sites in mRNA and may contribute to stabilizing the heterodimer complex and RNA binding, whereas METTL1 may mediate m⁷G methylation. In vitro methylation assay using tRNA and mRNA purified from METTL1-depleted cells revealed good m⁷G methylation activity of METTL1-WDR4 towards both tRNA and mRNA. This complex also mediates m⁷G methylation of a long synthetic RNA probe that mimics internal mRNA m⁷G sites; however, the complex failed to react with several RNA probes ranging from 12-mer to 17-mer, suggesting that this methylase complex may require tRNA-like structures with a stem loop of a minimum length to facilitate substrate recognition. Please note that our m⁷G-seq identified diverse sequence motifs around internal m⁷G sites in mRNA (Figure 5D), suggesting involvements of additional methyltransferase(s) for m⁷G installation.

Functional roles of m⁷G within mRNA could be diverse. We found that the knockdown of METTL1 decreases the translation efficiency of transcripts that contain METTL1-affected hypomethylated m⁷G sites. Because of the unique nature of positive charge, internal m⁷G could reorganize local RNA structure, modulate protein-RNA interactions, and be better recognized through aromatic stacking. It could be recognized by different reader proteins and modulate local RNA structure to affect protein-RNA interactions, which might lead to regulation of mRNA translation or other processing and functions of mRNA, such as nuclear pre-mRNA processing, export, and stability. Perhaps different methyltransferases install m⁷G to distinct mRNA sequences and/or secondary structural motifs to further modulate the fate of the target mRNA. The identification of more frequently m⁷G methylated mRNAs and the availability of the new m⁷G-seq method should allow future functional studies and mechanistic investigations.

Limitations

The current m⁷G-seq protocol can still be improved in the future. We noticed that the mild chemical reactions for selective m⁷G reduction and depurination could not achieve quantitative yields, with a portion of m⁷G sites converted to abasic sites and then a portion of these abasic sites labeled with biotin hydrazide. These non-stoichiometric labeling chemical reactions explain higher mutation rates after enrichment of biotin-labeled RNA than those before enrichment. This limitation will affect accurate measurements of the modification levels of individual m⁷G sites. We have performed reactions on calibration and real biological samples under consistent conditions that should normalize variations introduced by the depurination and labeling chemistry. However, further improvement of the reaction efficiency and robustness as well as calibration consistency would further improve the quantitative feature of m⁷G-seq.

In summary, we discovered and mapped internal m⁷G sites within mammalian mRNA. Our single-nucleotide-resolution m⁷G sequencing technology allows precise mapping of m⁷G transcriptome-wide. This method could be widely used to detect m⁷G not only in tRNA and rRNA, but also within mRNA and various non-coding RNAs. We provide the first base-resolution m⁷G maps in human cell lines. Importantly, we show that certain internal m⁷G sites accumulate relatively high m⁷G fractions, suggesting functional roles. We identified METTL1 as one writer that can install a subset of internal m⁷G within mRNA and show that m⁷G might facilitate translation promotion of METTL1-affected transcripts. We propose that m⁷G could exhibit diverse functions on target mRNA, perhaps through modulating structure and protein-RNA interactions. This work provides critical resources for future investigations of m⁷G that may impact diverse biological processes and human diseases.

CONTACT FOR REAGENT AND RESOURCE SHARING

Requests should be directed to and will be fulfilled by the Lead Contact, Chuan He (chuanhe@uchicago.edu), after execution of a Materials Transfer Agreement.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

METHOD DETAILS

QUANTIFICATION AND STATISTICAL ANALYSIS

DATA AND SOFTWARE AVAILABILITY

RNA sequencing data listed in Table S1 are available in the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo) under accession number GSE112276.