Transcriptome-wide mapping reveals widespread dynamic-regulated pseudouridylation of ncRNA and mRNA

Abstract

Pseudouridine is the most abundant RNA modification, yet except for a few well-studied cases, little is known about the modified positions and their function(s). Here, we develop Ψ-seq for transcriptome-wide quantitative mapping of pseudouridine. We validate Ψ-seq with spike-ins and de novo identification of previously reported positions and discover hundreds of unique sites in human and yeast mRNAs and snoRNAs. Perturbing pseudouridine synthases (PUS) uncovers which pseudouridine synthase modifies each site and their target sequence features. mRNA pseudouridinylation depends on both site-specific and snoRNA-guided pseudouridine synthases. Upon heat shock in yeast, Pus7p-mediated pseudouridylation is induced at >200 sites, and PUS7 deletion decreases the levels of otherwise pseudouridylated mRNA, suggesting a role in enhancing transcript stability. rRNA pseudouridine stoichiometries are conserved but reduced in cells from dyskeratosis congenita patients, where the PUS DKC1 is mutated. Our work identifies an enhanced, transcriptome-wide scope for pseudouridine and methods to dissect its underlying mechanisms and function.

Figure 1. Ψ-seq quantitatively measured transcriptome-wide pseuodouridylation profiles

(A) Ψ-seq procedure. Poly(A) selected RNA is treated with CMC, which covalently binds to U and residues, and, to limited extent, G (Ho and Gilham, 1971; Metz and Brown, 1969a, b), followed by incubation at alkaline pH, leading to hydrolysis of U-CMC adducts, which are less stable than Ψ-CMC counterparts. RNA is next fragmented to a size range of 80–150 nt, followed by adapter ligation to the 3′ end. Reverse transcription is primed off of the adapter, expected to lead to premature termination one base pair immediately downstream of pseudouridylated sites. A second adapter is ligated to the 3′ end of the cDNA and libraries are amplified and sequenced in paired end mode. (B) Scoring modified sites. Ψ-ratio (red, top) between reads terminating at a site and reads overlapping it, and Ψ-fc (blue, bottom), the fold change in Ψ-ratio in the CMC-treated sample over non-treated control. (C) Validation with a synthetic spike in. Ψ-fc (Y axis) at each position (X axis) across a synthetic spike-in harboring one pseudouridylated site at varying stoichiometries (color legend) at the indicated position (arrow). (D) Ψ-seq quantifies relative stoichiometry. Scatter plot comparing Ψ-ratios (Y axis) at the pseudouridylated site in the spike-in with pseudouridine stoichiometries (X axis). The Pearson correlation coefficient (r) and associated P value are denoted. (E) Ψ-seq detects known Ψ sites in rRNA. Ψ-fc values (Y axis, dark blue) at each position (X axis) in 18S rRNA, overlaid with all known Ψ sites in this subunit (grey vertical lines). (F) ROC curves for different metrics (color legend) for calling putative Ψ sites. Each classifier was trained based on Ψ-seq data from human rRNA, and tested on their performance in yeast rRNA. (G) Nucleotide distribution across all sites passing our detection thresholds in a sample in mid-log yeast. See also Figure S1.

Figure 2. Identification of PUS-dependent Ψ sites in yeast

(A) Ψ-seq correctly distinguishes targets of snoRNAs. A comparison of Ψ-fc scores across all known Ψ sites in three rRNAs (18S, 28S, 5.8S) measured in yeast strains in which either H/ACA box snR34 (X axis) or snR189 (Y axis) were deleted. snR34 is known to target positions 2880 and 2886 in 25S rRNA (red), snR189 targets positions 2735 in 25S and 466 in 18S rRNA (purple), all detected as outliers in these experiments. (B) Ψ-seq associated modified bases with their cognate PUS and snoRNAs. Heatmap depicting Ψ-ratios (Z scores per row, color bar) across all sites (rows, gene and position labeled on left) in the PUS-dependent collection that are dependent on a snoRNA or a non-essential site-specific PUS (columns). For tRNAs, the labeled position is with respect to a multiple alignment of tRNAs with manual correction at two sites, to allow comparison between sites. A check mark indicates that a site was both previously known to be pseudouridylated and that its catalysis was known to be mediated through the enzyme or snoRNA with the minimal Z score in the figure. (C) Ψ-seq identifies Cbf5-dependent Ψ sites. As in (B) for all sites (rows) dependent on Cbf5. (D) Distribution of Ψ sites in different RNA classes, in the Pan-Ψ collection (left) and PUS-dependent collection (right). CDS: mRNA coding sequence. (E) Ψ-ratios (Y axis) across the snoRNA snR47 (left), the mRNA GLK1 (middle), and the mRNA SPI1 (right), which depend on Cbf5, Pus1 and Pus7, respectively. Ψ-ratios are presented for CMC-treated (black) and non-treated (red) samples. Called Ψ sites are marked in vertical thick grey lines, and their ‘absence’ position (in the relevant deletion strain) is marked by a dashed grey box. See also Figure S2 and Table S1.

Figure 3. Pseudouridylation in yeast snoRNAs

(A) Complementarity regions between three C/D box snoRNAs in yeast (black) and target sites in rRNA (red). Red arrows: detected Ψ sites. (B) As in A, but in human. (C) As in A, but for H/ACA box snoRNAs. (D) SnoRNAs expression in WT strains (X axis) compared to Cbf5 knockdown strains (Y axis), showing decreased expression of H/ACA box snoRNA (red), but not C/D box snoRNAs (orange) or snRNAs (dark red), upon knockdown. (E) Cumulative distribution frequencies of minimal Z-score transformed free energies, predicted by cofolding each putative CBF5-dependent (or shuffled) site against all known H/ACA box targeting arms. Distributions are plotted separately for sites in the rRNA (green), all remaining sites in mRNA and snRNAs (blue), and their corresponding shuffled versions (yellow and red, respectively). (F) Predicted associations between three CBF5-dependent sites (red) and top-scoring H/ACA box snoRNAs (black).

Figure 4. Induction of Pus7-dependent sites in heat shock

(A) Bar plots of the number of differentially modified sites (Y axis), in each of three conditions (X axis) relative to mid-log yeast. Sites are color-coded based on RNA class (color legend). (B) Differential modification in heat shock is not merely due to changes in expression levels. Expression levels in heat shock (X axis) or non-heat shock (Y axis) conditions for all genes harboring heat-shock induced Ψ sites. (C) Sequence motif learned across all Ψ sites induced in heat shock, is strongly enriched for the TGTAR Pus7 motif. (D) Ψ-ratios (Y axis) in CMC-treated (black) and non-treated (red) samples in five genes (columns) acquiring Ψ sites in heat shock. Ψ-ratios are shown in WT in 30° (top) and heat-shock (45°, middle), and in heat-shock in Δpus7 (bottom). Called sites are depicted in vertical grey bars in the samples in which they are called, and dashed grey boxes in the remaining samples. (E) Box plots of Ψ-ratios in WT or Δpus7 strains, grown in heat shock or 30°, for 159 heat-shock induced sites harboring a Pus7 consensus (left) and 106 sites without the consensus (right). (F) Cumulative distribution plots of fold changes of expression levels between WT and Δpus7 strains (X axis), comparing genes harboring a Pus7-induced Ψ sites (blue) to genes lacking it (red). At 45° (left), genes harboring Ψ sites are expressed at higher levels in the WT strain than in Δpus7, whereas non-pseudouridylated genes are present at similar levels. At 30° (right), these differences are marginal. P-value of the Kolmogorov-Smirnov test is depicted. (G) Pus7 transcript (left; mean TMM-normalized FPKM values from three replicates) and protein (right, data in four replicates from (Nagaraj et al., 2012)) levels (Y axis) in heat shock and non-heat shock conditions. Error bars are standard error of the mean.(H, I) Pus7 localization. Fluorescent microscopy images (H) and quantification (I) during a heat shock time course show mostly nuclear localization at 30°C (time 0 in H, and red in I), which is significantly reduced after heat shock (time 60 in H, blue in I). P-value of a Mann-Whitney test is noted. See also Figure S3 and Table S2.

Figure 5. Ψ-seq of human RNA highlights conserved features and disease-relevance

(A–B) Sequence motifs from 70 sites harboring a Pus4 associated GUUC core (A) and 13 sites harboring a Pus7 associated UGUA core (B), of 396 modified sites in HEK293 RNA. (C) Difference in Ψ-ratios in DKC1 knockdown relative to control (Y axis) in HEK293 cells. Box plots from left to right: DKC1-dependent sites, sites harboring Pus4 (GUUC) or Pus7 (UGUAR) associated motifs, and all remaining sites. (D–F) Correlation of pseduoridylation stoichiometries in rRNA between homologous positions in S. cerevisiae and C. albicans rRNA (D), human and mouse (E), and yeast and human (F). Positions are color-coded based on rRNA subunit. (G) Distributions of Ψ-ratios across rRNA positions in fibroblasts from dsykeratosis patients (grey) and age-matched controls (pink). Left panel: 7 year old patient; Right: 11 year old patient. (H) Ψ-ratios in rRNA positions in DKC1 knockdown (grey) or control (pink) in HEK293 cells. (I) TERC modification is affected by DKC1 knockdown. Number of reads (top; indicative of RT termination) and Ψ-ratios (bottom) at each site in the TERC transcript (X axis). Reads are shown in both CMC-treated (black) and non-treated (red) samples; Ψ-ratios for the treated samples. The putative Ψ position is highlighted in a grey vertical bar. (J) Ψ-ratios for position 307 in 7-yr old patient and age matched control; Error bars: standard error. (K) RNA secondary structure of the CR4/CR5 region in TERC redrawn based on (Zhang et al., 2011). Red arrow: putative Ψ site. See also Figure S4 and Table S3.

Figure 6

(A) mRNA pseudouridylation in eukaryotes is mediated by at least four conserved PUSs, some site specific (Pus1, Pus4, Pus7) and others snoRNA mediated (yeast Cbf5/human DKC1). (B) Pus7 may orchestrate the yeast mRNA pseudouridylation program in heat shock. At 30°, Pus7 is primarily nuclear, and its localization into the cytoplasm upon heat shock may induce mRNAs pseudouridylation.