Transcriptome-wide profiling of multiple RNA modifications simultaneously at single-base resolution

Abstract

The breadth and importance of RNA modifications are growing rapidly as modified ribonucleotides can impact the sequence, structure, function, stability, and fate of RNAs and their interactions with other molecules. Therefore, knowing cellular RNA modifications at single-base resolution could provide important information regarding cell status and fate. A current major limitation is the lack of methods that allow the reproducible profiling of multiple modifications simultaneously, transcriptome-wide and at single-base resolution. Here we developed RBS-Seq, a modification of RNA bisulfite sequencing that enables the sensitive and simultaneous detection of m⁵C, Ψ, and m¹A at single-base resolution transcriptome-wide. With RBS-Seq, m⁵C and m¹A are accurately detected based on known signature base mismatches and are detected here simultaneously along with Ψ sites that show a 1-2 base deletion. Structural analyses revealed the mechanism underlying the deletion signature, which involves Ψ-monobisulfite adduction, heat-induced ribose ring opening, and Mg²⁺-assisted reorientation, causing base-skipping during cDNA synthesis. Detection of each of these modifications through a unique chemistry allows high-precision mapping of all three modifications within the same RNA molecule, enabling covariation studies. Application of RBS-Seq on HeLa RNA revealed almost all known m⁵C, m¹A, and ψ sites in tRNAs and rRNAs and provided hundreds of new m⁵C and Ψ sites in noncoding RNAs and mRNAs. However, our results diverge greatly from earlier work, suggesting ∼10-fold fewer m⁵C sites in noncoding and coding RNAs and the absence of substantial m¹A in mRNAs. Taken together, the approaches and refined datasets in this work will greatly enable future epitranscriptome studies.

Keywords: RNA methylation; RNA modification; m1A; methyl adenosine; pseudouridine.

Fig. 1.

RBS-Seq enables simultaneous base-pair resolution transcriptome-wide mapping of m⁵C, m¹A, and Ψ. (A) Schematic of reactivity of modified nucleotides to bisulfite (Top) and the principle of simultaneous mapping of modified nucleotides (Bottom). For m⁵C, bisulfite treatment deaminates Cs (converting to Us; Ts upon cDNA sequencing), whereas m⁵Cs resist bisulfite (remain Cs upon cDNA sequencing) establishing sites of cytosine methylation. For m¹A, cDNA synthesis at sites with m¹A often confers misincorporation/mismatch in the NBS sample. In contrast, during treatment, m¹A becomes converted to m⁶A via Dimroth rearrangement (methyl passes from N¹ to N⁶), which faithfully templates thymine (cDNA remains adenosine) in the BS sample. Thus, comparison of NBS and BS samples identifies m¹A sites. For Ψ, Ψ nucleotides upon bisulfite treatment form a stable monobisulfite adduct (Fig. 4) causing frequent bypass with reverse transcriptase, leaving a deletion signature at the exact modified sites, evident exclusively in the BS samples. (B) Schematic representation of tRNA^Gly, indicating its well-known multiple m⁵C and single m¹A and Ψ modified sites. (C) Bar graph summarizing the actual RBS-Seq results from HeLa cell line for a tRNA^Gly locus indicating the exact locations of the modified nucleotides and their levels. The low levels of m⁵C shown at positions 40 and 60–66 have been shown previously for a subset of tRNA types (31, 32).

Fig. 2.

Transcriptome-wide RBS-Seq analysis in HeLa cells shows widespread m⁵C and Ψ in mRNAs and ncRNAs and m¹A almost exclusively in rRNAs and tRNAs. (A) HDGF mRNA bears a single m⁵C base, within its 3′ UTR. (B) Distribution and dynamic range of m⁵C sites in coding and noncoding RNAs, stratified by the nonconversion rate (percent of reads that contain a single cytosine that remains cytosine after bisulfite treatment). (C) Mismatch rate (percent of reads with a base mismatch) at human 28S rRNA transcript reveals a single well-known m¹A nucleotide at position 1322. (D) Mismatch rate at the conserved well-known A58 sites of human tRNAs (from HeLa cells) in NBS and BS samples. Representatives of each tRNA type with the highest mismatch rate are shown. (E) Deletion rate (percent of reads bearing a 1–2 base deletion) along CDC6 mRNA reveals a single Ψ nucleotide. (F) Distribution and dynamic range for Ψ candidate sites in coding and noncoding RNAs.

Fig. 3.

(A) Distribution of candidate m⁵C sites in RNA species, with expanded mRNA annotations. (B) Distribution of candidate Ψ sites in different RNA species, with expanded mRNA annotations. See Datasets S1 and S5 for all sites and their corresponding annotations.

Fig. 4.

Characterization of a pseudouridine-bisulfite adduct and heat/Mg²⁺-induced rearrangement to elicit reverse transcriptase bypass. (A) Sequence and intramolecular folding of pseudouridylated 70-mer RNA oligonucleotide used in the downstream experiments. The two Ψ sites (in red) are indicated by arrowheads. (B) Flowchart of oligo treatments, RT-PCR, TA cloning, and Sanger sequencing of individual colonies. (C) Summary of the deletion signatures obtained from oligonucleotide experiments with the reference sequence and the two Ψ sites at the bottom, showing the insufficiency of the bisulfite step, and requirement for the subsequent heat + MgCl₂ step to generate the deletion signatures. (D) Sequence and calculated mass for 12-mer control (12-U) and pseudouridylated (12-Ψ) oligomers used in the downstream experiments. (E) Reaction sequence and methods used for Ψ reactivity studies with 12-mers. (F) Mass spectrum for 12-Ψ after bisulfite and subsequent heat + MgCl₂ treatments shows formation of a stable monobisulfite adduct. Mass spectrum for the 12-U and 12-Ψ with only bisulfite treatment is provided in SI Appendix, Fig. S24. (G) A proposed model showing that during cDNA synthesis, ribose ring-opened Ψ-monobisulfite is oriented away from the polymerization site, reinforced by Mg²⁺, explaining base skipping.