Enzyme-mediated alkynylation enables transcriptome-wide identification of pseudouridine modifications

Domain: biology
Published: 2026-03-23
Editorial outcome: major_revision

Abstract
Pseudouridine (Ψ) is one of the most abundant chemical modifications and plays important roles in RNA function. Advances in our understanding of Ψ have been hindered by a limit of robust methods to precisely and sensitively map their distributions in cellular RNAs. Here, we present ELAP-seq (Enzymatic Labeling and Pull-down for Sequencing) for Ψ detection, which leverages a naturally occurring N1-methyl pseudouridine methyltransferase from Methanocaldococcus jannaschii (Mj1640). This enzyme promiscuously converts Ψ to N1-methyl-Ψ (m1Ψ) or installs a propargyl group at the same location in vitro under a mild condition, exhibiting high sensitivity and specificity, and is also functional inside cells. ELAP-seq enriches Ψ-containing RNA fragments and enables single-nucleotide-resolution Ψ detection with markedly enhanced signal-to-noise ratio and reduced sequencing and computational demands. Using ELAP-seq, we identify over 5,000 candidate Ψ sites in human HeLa and HEK 293T transcriptomes, validating many previously identified sites as well as reporting additional ones. This versatile enzymatic platform expands the toolkit for sensitive labeling and detection of Ψ, advancing the study of RNA modification biology.

Introduction
Chemical modifications of RNA are critical to most RNA functions1. Among the over 150 types of RNA chemical modifications, pseuoduridine (Ψ) is one of the most abundant, and exists in multiple types of RNA species2,3. Ψ is formed via the isomerization of uridine, a reaction catalyzed either by various stand-alone Ψ synthases (PUS) or by the DKC1 complex guided by H/ACA box small nucleolar RNAs (snoRNAs)3,12. The enzymatic cores of the Ψ synthases are conserved, and their mutations have been associated with human diseases, including autoimmune digestive disorders and intellectual disability 13-15. The presence of Ψ in RNA can alter RNA properties including structure, stability, alternative splicing and translation4-7, and Ψ modifications can have a profound effect on cellular activities such as cell proliferation and stem cell differentiation 8,9. Ψ has also been found to exist in viral RNAs and may play roles in viral infection and antiviral defense pathways10. Furthermore, Ψ and its derivative N1-methylpseudouridine (m1Ψ) are among the most critical nucleotide substitutions that notably promote efficacies of mRNA vaccines11. The correlation of Ψ with human diseases, along with its therapeutic potentials in developing RNA-based drugs, have inspired enormous interest to investigate their natural presence and biological functions in cells.

While Ψ has been extensively studied in the context of abundant rRNA, tRNAs and ncRNAs for decades, its significance in less abundant mRNAs remains poorly understood, largely due to a lack of robust methods for sensitive and precise detection of Ψ modifications in mRNA species. Ψ shares the same base pairing property as uridine, making it indistinguishable in Sanger and next-generation sequencing. Existing detection strategies largely rely on chemical treatments that convert Ψ into forms generating stop, deletion or mutation signatures at or near the modification sites during reverse transcription 16-20. Notable recent examples include RBS-seq, BID-seq and PRAISE, which utilize bisulfite treatment, and BACS, which employs N1-specific Michael addition of acrylonitrile, for qualitative or quantitative transcriptome-wide Ψ mapping 16,20-22. While accurate, these methods require deep sequencing due to a lack of enrichment, which leads to high sequencing cost and demanding computing capacity and may also result in missing sites in lowly expressed genes and/or with lower modification fractions. In addition, methods based on bisulfite treatment rely on nucleotide deletion signatures, making it difficult to determine the exact location of the modification among consecutive uridines. An enrichment-based method using CMC addition has been developed, but prolonged chemical treatment led to severe RNA degradation therefore requiring a high level of RNA input, and incomplete CMC release from non-specific addition to unmodified nucleobases could lead to a non-negligible rate of noise signals especially under the enrichment conditions19. Meanwhile, direct Ψ detection via nanopore sequencing has emerged, yet it remains limited by relatively low resolution, high RNA input requirement, and complex data analysis pipelines 23,24. Current transcriptome-wide studies have reported hundreds to a few thousands of modification sites. However, LC-MS/MS experiments suggest that Ψ abundance in mRNA is comparable to N6-methyladenosine (m6A, 0.2-0.4% in HeLa and HEK 293T cell lines)19, implying that potentially more Ψ sites are to be discovered. Thus, new methods using independent strategies are highly desired to cross-validate previously reported modification sites as well as identify additional modification sites. Furthermore, existing labeling strategies are restricted to isolated RNAs, meaning that there is no feasible way to label Ψ inside cells so far.

To provide additional approaches for both in vitro and potential cell-based manipulation of Ψ, in this study, we examine several natural Ψ-specific methyltransferases for in vitro reactivity and focus on a N1-methyltransferase from Methanocaldococcus jannaschii that naturally targets Ψ in tRNA26. We find that this enzyme can install a methyl or propargyl group onto the N1 position of Ψ in various sequence contexts independent of RNA structures. The labeling strategy also exhibits higher specificity and sensitivity compared to currently available antibodies. This discovery provides a simple and mild one-step reaction to specifically yet broadly label Ψ under physiological conditions, and can potentially be used to modulate Ψ inside human cells, making it a highly promising and versatile strategy to detect and modulate Ψ. To showcase the utility of this labeling strategy, we leverage it to develop a method that enriches Ψ-containing transcripts for transcriptome-wide Ψ profiling. We validated the method by applying it to detect known Ψ modifications in rRNA and then used it to identify candidate Ψ modification sites in the transcriptomes of HeLa and HEK 293T cell lines. Thousands of candidate modification sites are identified for each cell line, which validate many previously reported sites while also revealing additional ones, providing insights into features of Ψ distribution and function in human cells. Altogether, our study provides a tool that allows studies of Ψ modification both in vitro and potentially in cells.

Results
A Ψ-specific methyltransferase that labels Ψ modifications
Currently, labeling tools to investigate RNA Ψ modification are limited. To address this gap, we explored natural enzymes that specifically act on Ψ. Natural enzymes often serve as valuable starting points to develop biotechnologies thanks to their exceptional substrate specificity, catalytic efficiency and ability to function inside cells 25. A plethora of proteins and enzymes exist that selectively recognize and act on RNA modifications; exploiting these biomolecules offers a promising direction for the modulation of RNA modifications 25.

One major physicochemical difference between uridine and Ψ is an extra hydrogen donor at the N1 position of the Ψ nucleobase (Figure 1a). We reason natural enzymes specifically reacting at this position would facilitate selective labeling of Ψ over uridine. N1-methylpseudouridine (m1Ψ) is a modification naturally found in tRNA of most archaea and is formed via methylation at the N1 position of Ψ, a reaction catalyzed by SPOUT-class methyltransferases27. We screened a few SPOUT-class methyltransferases that catalyze this conversion and focused our attention on an enzyme named Mj1640 from Methanocaldococcus jannaschii, which is responsible for methylation of Ψ at position 54 in the T-arm of tRNAs (Supplementary Figure 1a)26. Although the methylation was reported to occur naturally only on Ψ54 in tRNA, Wurm et al. reported that the enzymatic catalysis is highly efficient in vitro on a 17-mer RNA probe resembling the T-arm of Methanocaldococcus jannaschii tRNATrp at 65 °C26. We reasoned that the relatively high reaction temperature may have disrupted the secondary structure of the 17-mer RNA probe, thus the enzyme might be more promiscuous for labeling Ψ independent of local RNA structure.

We proceeded and tested activity of recombinant Mj1640 protein (Supplementary Figure 1b-e) in vitro on the same 17-mer synthetic RNA oligonucleotide as by Wurm et al.26. Encouragingly, the Ψ-RNA was methylated to > 90% in five minutes as evaluated by MALDI-TOF MS, whereas no reaction was observed for a control RNA with a U replacing the Ψ after over a Ψ after over a one-hour incubation (Figure 1b), validating the excellent specificity of the reaction towards Ψ. Importantly, a synthetic 8-mer oligo containing a Ψ site was also efficiently methylated (Supplementary Figure 1d) whereas the control oligo with U was not, demonstrating that the enzymatic catalysis does not rely on RNA secondary structure. Sequence preference of this enzyme was also studied and will be shown later.

In addition, we found that when expressed in human HEK 293T cells, the Mj1640 enzyme is functional to convert Ψ in endogenous transcripts to m1Ψ (Supplementary Figure 3b-d). Altogether, our study presents an enzymatic strategy for labeling Ψ both in isolated RNAs and inside cells.

Mj1640 can transfer a propargyl group to Ψ in RNA
Next, we tested whether the Mj1640 enzyme could use the SAM derivatives allylic-SAM or propargylic Se-adenosyl-L-selenomethionine (propargyl-SeAm) to label Ψ sites with chemical groups suitable for further manipulation. To our delight, Mj1640 efficiently utilized both as cofactors to label the 17-mer Ψ-RNA (> 80% in 1 hr) and the 8-mer Ψ-RNA, but not the control U-RNAs (Figure 1c and Supplementary Figure 1d-Ae).

To evaluate the ability of the enzymatic reaction to label Ψ modifications in biological RNA, we assessed the labeling efficacy on RNA extracted from HeLa and HEK293T cell lines using quantitative LC-MS/MS. We detected a dramatic decrease of Ψ signal and the emergence of N1-methylpseudourisine (m1Ψ) or N1-porpargyl-pseudouridine (p1Ψ) signal after the enzymatic treatment with either with SAM or propargyl-SeAm as the co-factor (Figure 1d and Supplementary Figure 2a). With optimized protein purification and enzymatic reaction conditions (see the Methods section), the Ψ in vitro methylation efficiency was around 50% for total RNA, and 33-48% and 50-60% for poly-A RNA of HeLa and HEK293T cells, respectively (Figure 1e, Supplementary Figure 2a-c). The overall 33-60% m1Ψ/Ψ ratio indicates that a considerable proportion of Ψ modifications were successfully labeled.

Additionally, total RNA and poly-A RNA were labeled with the propargyl group at ~25-30% of all Ψ (Figure 1e and Supplementary Figure 2d-g). This slightly lower labeling percentage compared to methylation is likely due to a weaker binding affinity of the enzyme to propargyl-SeAm than to SAM rather than a biased labeling. These results demonstrated the feasibility of using Mj1640 to label Ψ in biological RNA. It is important to note that the enzymatic reaction condition is mild, resulting in largely intact mRNAs with ~80% recovery rate (Supplementary Figure 3a). This observation highlights one advantage of this enzymatic labeling strategy over chemical labeling approaches that typically cause RNA degradation.

Mj1640 labeling assay exhibits superior sensitivity and specificity compared to antibodies
Pseudouridine has clear relationship with several diseases; however, the potential of both PUS activity and global modification signatures as biomarkers for disease diagnosis has yet to be effectively evaluated. The click chemistry enabled by our enzymatic assay allows flexible attachment of different functional groups for signal readout, suggesting its potential suitability for diagnostic or detection applications.

To evaluate this potential, we characterized the specificity, sensitivity, and quantitative performance of the enzymatic labeling strategy by clicking either a fluorophore or biotin onto the modification for fluorescence- or chemiluminescence-based dot blot analyses. We observed a strong linear correlation between the amount of modified RNA and the detected signal (Supplementary Figure 4a-c), and high specificity was achieved in both synthetic and total RNA samples (Supplementary Figure 4a-c). Importantly, the chemiluminescence-based dot blot assay can detect as little as 0.05 ng of synthetic RNA oligos (one modification in 33-mer RNA oligonucleotides) (Supplementary Figure 4c), demonstrating substantially higher specificity and sensitivity than currently available antibodies (Supplementary Figure 4d-e). Together, these preliminary results indicate that this enzymatic labeling strategy has strong potential to serve as the basis for diagnostic and detection tool development.

Development of the Enzymatic Labeling and Pull-down for sequencing (ELAP-seq) method for detecting Ψ in the transcriptome
To apply the enzymatic labeling strategy to study Ψ modification in the transcriptome, we developed a method coupling Enzymatic Labeling and Pull-down for Sequencing (ELAP-seq) (Figure 2a). Compared to methods using direct sequencing, ELAP-seq enriches the modified RNA for sequencing, which provides the advantages of effectively increasing signal-to-noise ratios, reducing sequencing and computational cost, and allowing for more comprehensive mapping of Ψ even in lowly expressed mRNA species or those with low stoichiometry. The enzyme’s exclusive specificity towards Ψ is critical to ensure the specific enrichment of Ψ-modified RNA, thereby reliable identification of the modification.

In this protocol, biological RNA such as total RNA or mRNA is first fragmented to 60-200 nucleotides, followed by enzymatic treatment to label the Ψ with a propargyl group (Figure 2a). An aliquot of fragmented RNA set aside prior to enzymatic treatment is used to directly prepare a sequencing library to control for signals from factors other than Ψ existing in the starting sample (input library) (Figure 2a). The enzymatically labeled RNA is then split into two parts, with one part undergoing a copper-catalyzed azide–alkyne cycloaddition (CuAAC or click) reaction to add on a biotin molecule while the other part undergoing a mock click reaction (control library). Pull-down is then performed with streptavidin-coated magnetic beads to enrich RNA fragments containing Ψ modifications, followed by library construction and Illumina sequencing (pull-down library) (Figure 2a). The control library is expected to yield no detectable DNA product under the same conditions used for pull-down library construction and is not subjected to sequencing.

To evaluate the signal-to-noise ratio of the ELAP strategy, we treated synthetic RNA oligonucleotides bearing either a Ψ or U in the center with enzymatic labeling, biotin addition, and pull-down sequentially, and then quantified the amount of the resulting RNA material via RT-qPCR. With an optimized pull-down buffer condition, the Ψ-RNA was enriched 4~450-fold relative to the U-RNA with input ranging from 0.01 fmol to 10 pmol (Figure 2b), demonstrating an excellent signal-to-noise ratio.

ELAP-seq reveals sequence preference of Mj1640 and detects Ψ at single base resolution
Before applying ELAP-seq to biological samples, we first used it to comprehensively evaluate the sequence context preference of the Mj1640 enzyme by building sequencing libraries from a mixture of 33-mer synthetic standard RNA oligonucleotides containing NNΨNN sequences (N signifies any nucleoside) which cover 256 sequence contexts surrounding Ψ. As expected, the control library yielded no sufficient material to produce a sequencing library in the end unless being over-amplified, highlighting the low background. Interestingly, sequencing data from the pull-down libraries showed that RNAs with 5’-ΨU-3’ context were preferentially enriched while marginal sequence preferences were observed at other nearby sites (Figure 2c). In contrast, control RNA oligonucleotides containing NNUNN sequences, whose sequencing libraries were built with overamplification, showed no such specific enrichment of any RNA (Supplementary Figure 5), further demonstrating the high specificity of this method.

To further validate the observed sequence preference, we synthesized 10-mer RNA oligonucleotides containing degenerate A, U, C and G (denoted as N) at one of the four positions within two nucleotides upstream or downstream of the target Ψ site, and evaluated the enzymatic reaction efficiency on these oligonucleotides using MALDI-TOF MS. In agreement with the ELAP-seq results, Mj1640 does not show obvious selectivity for nucleobases at the -2, -1 and +2 positions but exhibits preference for uridine over other nucleobases at the +1 position (Supplementary Figure 6). Importantly, the strong concordance between the ELAP-seq result and the MALDI data underscores the reliability of the ELAP-seq approach.

These results are consistent with a previous report that Mj1640 requires 5’-Ψ (U/ Ψ)-3’context for optimal catalysis 27. Earlier research indicates that ~30% of Ψ sites in mRNA are found in the 5’-ΨU-3’ context 18,19, the preferred sequence context of Mj1640, as supported by our observation that the Mj1640 converted Ψ in cellular mRNA to m1Ψ at a level of 30-60% (Figure 1e and Supplementary Figure 2). Despite this sequence preference, our enrichment protocol allows pulling down oligonucleotides of almost all sequence contexts to various extents (Figure 2c). Meanwhile, two independent pull-down experiments of the standard oligos resulted in highly consistent results (Linear correlation with R2=0.98) (Figure 2d), demonstrating excellent reproducibility of the pull-down protocol.

Pull-down of fragmented RNA for sequencing typically yields Ψ signals at a resolution of 60-200 nt. To achieve detection of Ψ at single base resolution, we found that a 500x lower dNTP concentration than recommended in the commercial manual led to significant arrest of superscript III or M-MLV at the alkynylated Ψ site with or without the addition of a biotin molecule (Figure 2e). Collectively, the enrichment of RNA fragments containing Ψ in various sequence contexts, coupled with the reverse transcription (RT) stop signature at the labeled sites, would enable ELAP-seq to specifically identify Ψ modification sites at single base resolution.

Validating ELAP-seq by detecting known Ψ modifications in rRNA
Next, we tested the efficacy of ELAP-seq in the context of biological samples. We first evaluated its ability to detect known Ψ modifications in ribosomal RNA 32. We used superscript III with the low dNTP concentration that facilitates the generation of stop signatures as shown in Figure 1e. As expected, reverse transcription mediated by superscript III was substantially arrested after putting an adenosine opposite the Ψ, generating a drastic stop signature right at the site of the modification (Figure 3a). Importantly, Ψ modifications were specifically identified based on the stop signatures uniquely generated after enzymatic treatment (Figure 3a and Supplementary Figure 7b), whereas Nm sites in rRNA were effectively filtered away based on the same stop signature in the input library (Figure 3a and Supplementary Figure 7b-c).

Our method successfully identified 72 out of 104 known Ψ sites in rRNA (69.2% sensitivity) at single base resolution, covering 29/29 sites in the 5’ΨU3’ context and 41/75 sites in other contexts (Figure 3b-c). Our method also led to stop signatures at U1615 of 18S, whose originalities are unclear, as well as U31 and U630 of 18S and U1603 of 28S, which likely originated from off-target mapping of short intronic reads to rRNA genes (95% accuracy) (Supplementary Data 5). The known Ψ sites undetected by our method were primarily due to the clustering of these sites in rRNA (Supplementary Figure 7d), leading to low stop signature below the detection limit of our method.

ELAP-seq identifies thousands of candidate Ψ sites from the human transcriptome
We next applied ELAP-seq to poly-A RNAs purified from HeLa and HEK 293T cell lines (Supplementary Figure 9a). Supporting the specific enrichment of modified RNAs, we found selective enrichment of peaks in the pull-down samples, as exemplified by the maps on the mitochondrial chromosome (Figure 4a). In addition, candidate Ψ modification sites are clearly indicated by prominent stop signatures uniquely existing in the pull-down samples and absent in the input samples (Figure 4b).

We applied the same stop ratio and read coverage parameters that were successfully used for detecting rRNA modification sites to identify candidate modification sites in mRNA (See Supplementary Figure 8 and the Methods section). With our pipeline, we identified 3,738 RT stop sites from HeLa poly-A RNA libraries constructed with superscript III, among which 3,600 sites were located at U (96.3%), 138 at non-U (3.6%) (Figure 4c). Similarly, in HEK 293T cells, 3,473 RT stop sites were identified, with 3,331 at U (95.9%), 142 at non-U (Supplementary Figure 9b). The U sites showing stop signatures were further defined as candidate modification sites. Here we evaluated stops across non-U nucleotides to provide a measure of specificity and potential false positives. This is based on the rationale that, if no specific labeling and enrichment of pseudouridine occurs, RT stop would occur randomly with similar frequencies at A, U, C and G. Thus, the small fraction of non-U stops provides an empirical estimate of background RT termination events and the predominant enrichment of the stop signatures specifically at U compared to any other nucleotide (64-fold), which was not seen in the control input samples (Figure 4e), strongly supports the specific detection of pseudouridine33.

While Superscript III with low dNTP produces prominent RT stop signature, we noticed that the coverage at some sites in mRNA is below the reads coverage cutoff although both replicates share strong stop pattern. To achieve higher read coverage, we also used superscript IV with the commercially recommended concentration of dNTP (0.5 mM final) to generate libraries in parallel from the same pull-down product. Interestingly, superscript IV also generated stop signatures at Ψ, although to a lesser extent (Supplementary Figure 7a), suggesting that the on-bead RT reaction used in our library construction protocol further facilitated the generation of the stop signature (Supplementary Figure 7a). To rescue additional sites by increasing the reads coverage, we combined data from superscript IV and III for further analysis. We first tested this analysis strategy on rRNA and verified that, with slightly more stringent cutoff in stop ratio (See Methods section), this strategy did not introduce additional false positives. Additionally, this analysis enabled the identification of the two modification sites in 5.8S rRNA, which were previously missed due to low coverage (Supplementary Figure 7e).

Applying the same analysis on mRNA identified an additional 1,628 and 2,410 candidate sites from HeLa and HEK293T poly(A) RNA, respectively, yielding a total of 5,366 candidate sites for HeLa and 5,741 for HEK (Figure 4d, Supplementary Figure 9b). Of these, 2, 531 candidate Ψ sites are shared between HEK293T cells and HeLa cells, accounting for ~46% of all candidate sites in each cell line (Figure 4f). We note that some sites passing the filter in one cell line but found in only one replicate in the other cell line were not included. Therefore, the actual number of shared sites between the two cell lines may be higher than reported here.

The Pearson correlations for normalized read coverage RPM of identified Ψ sites between two pull-down replicates is higher than 0.98 (Supplementary Figure 10a). We performed a third biological replicate, and the three replicates repeatedly identified 3,820 modification sites in HEK293T (and a similar number identified for HeLa cell line) (Supplementary Figure 10b-c). Notably, RPM values for 81% of the 3,820 sites shared by all three replicates exhibit excellent coefficient variation values of between 0 and 0.25 across the three pull-down experiments (Supplementary Figure 10e), suggesting high consistency in pull-down efficiency across biological replicates. The sites identified by only one replicate exhibited significantly lower numbers of stopped reads, with most of them falling below the detection threshold in the other replicates (Supplementary Figure 10d). This indicates that the discrepancies mainly arise from sites with near-threshold stopped reads counts. In this study, we define candidate modification sites as those detected in at least two independent replicates, which effectively enriches for sites with robust signals.

Validation of modification sites
To validate the authenticity of the candidate Ψ sites, we used a modification-free HEK mRNA transcriptome prepared via in vitro transcription (IVT) of cDNA of poly-A RNA as a negative control. We first confirmed that the IVT libraries covered gene bodies similar to the native biological poly-A RNA (Supplementary Figure 11a-b). With enzymatic labeling and pull-down, ELAP-seq could only identify 412 sites from the IVT RNA libraries, with 101 sites located at U (Supplementary Figure 11c). Of the 101 candidate Ψ sites detected from the IVT RNA libraries, 11 overlapped with the 5,741 candidate Ψ sites detected in native biological poly-A RNA samples (Figure 4g) and were excluded from the final list. This demonstrates that the signals we identified in native RNAs are not artifacts from the labeling or sequencing processes but true modification sites.

We also performed siRNA knock down for one of the pseudouridine synthases, DKC1, whose expression perturbation has been found to affect a large fraction of modification sites on mRNA20. We confirmed the successful RNA and protein knockdown by RNA-seq and western blot, respectively (Supplementary Figure 12a). Enrichment levels, which were used to compare differential modifications between conditions (see the ‘semi-quantification capacity of ELAP-seq’ section below), showed great correlations between the two biological replicates under both conditions (Supplementary Figure 12b-c). We identified 654 sites downregulated upon DKC1 depletion (p-value < 0.05 and Log (Fold Change) < -1) (Figure 4h and Supplementary Data 6), including 524 sites newly reported in this study. The changes in expression of genes harboring these downregulated sites are mild and unlikely to substantially affect the efficiency of pull-down of this modification or the reliability of the detected Ψ downregulation events (Supplementary Figure 12d-g). Notably, the DKC1-dependent sites determined by PRAISE are also predominantly downregulated in our dataset, supporting the reliability of the results (Supplementary Figure 12h). More studies are required to identify modification sites responsive to knockdown of other pseudouridine synthases.

Lastly, we set out to validate individual candidate Ψ sites using a low-throughput CMC-based RT-qPCR method (Supplementary Figure 13a). Among the 20 sites producing reliable qPCR signals, 17 were successfully validated, including 9 sites not previously reported (Supplementary Figure 13b). The ones not successfully validated likely have low levels of modification or labeling that are below the detection limit of the CMC RT-qPCR method without enrichment.

To better evaluate the reliability of individual candidate Ψ sites, we aimed to categorize them into different confidence levels. The non-U sites, which represent noises, exhibited markedly lower stop ratios compared to the U sites (Supplementary Figure 14a). To define higher confidence sites, we set a higher threshold for stop ratio (0.3, the median stop ratio of non-U sites) to (Supplementary Figure 14b). Additionally, sites detected in a greater number of replicates exhibited a reduced rate of false positives (Supplementary Figure 14c). We thereby classified candidate Ψ sites into the following categories: 1) highest-confidence (stop ratio > 0.3 in every replicate in which it is detected and detected in all three replicates, n = 3428 sites); 2) higher-confidence (stop ratio > 0.3 in every replicate it is detected, or detected in all three replicates, n = 1704 sites); and 3) lower-confidence (stop ratio < 0.3 in at least one replicated it is detected, and detected in only two replicates, n = 618 sites) (Supplementary Figure 14d and Supplementary Data 1 and 2).

Semi-quantification capacity of ELAP-seq
To evaluate the capacity of the method to probe differential modification at a site across different conditions, we calculated the enrichment level at each site by normalizing the local read coverage per million mapped reads (RPM) at the site in the pull-down sample to that in the input sample to account for gene expression difference. We found that these parameters exhibited excellent Pearson correlations greater than 0.95 between two biological replicates for both HeLa and HEK293T cell lines, demonstrating their reliability for comparing differential modification across conditions (Supplementary Figure 15a).

We further investigated intra-group modification level variations across different sites. We reasoned that the enrichment level is influenced by not only the modification level, but also the enzyme’s preference for the sequence context a site is in. To account for this, we adjusted the enrichment level by the sequence preference, which was determined by the corresponding synthetic oligonucleotide (Figure 2c). Due to the large discrepancy in enzymatic preference between the 5’ ΨU 3’ context and other sequence contexts, we categorized the sites into ΨU and non-ΨU groups for more accurate intra-group comparisons. Additionally, pull-down and qPCR on a synthetic oligonucleotide modified to various levels revealed an exponential correlation between the modification level and the enrichment level (Supplementary Figure 15b). We thus used the logarithmic form of the adjusted enrichment level as an estimate of the relative degree of modification across different sites.

Supporting the semi-quantification capability of the method, Ψ sites in rRNAs, which are known to be modified to nearly 100%32, are distributed at the higher end of the relative modification level spectrum in both groups (Figure 4i and Supplementary Figure 15c). Furthermore, we compared the relative modification levels estimated by our method to the absolute modification levels of 10 mRNA Ψ sites quantified by CLAP so far, either previously16, 34 or in the current study (Supplementary Figure 15d). Encouragingly, we found a correlation coefficient of 0.82 (Figure 4j). Reinforcing the semi-quantitative nature of the method, there exists a moderate-to-strong correlation between our quantification and that obtained via PRAISE (correlation coefficient of 0.45-0.79) (Supplementary Figure 15e). These results suggest that our method is reasonably accurate in distinguishing relative modification levels between sites.

Based on the correlation between our method and the quantification by CLAP (for Ψ in non-ΨU context, r = 0.82 (Figure 4j)) or PRAISE (for Ψ in ΨU context, r = 0.79 (Supplementary Figure 15e)), we further categorized the modification sites into highly modified, moderately modified and lowly modified groups (Supplementary Figure 15f and Supplementary Data 1 and 2). This shows that 72% of the modification sites are defined as lowly modified (estimated modification ratio < 30%), which is generally consistent with previous reports. Note that quantification may not be valid for clustered sites, as RT arrests could skew the actual enrichment level, leading to weaker correlation between enrichment and modification levels.

Comparisons to other methods
Our dataset overlaps with 1168 of 3818 sites (30%) collectively identified by PRAISE, BID-seq and CeU-seq (Supplementary Figure 16a) for HEK 293T cells. Specifically, our method successfully identified 46% (961) of the sites reported by PRAISE and 35% (192) of the sites reported by BID-seq (Supplementary Figure 16a), providing cross-validation for both strategies, while also determining the exact modification locations in multiple Us, a task that is challenging for PRAISE and BID-seq due to the RT deletion signature. Importantly, sites well overlapping with other methods were mainly found in the highest- and higher-confidence categories, supporting the validity of the confidence classification (Supplementary Figure 16b). We suggest caution when approaching sites in the lower confidence category (constituting 10% of all sites) and recommend additional validation or referring to Integrated Genome Viewer (IGV) to visually assess their reliability.

To address the discrepancy between our method and the other methods, we identified an additional 40.1% PRAISE-only sites (sites not overlapping with ELAP-seq) in our dataset that exhibit prominent RT stop signatures but either appear in only one replicate or have read coverage below our filter cutoff (Supplementary Figure 16c-d). Only 13.9% of PRAISE-only sites are entirely undetected by our method, and these sites exhibited lower deletion ratios and read coverage in PRAISE compared to the overlapping sites (Supplementary Figure 16c, e), suggesting that they are intrinsically more difficult to detect. In addition, our method detected 59% (191) of the TruB1-dependent modification sites reported by PRAISE, which are primarily located in the 5’ ΨC sequence context, demonstrating the capacity of ELAP-seq to detect modifications in diverse sequence contexts (Supplementary Figure 16f). The remaining TruB1-dependent sites undetected by our method exhibited lower deletion or stop ratios and read coverage in both ours and PRAISE dataset (Supplementary Figure 16g).

ELAP-seq-only sites (sites only detected by ELAP-seq) tend to have lower modification levels and read coverage than sites overlapping with other methods (Supplementary Figure 16h), suggesting that they could have been missed by other methods due to weak signals. We examined the raw data for 4060 ELAP-seq-only higher- or highest-confidence sites in the BID-seq dataset. Of these, 2567 sites (59%) exhibited deletion ratios above 0.05 in at least one replicate of BID-seq (Supplementary Figure 16i-j), demonstrating that they are modification signals in BID-seq as well, despite falling below the filtration cutoff. There are an additional 842 sites showing deletion ratios above 0 but below 0.05, with 89% of them having average deletion ratios of between 0.02 to 0.05 across three replicates (Supplementary Figure 16i-j). Meanwhile, 651 sites had a deletion ratio of 0 (Supplementary Figure 16i), but they exhibited significantly lower read coverage than the overlapping sites in the BID-seq dataset (Supplementary Figure 16k), suggesting that they were likely missed due to low gene expression.

Our method identified 17 candidate Ψ sites in HeLa mitochondria mRNA and 24 in HEK mitochondrial mRNA. This set of sites encompassed all five mitochondria mRNA sites for HEK 293T (chrM: 6293, 7329, 9905, 11154 and 14868) and three for HeLa (chrM: 6293, 7329, and 11154) collectively identified by BID-seq and PRAISE. Notably, our data shows that 15 mRNA modification sites were shared between the two cell lines, including sites 9905 and 14868, which were previously identified only for HEK 293T cell line (Supplementary Figure 17). It's worth noting that most of the mitochondria mRNA sites uniquely detected by ELAP-seq are characterized as lowly expressed or lowly modified (Supplementary Data 1- 4), suggesting potential reasons why they were missed by the other methods.

Altogether, our findings indicate that most discrepancies stem from sites that are intrinsically more difficult to detect due to their low modification or expression levels. This highlights the value of our method for sensitive detection and validation of weak signals overlooked by other methods.

Insights into Ψ biosynthesis and function
The 5,366 candidate Ψ sites identified from HeLa poly-A RNAs are distributed in 2,616 transcripts, of which 82% (4,403 sites) are located in messenger RNAs (Figure 5a), with enrichment in CDS and 3’ UTR and depletion in 5’ UTR, which agrees with previous studies (Figure 5b)16. Most modified transcripts (1,621) harbor only one Ψ, although some transcripts are modified at multiple locations, including FLNA (19 sites), A2M (16 sites), ATP1A1 (16 sites), GANAB (8 sites) and PRKDC (10 sites) (Figure 5c). Importantly, we found that Ψ exhibits higher trend to locate in mismatched and UG or GU wobble paired positions compared to unmodified U sites (Figure 5d and Supplementary Figure 18a), which is consistent with the previous report that Ψ installed by PUS7 tends to occur in less structured RNA35. This observation is unlikely due to a structural preference of the enzymatic labeling because: 1) we used a protocol that disrupts any remaining RNA secondary structures prior to enzymatic labeling; 2) sites in matched positions in rRNA were as effectively detected as sites in mismatched and wobble positions (Supplementary Figure 18b); and 3) the enzymatic labeling efficiency on isolated total RNA at 37 °C is the same as that at 65 °C (Supplementary Figure 18c), although the enzyme activity is higher at 65 °C.

We focused on moderately and highly modified mRNA transcripts for functional analysis. Gene ontology analysis of biological processes and molecular functions suggests a notable correlation of Ψ to cytoplasmic translation, mitochondrial electron transport and ATP synthesis (Figure 5e and Supplementary Figure 18d-f). Mutations in Ψ synthase PUS1 have been linked to mitochondrial myopathy. The identification of Ψ modifications on mRNAs encoding proteins critical for mitochondrial functions, such as UQCRC2, PMPCB, SLC25A5-6, and BNIP3, suggests a direct correlation between pseudouridylation and mitochondrial diseases. Additionally, the link of Ψ to double stranded DNA break repair is consistent with the previous observation that TruB1 depletion could retard cell growth and cause cell arrest in G1 phase5. More studies are required to thoroughly interrogate the determinants of the biosynthesis of Ψ and its biological function.

Methods
Protein purification
The Mj1640 plasmid (deposited in Addgene with deposit ID 87044) for overexpression in E. coli was a kind gift of the Jens Wöhnert laboratory at the Goethe-University Frankfurt. In the plasmid, the Mj1640 gene with an N-terminal His-tag optimized for overexpression in E. coli was cloned into the pET11a overexpression vector. The plasmid was expressed in T7 express competent E. coli cells (NEB) grown in LB medium and overexpression of the protein was induced by adding IPTG to 1 mM at OD6 around 1 and incubating for 3 h at 37 °C. Two liters of induced cell cultures were pelleted and resuspended in 60 mL lysis buffer composed of 50 mM sodium phosphate (pH 7.5), 500 mM NaCl, 10 mM imidazole, 5 mM beta-mercaptoethanol, and 1 tablet of complete protease inhibitor EDTA-free (Roche). Cells were lysed in French Press with 1-2 rounds of lysis under pressures between 15 000 psi and 20 000 psi. The purified lysate was run through a Ni-NTA column manually, washed with 300 mL wash buffer containing 2.5 mM sodium phosphate (pH 7.5), 12.5 mM NaCl and 30 mM imidazole, and eluted sequentially with three elution buffers containing 25 mM sodium phosphate (pH 7.5), 125 mM NaCl, and 100 mM, 200 mM and 300 mM imidazole, respectively. To thoroughly remove bound nucleic acids and contaminating proteins, the eluted protein was further purified sequentially with the anion exchange column twice and Heparin column once. Specifically, proteins were diluted with buffer A (20 mM Tris-HCl, pH 7.4) to lower NaCl concentration to below 75 mM and then loaded onto a FPLC- mono Q or -Heparin column. After washing with buffer A, the proteins were eluted with a gradient mixture of buffer A and buffer B (20 mM Tris-HCl pH 7.4 and 1M (for the Q15 column) or 1.5 M (for the heparin column) NaCl). The proteins eventually eluted from the Heparin column were buffer exchanged to the storage buffer (20 mM Na2HPO4 pH 7.5, 150 mM KCl), concentrated to 20 µL, flash frozen in liquid nitrogen and stored at -70 °C. The resulting protein was > 95% pure.

In vitro methylation reaction
The condition used for testing the reactivity on synthetic RNA oligos (ordered from Dharmacon) was 15 mM sodium phosphate (pH 7.5), 100 mM KCl, 500 µM SAM or pro-SeAm, 10 µM RNA substrate and 220 µM enzyme in a 50 µL scale reaction for various time periods at 65 °C. The SAM cofactor was purchased from New England Biolabs, and peo-SeAm was synthesized and HPLC purified by Dr. Mingkui Luo’s lab at Memorial Sloan Kettering Cancer Center, New York. After methylation, synthetic RNAs were purified via phenol-chloroform extraction and ethanol precipitation. The condition used for labeling isolated cellular RNA for sequencing library preparation was 15 mM sodium phosphate (pH 7.5), 100 mM KCl, 321 µM pro-SeAm (added from 1.5 mM stock of Pro-SeAm dissolved in water containing 0.5% TFA (trifluoric acid), with pH adjusted to 7 using Tris-HCl pH 8.3), 12.5% PEG8000, 0.8 U/µl RNaseOut, 0.8 U/µl SUPERase•In™ RNase inhibitor, around 200 ng mRNA (around 70 nM for 200 nt average size) and 64 µM enzyme in 40 µL scale reaction for 2 h at 37 °C. Cellular RNAs were purified using the RNA Clean and Concentrator kit (Zymo Research) and stored in RNase-free H2O.

Dot blot assay
RNA samples after enzymatic labeling and click reaction were in 2 µL water, boiled at 93 ˚C for 3 min to denature local structure, and then promptly put on ice to avoid reformation of structure. Dotted the desired amount of RNA from each sample on the Amersham hybond N+ membrane, with 1 µL each round and wait until totally dry before dot again. Crosslinked in auto-crosslinker twice under autocrosslinking mode (70,0000 microJ /cm2, 254 nm). Fluorescence signals were directly read with a bioimager. Chemiluminescence signals were detected by the chemiluminescence nucleic acids detection kit. For antibody tests, RNAs without labeling were denatured, dotted, and crosslinked. Then western blots were performed following the manual of the antibodies.

MALDI-TOF
Purified RNA oligonucleotides from the enzymatic labeling reaction were resuspended in 5 μL water, mixed with 40 μL ion exchange resin (Bio-Rad), and left at room temperature for 30 min to remove salts. We loaded 1 μL of supernatant of the mixture to a MALDI sample plate, mixed with 1 μL THAP matrix freshly added with salt, and pipetted to mix until crystals were formed. The data was acquired using a Bruker Ultraflextreme MALDI-TOF-TOF mass spectrometer using a negative reflector mode and analyzed using the Bruker FlexAnalysis 3.4 software. More than two biological replicates were analyzed for each sample type.

Mammalian cell culture and RNA preparations
HEK 293T cells (purchased from ATCC (cat# CRL-3216)) and HeLa cell lines (purchased from ATCC (cat# CCL-2)) were grown and maintained in DMEM supplemented with 10% FBS and 1% penicillin-streptomycin at 37 °C with 5% CO2. Total RNA was extracted from cells using Trizol reagent and poly-A RNA was extracted using Dynabeads® mRNA DIRECT™ Kit for two rounds with an additional ribosomal RNA removal step using the RiboMinus™ Eukaryote Kit v2.

LC-MS/MS
RNA was digested with 3 μl 10× Reaction Buffer (100 mM NH4OAc, pH = 5.3) and 1 μl Nuclease P1 (>1U/ μl) in 25 μl solution at 42 °C for 2 h to break the phosphodiester bond, after which phosphate was removed by adding 1 μl FastAP (1U/μl) (ThermoFisher, EF0654) and 3 μl 1M NH4HCO3 (pH unadjusted), and incubating at 37 °C for 2 h. Upon the completion of incubation, we diluted the reaction mixture to 50 μL with water and filtered the samples through 0.22 μm Millex-GV PVDF filters (Millipore). Five microliters of sample were injected into a ZORBAX SB-Aq 4.6 x 50 mm column (Agilent) on UHPLC (Agilent) coupled to a SCIEX 6500+ Triple Quadrupole Mass Spectrometer in a positive electrospray ionization mode. The nucleosides were quantified based on the nucleoside to base transitions: 268.0 to 136.0 (A), 284.0 to 152.1 (U), 245.1 to 125.0 (Ψ), 259.1 to 139.0 (m1Ψ) and 283.1 to 163.0 (p1Ψ). At least two independent biological replicates were used for Ψ level quantification, and each sample was injected 2 times.

Primer extension assay
In this assay, a 33-mer synthetic RNA with a Ψ modification in the middle was used as the template. A fluorescently labeled RT primer was designed to target the position immediately downstream of the modification site. The RT primer and template were mixed at a 1:10 ratio to drive the majority of primers to hybridize with the template. In the absence of any treatment, the RT primer reads through the modification site, producing the full-length RT product (upper FL band). Conversely, when encountering the propargyl-Ψ adduct (with or without subsequent biotinylation at a low dNTP concentration), the elongation of the RT primer is impeded, resulting in a significant amount of un-extended RT primer (lower Primer band in Figure 2e). The limited presence of un-extended RT primer bands in untreated RNAs is likely attributed to a low level of unhybridized RT primer.

ELAP-seq
RNA fragmentation, enzymatic labeling, end repair and click reaction
RNA was first fragmented using 10x fragmentation buffer (Thermos Fisher) in a 10 µL scale reaction containing 1x fragmentation buffer and 200-1000 ng mRNA. The fragmentation reaction was incubated at 70 °C for 12 min to fragment the RNA to 60-200 nt. RNA was recovered using the Oligo Clean and Concentrator kit (Zymo Research) and eluted in water. For downstream library preparation, the same amounts of mRNA (200 ng) were used in library construction processes. The recovered RNA was denatured at 90 °C for 3 mins and quickly set on ice to break RNA secondary structures. One tenth of each RNA sample was aliquoted as input.

The remaining RNA was then labeled in an in vitro methylation reaction using pro-SeAm as co-factor following the protocol mentioned above. RNA was purified using the RNA Clean and Concentrator kit (Zymo Research), eluted in 9 µL of water, and all of it was used in a 50 µL scale end repair reaction containing 1x T4 PNK buffer (NEB B0201S) and 5 µL T4 PNK (Thermos Fisher) with incubation at 37 °C for 1h. The RNA was then purified using the Oligo Clean and Concentrator kit and eluted in 15 µL water.

We used half of the samples in a 100 µL scale click reaction containing 250 µM biotin-azide, 0.4 mM CuSO4, 4mM THPTA and 4 mM sodium ascorbate which was added last. We used the other half in a control reaction where THPTA and CuSO4 were omitted. We incubated the reactions at room temperature with shaking at 600 rpm for 1.5 h. RNA was recovered by ethanol precipitation at -80 °C overnight. The input RNA was treated the same way with the end repair reaction above and used directly in RNA ligation.

RNA ligation, pull-down and reverse transcription
RNA recovered from ethanol precipitation was resuspended in 9 µL water and mixed with 1 µL of 30 µM 3’ RNA linker (5′rApp-NNNNNAGATCGGAAGAGCGTCGTG-3SpC3). The mixture was denatured at 72 °C for 2 mins and then quickly set on ice. Added to the mixture T4 ligation buffer to 1X, PEG8000 to 15%, 0.8 U SUPERase• In™ RNase Inhibitor, and water to a total of 23 µL. Mixed thoroughly, added two µL T4 RNA ligase 2 truncated KQ and mixed again. Incubated the mixture at 25 °C for 2 h and then 16 °C for 12 h. RNA was then purified with the RNA Clean and Concentrator kit.

Aliquoted MyOne C1 magnetic beads suspension (5 µL per sample) and washed the beads following the manual. Blocked the beads with 2 volumes of 1x W&B buffer containing 0.5 mg/ml BSA and 50 µg/ml tRNA for 1 h at 4 °C. Resuspended the beads in 2 volumes of pull-down buffer containing 10 mM Tris-HCl (pH 7.5), 2 M NaCl, 1 mM EDTA, 8.8 µM biotin-azide, 44 µg/mL yeast tRNA, 0.44 mg/mL BSA and 0.1% Tween-20. When working with beads, use low adhesion tubes and tips to minimize loss of beads due to sticking to the tube wall or tips.

Mixed 17 µL purified ligation product with 10 µL bead suspension, and added 15 µL of 0.73x IP buffer (without biotin-azide) to decrease NaCl concentration to 1 M to facilitate biotin-streptavidin binding and decrease PEG concentration to below 5% to avoid bead aggregation. Beads aggregation problem, if met, could be mitigated by adding Tween-20 or DTT when applicable. Incubated the mixture for 1 h at 4°C and then 30 min at room temperature. Wash the beads 6 times with 100 µL 1x W&B buffer supplemented with 0.05% Tween-20, and then 2 times with 50 µL 1x W&B buffer with each wash incubated at 4°C for 5 mins. Changed tubes every three washes. Changed tubes again after the final wash and resuspended beads in 12 µL 10 mM Tris-HCl buffer (pH 7.5). Beads after pull-down were divided into two parts for the downstream reverse transcription step, with three-quarters (corresponding to 150 ng starting RNA) used for SuperScript III RT and one-quarter (corresponding to 50 ng starting RNA) for SuperScript IV RT.

For input RNA samples, diluted each ligation reaction to 47 µL with water. To remove un-ligated adapters, first add 2 µL 5’ de-adenylase and incubate at 30 °C for 1 h, then add 1 µL RecJ and incubate at 37 °C for 1 h. Purified the RNA using RNA Clean and Concentrator kit and eluted the product in 12 µL water.

Beads resuspended in six µL Tris-HCl buffer from the pull-down experiment or 6 µL purified RNA from input was used in reverse transcription (RT). Set up RT reaction by adding 1 µL 2 mM RT primer (5’-ACACGACGCTCTTCCGATCT-3’), 1 µL 20 µM dNTP and water to 13 µL. For Superscript IV, used 10 mM dNTP instead of 20 µM dNTP. Denatured the mixture in a shaking incubator at 65 °C for 5 min with shaking at 900 rpm and then quickly set on ice for at least 1 min. Added 4 µL 5x SSIV buffer (for the Superscript IV) or 5x first strand synthesis buffer (for the Superscript III), 1 µL 100 mM DTT, 1 µL RNaseOut, and 1 µL Superscript IV or Superscript III reverse transcriptase, and incubated at 50-55 °C for 10 min (for Superscript IV) or 50 min (for Superscript III) with shaking at 900 rpm. After the RT reactions, heat denatured the RT proteins. Then washed beads once with 100 µL high salt buffer composed of 20 mM Tris-HCl pH 7.4, 1 M NaCl and 0.1% tween-20, and once with 50 µL low salt buffer composed of 20 mM Tris-HCl pH 7.4, 100 mM NaCl and 0.1% tween-20. Resuspend beads in 18 µL buffer composed of 50 mM Tris-HCl pH 8.3, 75 mM KCl, 3 mM MgCl2 and 0.1% tween-20 and denatured the RT enzymes by heating at 80 °C for 15 mins. After cooling the reaction on ice, added 1 µL 200 mM DTT and 1 µL RNase H and digested at 37 °C for 30 min with shaking at 900 rpm. Heated the mixture at 72 °C for 2 min and quickly set on the magnet to collect the supernatant. cDNA was purified from the supernatant using Oligo Clean and Concentrator and eluted in 8 µL water.

cDNA ligation and PCR
All cDNA was used to mix with 1x T4 RNA ligase buffer, 0.81 µM adapter (5′Phos-NNNNNAGATCGGAAGAGCACACGTCTG-3SpC3), 1 mM ATP, 24% PEG8000, 7.5% DMSO and water in a 22 µL scale reaction, and incubated at 72 °C for 2 min before quickly chilling on ice for 2 min. One µL of 50 mM Co(NH3)6Cl3 and 1.8 µL T4 RNA ligase 1 (high concentration) were added and the ligation reaction was allowed to proceed at 25 °C for 12 h, followed by protein denaturation at 65°C for 5 min. The ligation product was purified with 1.9X Ampure beads and eluted in 10 µL water.

Two µL ligation product was used in qPCR to estimate the number of cycles needed in library amplification PCR. The rest of ligation product was amplified in a 25 µL scale PCR reaction containing LongAmp Taq 1x mix and 0.24 mM forward and reverse primers. To minimize adapter dimers, PCR products of the pull-down samples reverse transcribed using Superscript III were first purified with 0.9X Ampure beads, further amplified for 6 cycles, and then loaded on a 3% low melting-point agarose gel to remove primer dimers. For input samples and pull-down samples built from superscript IV, PCR products were directly run on a 3% low-melting point agarose gel. Control libraries are expected to yield no DNA product under the same PCR cycle numbers as the experimental group. Excise the bands in the desired range and perform gel extraction to recover pure DNA for Illumina sequencing PE150. 25M reads were obtained from each sample.

Library preparation for synthetic RNA oligos containing 5’-NNΨNN-3’ and 5’-NNUNN-3’ sequence
Library preparation for synthetic RNA oligo was performed as above with the following modifications. Synthetic RNA (20 ng) was used in in vitro methylation reaction using pro-SeAm as co-factor following the protocol mentioned above. RNA was purified using ethanol precipitation, resuspended in 9 µL of water, and entirely used in a 50 µL scale end repair reaction for incubation at 37 °C for 1h. The RNA was then purified using ethanol precipitation and resuspended in 9 µL water. One microliter RNA was aliquoted as input.

Both input RNA and RNA for IP were end-repaired and ligated with RNA adapter. After on-bead pull down, reverse transcription was performed on-bead using superscript IV with 0.5 mM dNTP. The on-bead reaction and wash substantially removed adapter-dimers from the libraries.

ELAP-seq Data analysis
Synthetic standard oligo data analysis
Paired end reads were merged using bbmap and deduplicated. 5nt and 6nt were trimmed respectively on the 5’ and 3’ ends to remove UMI with the option of -q 10, 10 to further trim both ends with quality cutoff of 10. Reads of synthetic oligonucleotides were extracted and sorted using the command grep -o "GCTCGCAG.*A" data.fq.gz > output.txt <output.txt sort | uniq -c > output1.txt . Subtotals of number of reads for sequences with identities varying at 2 nucleotides upstream and downstream of the Ψ site were summarized. Enrichment level was determined by dividing the RPM in the pull-down sample to the RPM in the input sample. The heatmap for enrichment levels of all 256 sequences was plotted using R.

rRNA data analysis
Reads R2 are used for analyzing data of biological RNA. Poly-A RNA ELAP-seq data also contain rRNA reads thus was used for analyzing Ψ modifications in rRNA. Adapter sequence of “AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT” was first trimmed from the 3’ end of library R2 of HEK poly-A built with superscript III, followed by de-duplication and removing UMI of 6 nt and 5 nt on the 5’ and 3’ end respectively. By design, after trimming adapter and UMI sequences, the first nucleotide of R2 is the stop site of reverse transcription. The trimmed reads were mapped to ribosome RNA genome using hisat2 with the options of -k 1 --no-softclip --mp 1000,999 to allow only one distinct alignment, disallow soft-clipping and allow no mismatch. RT arrest was called using JACUSA2 with the options of ‘rt-arrest -m 0 -p 2 -c 1 -P FR_SECONDSTRAND’ to filter off positions with mapping quality < 0 or reads coverage < 1 and specify second strand sequenced. The resulting output file was processed sequentially with in house scripts rRNA_filter_p1.sh, Stutter_removal_rRNA.py and rRNA_filter_p2.sh to calculate stop ratios, remove false positive sites due to the stutter effect, and select for sites whose stop ratio is > 0.05 in the pull-down sample. In Stutter_removal_rRNA.py, we required that the stopped reads at a candidate Ψ site must be at least 15% higher than its neighboring sites. Sites called from two replicates were overlapped for stop patterns between independent biological replicates. The resulting sites were further filtered in excel for final candidate sites.

A candidate site in rRNA is defined when meeting the following criteria: 1) the stop ratio at the site in the pull-down library is > 0.05 in both replicates; 2) stop ratio in the pull-down library – stop ratio in the input library is > 0.05 in both replicates; 3) the number of stopped reads in the pull-down library is >=5 in both replicates; 4) the average value of stop ratio * stopped reads between the two pull-down replicates is >=1.5 (this imposes a more stringent requirement on the number of stopped reads when the stop ratio in the pull-down sample is lower than 30%); and 5) the stop event occurs at uridine (U). Finally, we further remove sites meeting the following criteria: 1) the stop ratio in the input library is > 0.1 and > 1/3 of the stop ratio in the pull-down library, and 2) number of stopped reads in input is >=3 in any replicate. When using superscript III and IV data together, we follow the same criteria except for requiring that the stop ratio in the pull-down library- stop ratio in the input library is > 0.1.

Poly-A RNA data analysis
Reads R2 from samples generated with superscript III were trimmed 6 nt on the 5’ end and 5 nt on the 3’ end, whereas reads from samples generated with superscript IV was trimmed 7 nt on the 5’ end and 5 nt on the 3’ end as superscript IV tends to elongate for one more nucleotide passing Ψ. Reads after adapter trimming, deduplication and UMI trimming were mapped to hg38 using hisat2 with options of “--known-splicesite-infile --rna-strandness F --no-softclip”. To combine samples built from superscript III and IV, bam files for the same sample generated from Superscript III and IV were combined to generate a new bam file. We then used MACS2 to call highly enriched peaks and categorized all mapped sites in the pull-down samples into groups inside (in-group) and outside (out-group) of highly enriched peaks. We used JACUSA2.0 to calculate RT-arrest at each site, with parameters ‘rt-arrest -m 0 -p 2 -c 1 -P FR_SECONDSTRAND’.

We pre-processed the resulting files with the following steps: 1) we filtered away sites in any region covered by reads that all share exactly the same start and end locations, which is likely caused by off-target mapping of reads whose originality is another region in the genome containing a same sequence as the current region; 2) we removed RT stop sites within 50 nt upstream (when reads are mapped to the - strand) or downstream (when reads are mapped to the + strand) of a main RT stop site to decrease false positives due to the effect of diminished processivity of the RT enzyme reading beyond the Ψ modification; and 3) to remove the noise caused by off-target mapping to repetitive regions, we only select for sites that are covered by at least one uniquely mapped read for further analysis.

Next, we determined candidate Ψ sites based on the same principle as used for detecting modifications in rRNA as described above. We raise the stop ratio in pull-down to be higher than 0.1 to decrease potential false positives, considering less crowdedness of modification sites in mRNA compared to rRNA. Sites identified from superscript III alone and new sites identified from combined reads of Superscript III and IV were combined.

Metagene analysis was performed using MetaPlotR48. Gene ontology enrichment analysis was done using clusterProfiler49.

IVT RNA preparation and processing
The protocol reported by Xiao et al. was followed41, with the exception that RNase A/T1 digestion was performed in addition to RNaseH digestion after reverse transcription to remove any residual RNA template. The library construction and data analysis followed the same procedures as those used for the biological RNAs.

Validation using CMC-RT-qPCR
Total RNA was treated with Cyclohexyl-N’-(2-morpholinoethyl)carbodiimide (CMC) as previously described34. Briefly, 15 or 10 µg total RNA in 12 µl sterile H2O was used for the treatment (CMC +) and the control reaction (CMC-). 24 µl TEU buffer (50 mM Tris-HCl (pH 8.3), 4 mM EDTA, 7 M urea,) and 4 µl freshly prepared 1M CMC (C106402-5G, Sigma) in TEU buffer or 4 µl TEU buffer was added to the CMC+ or CMC- tube respectively. Both samples were incubated at 30 ºC for 16 h to ensure complete conversion of Ψ to Ψ-CMC adduct. RNA was then purified with ethanol precipitation. The RNA pellet was then air dried and resuspended in 40 µl of 50 mM sodium carbonate and 2 mM EDTA (pH 10.4) and incubated at 37 °C for 6 h to remove other CMC adducts. RNA was purified again with ethanol precipitation, air dried and resuspended in sterile deionized H2O and stored at -80 ºC before use. The RNAs were used in RT-qPCR using primers as listed in Supplementary Data 7.

Estimation of modification levels
For quantification of modification levels, we used enrichment levels at the modification sites. Enrichment level is defined as the local read coverage per million mapped reads (RPM) value at the site in the pull-down sample divided by that in the input sample. For calculating enrichment levels, two input samples can be combined into one to increase read coverage and to minimize technical variation. For more reliable quantification, we suggest the readers apply a cutoff of 1 for Input-RPM, which corresponds to at least 10 local read coverage under the current sequencing depth (10-20 M mapped reads). Such cutoff can be user-defined based on their sequencing depth and research needs.

We used a synthetic oligonucleotide with sequence of 5'-GGA GAC GGU CGG GΨUC CAG AUA UUC GUA UCU GUC to study the relationship between the enrichment level and the modification percentage. Different ratios of biotin-propargyl-Ψ -modified and unmodified oligos are mixed to make a total of 0.1 pmol comprising 1%, 9%, 25%, 50% and 100% of modified oligos. Pull-down was performed followed by RT-qPCR to determine relative amount of pull-down oligos. By normalizing the enrichment level to the data of modification of 1%, an exponential correlation was found between the modification percentage and the normalized enrichment level: Y = 1.879 e0.0469X , where Y is the normalized enrichment level and X is the apparent modification level. Based on this result, we propose to estimate the relative modification levels of each modification site using the natural logarithm of the enrichment level at the site adjusted by the preference for the specific sequence context by Mj1640 as described below.

Relative modification level= ln (Enrichment at the modification site / Enrichment of the oligo in the sequence context)

Based on the correlations in Figure 4j (for sites in non-ΨU contexts) and Supplementary Figure 15e (for sites in ΨU contexts), sites are defined as highly modified when relative modification levels are >= 2.3 for sites in ΨU context and >= 3.5 for sites in non-ΨU context. Sites are defined as moderately modified when relative modification levels are >= 1.6 for sites in ΨU context and >= 1.8 for sites in non-ΨU context. Sites are defined as lowly modified when relative modification levels are < 1.6 for sites in ΨU context and < 1.8 for sites in non-ΨU context.

DKC1 knockdown
The day before transfection, 2 × 10⁶ HEK293T cells were seeded onto each 15 cm cell culture plate. Transfection was performed using RNAiMax transfection reagent following manufacturer's instructions. For one 15 cm plate, 20 μL of 20 μM siDKC1 (Qiagen, Hs_DKC1_6 FlexiTube siRNA) or siCtrl (Sigma Aldrich, catalog # sic001) was used. Cells were harvested 72–96 hours post-transfection. Library construction was performed as previously described, with the following modifications. To minimize variability during the pull-down and downstream processing steps, we designed four RNA adapters containing unique 3-letter internal barcodes (Supplementary Data 8). These adapters were individually ligated to the two siCtrl and two siDKC1 samples, respectively. Following the RNA ligation step, the four input or enzymatically labeled samples were pooled together to ensure uniform handling in subsequent steps.

Paired-end sequencing reads (R1 and R2) were first merged, and samples were then demultiplexed based on their internal barcodes. After deduplication to remove PCR duplicates, the internal barcode and a UMI (8 nt) were trimmed from the 5′ end, and a second UMI (5 nt) was trimmed from the 3′ end. We used the known modification sites identified from HEK293T cell line as a reference and local read coverage at these sites and on the expected strand (noting strand reversal when applying results from R2-only mapping to the analysis using merged-read mapping) was obtained using JACUSA2 for downstream analysis.

Differential modification was evaluated using changes in the enrichment level defined by local RPM at the site in the pull-down normalized by that in the input. We first selected sites whose enrichment fold change is greater than 2 with p-value < 0.05. For downregulated sites, we further excluded sites with average pull-down RPM values < 1 before treatment (corresponding to approximately 10× local read coverage given the current total mapped reads of 6-10M), as their enrichment efficiencies are more susceptible to technical variation. For more reliable fold-change quantification, readers can apply a cutoff of 1 for Input-RPM under both siCtrl and siDKC1 conditions, which corresponds to at least 10x local read coverage given the current total mapped reads of 10-20 M. Such cutoff can be user-defined based on their sequencing depth and research needs. A more detailed analysis procedure can be found in the GitHub repository associated with this paper.

Differential gene expression analysis was performed on RNA-seq data of input samples for both siCtrl and siDKC1 treated cells using DESeq2.

RNA structural analysis
RNA sequences of 100 nt surrounding all Ψ sites located in different mRNA subregions (5’ UTR, 3’ UTR and CDS) with 50 nt upstream and 50 nt downstream were used for RNA secondary structure prediction. To achieve this, the program viennaRNA was used with the command of ‘RNAfold -p < input.fa > output.res’. The resulting file was cleaned up to keep only the information of the RNA sequence and the most confident prediction (the one with lowest ΔG values) and saved as an excel file. Next, an in-house python script was used to extract the information of the base pairing status of each modification site as well as the location of the matched base if it is base paired. The identity of the matched base was further determined in the excel using the ‘MID’ function. Next, we calculated the frequency for all modification sites to reside in a mismatch/UG wobble pair versus in UA base pair. Control U sites were chosen at both the -100 position and +100 position of Ψ sites and analysis were performed the same way as for Ψ sites. mfold were used to predict secondary structures of 100 nt surrounding the chosen Ψ sites.

Statistics and reproducibility
Statistical analysis for siDKC1 knockdown data used one-sided student t-test and p-values smaller than 0.05 were considered significant. For DKC1-dependent and -independent sites determined by PRAISE, the statistical analysis of these sites in our dataset was performed using a two-sided Wilcoxon rank-sum test. In this paper, error bars represent mean ± standard error of the mean (s.e.m), n>=3. No statistical method was used to predetermine sample size, and no data were excluded from the analysis.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Pseudouridine (Ψ) is an abundant and critical RNA modification, yet the full extent of its biological roles in the cell have been difficult to understand due to incomplete transcriptome-wide maps of modification sites. One of the major bottlenecks in generating comprehensive maps with current mapping methods is the extensive sequence depth required to confidently assign modifications globally, due the sparse distribution and general low occupancy of Ψ. In this manuscript, Wang et al., present an innovative method to enrich for RNA fragments containing Ψ by using an archaeal methyltransferase, which they claim can specifically add click-chemistry compatible propargyl groups to Ψ to facilitate biotin conjugation. This addition allows them to both enrich for modified RNA using streptavidin, and obtain single nucleotide resolution of modification position using RT-stop signatures induced by the biotin adduct. They apply this method to correctly identify most known rRNA pseudouridine sites and generate a map of pseudouridylation in human mRNA across three cell types and make several claims about global Ψ distribution.\nWhile we find this technique to be a potentially powerful new tool for profiling Ψ, the lack of overlap with existing methods and little validation for the many new sites uncovered is cause for concern. Additionally, several of the claims made about Ψ occupancy and distribution need more supporting evidence or explanation. We do not recommend acceptance of this manuscript until the following points are addressed:\nMajor Concerns:\n1. According to figure S11, only 4% of sites identified in HeLa cells overlap with previous mapping datasets, and 8% overlap in 293T. The authors do not address this discrepancy beyond this figure, except to suggest it is due to poor performance of CMC-based methods. More evidence is needed to demonstrate that the mRNA data is capturing mostly genuine Ψ modifications, including performing PUS enzyme depletions to demonstrate PUS-dependence of signal, and far more targeted validations (ex. using CLAP), specifically of sites not recapitulated by other datasets.\n2. The authors claim that 1/3 of sites have over a 50% modification rate, much higher than the average mRNA modification rates estimated from other methods such as BID-Seq, and inconsistent with their claim that the massive number of new sites are from enrichment of lowly-modified positions. The authors must more rigorously demonstrate that their estimation of modification occupancy is accurate, by correlating their estimates with estimates from mutational-profiling methods such as BID-Seq, and by demonstrating correlation with many more than 2 sites estimated by CLAP. Also, the CLAP vs ELAP-Seq comparison is difficult to see in figure 4h and should have its own figure panel (with the additional comparisons).\n3. The language explaining the rational for cutoffs and filtering is vague. Please provide specific details with statistical justification for the following:\n“We optimized cutoffs for stop ratio and reads coverage and allow for gradients of coverage for different ranges of stop ratios” [275-276]\n“We required higher reads coverage for lower levels of stop ratios to decrease false positives” [685-686]\n“We allowed slightly lower reads coverage for sites in the out-group” [688-689]\n“…we screened cutoffs for stop ratios and coverage to find a balance between the number of sites identified and percentage of U” [689-691]\n4. For the IFN dataset it is mentioned that “the interferon stimulated samples (IFN+) resulted in ~40% less reads compared to untreated samples (IFN-), thus we allowed lower reads coverage cutoff for the IFN+ samples.”[715-719] How does this not increase the false positive rate for the INF+ samples – if the read cutoffs were set to minimize false positive rate? Intuitively this cutoff should not be determined by read depth. This comparison would be more accurate with the INF- samples randomly sub-sampled down to the same depth as INF+, using identical read cutoffs.\nMinor critiques:\n“a significant proportion of Ψ sites in mRNA appear in the 5’-ΨU-3’ context…” [203-204] what proportion exactly?\n“30 – 60% of Ψ in cellular mRNA can be labeled by this enzyme” [205-206] Do you expect 30-60% of Ψ positions are being labeled – or each site is being labeled at a 30-60% rate?\nWhy is there a stop in untreated RNA at the expected psi position in figure 2e?\nHow do the NNpsiNN pools perform in the standard library prep? E.g. are there any RT-stop or ligation biases due to sequence context?\nIt is not described how the wobble pair analysis was performed, please add more detail in the methods or main text.\nEnrichment of Ψ at splice junctions is only informative with a comparison against the background distribution of U. Plot both distributions from your dataset and demonstrate that the difference in Ψ distribution is significant.", "reviewer_3": "The manuscript by Wang et al developed ELAP-seq to enrich and map Ψ at single-base resolution. This new Ψ sequencing technology leveraged enzymatic strategy to label and enrich Ψ-modified transcripts. They have identified approximately 10, 000 candidate Ψ sites in HEK293T and HeLa cell lines, a significantly larger number compared to all previous CMC-based and bisulfite-based methods. Based on this approach, the authors found differential Ψ modification upon type-I interferon stimulation in the glioblastoma cell line. Overall, the manuscript provides a new sequencing method and applies it to detect Ψ modification in the human transcriptome. Unfortunately, the major concerns remain non-specificity, unclear false positives, and unclear false negatives. Some weaknesses need to be addressed before this manuscript would be suitable for publication:\nMajor:\n1. The authors reported around ten thousand candidate Ψ sites from HeLa and HEK 293T cell lines. However, more than 9000 sites currently lack validation through existing datasets. Before asserting their claims of high sensitivity, the authors have to verify the specificity of ELAP-seq instead of just providing the ambiguous Ψ lists. The authors should perform perturbation of PUS proteins, such as employing PUS knockdowns or knockouts in these cells, to assign the majority of the newly identified Ψ sites to specific enzymes. This assignment would validate the specificity of ELAP-seq and contribute to functional insights.\n2. A more comprehensive comparison should be performed between the Ψ dataset sourced from HEK293T and HeLa cells in ELAP-seq and their previously published human mRNA pseudouridine maps (Dai, Q. et al. Nat Biotechnol. 2023). They should analyze the raw deletion signal for ELAP-seq-unique Ψ sites in BID-seq classify the authenticity of these sites and provide a clear explanation for their lack of detection in BID-seq.\n3. In L187, the strong bias for ΨU context might even underestimate the GUΨC context modified by TRUB1 protein which has the highest modification activity, which indicates the severe false negative of ELAP-seq. The authors need to demonstrate they could detect all these TRUB1-dependent sites from previous Ψ datasets.\n4. The authors claim the high sensitivity and low background of ELAP-seq; however, they have yet to provide substantial supporting evidence for this assertion. (1) For Ψ in rRNA, they only identified 70 out of 104 known Ψ sites in rRNA (L250), since Mj1640 the strong preference for 5’-ΨU-3’ sequence context in oligo RNA, they should analyze these 34 missed known Ψ sites, and then address whether these omissions were a result of the sequence context instead of being clustered within specific rRNA regions. (2) For Ψ in mRNA, although the author can filter out these highly modified Nm sites through bioinformatic cutoff, they cannot exclude the existence of false-positive Ψ sites derived from Nm modifications that vary with differential modification levels in mRNA, hence, the reliability of ten thousand of Ψ is questionable. The comparison and overlap analysis of Ψ sites identified via ELAP-seq with Nm sites identified through Nm-seq should be performed (Dai, Q., et al Nat Methods, 2017).\n5. In L267 and Figure 4a, how does the mitochondrial mRNA Ψ result compare to the previous Ψ modification map obtained by BID-seq? The comparison of mitochondrial Ψ sites should be performed for both HEK293T and HeLa cell lines.\n6. In L328, the authors claim “excellent” repeatability for Ψ detection in mRNA. They should provide mRNA Ψ site numbers for three biological replicates and hence the ratio of overlapping sites.\n7. In L432, the authors need to provide more experimental and bioinformatic analysis for 92 downregulated Ψ sites upon interferon-α treatment, for instance, they can analyze the motif and investigate whether they were modified by PUS3 protein. The authors could also detect the protein level of relevant PUS enzymes between IFN- and IFN+ samples to explain why Ψ was downregulated.\nMinor:\n1. In L55, the statement regarding quantitative mapping methods is incorrect, RBS-seq is not a quantitative approach.\n2. In L130, data was not shown for the control U-RNA.\n3. In L146, the authors should provide experimental data to demonstrate the capacity of the Mj1640 enzyme to label the synthetic 8-mer oligo using propargyl-SeAm.\n4. In L158, the size distribution of untreated mRNA should be provided."}, "author_response": "We thank the reviewers for the detailed analysis and constructive comments of our manuscript.\nBelow, we have positively addressed reviewers’ concerns regarding the authenticity, the\nreproducibility, and the quantitative capability of the ELAP-seq method by performing additional\nexperiments, conducting thorough data analysis, as well as re-writing part of the manuscript. We\nhave also improved the data analysis pipeline and provided a one-step command. We appreciate\nthe suggestions from the reviewers to perform additional work which greatly enhanced the\nclarity and usability of our method.\nREVIEWER COMMENTS\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nPseudouridine (Ψ) is an abundant and critical RNA modification, yet the full extent of its\nbiological roles in the cell have been difficult to understand due to incomplete transcriptome-\nwide maps of modification sites. One of the major bottlenecks in generating comprehensive\nmaps with current mapping methods is the extensive sequence depth required to confidently\nassign modifications globally, due the sparse distribution and general low occupancy of Ψ. In\nthis manuscript, Wang et al., present an innovative method to enrich for RNA fragments\ncontaining Ψ by using an archaeal methyltransferase, which they claim can specifically add\nclick-chemistry compatible propargyl groups to Ψ to facilitate biotin conjugation. This addition\nallows them to both enrich for modified RNA using streptavidin, and obtain single nucleotide\nresolution of modification position using RT-stop signatures induced by the biotin adduct. They\napply this method to correctly identify most known rRNA pseudouridine sites and generate a\nmap of pseudouridylation in human mRNA across three cell types and make several claims about\nglobal Ψ distribution.\nWhile we find this technique to be a potentially powerful new tool for profiling Ψ, the lack of\noverlap with existing methods and little validation for the many new sites uncovered is cause for\nconcern. Additionally, several of the claims made about Ψ occupancy and distribution need more\nsupporting evidence or explanation. We do not recommend acceptance of this manuscript until\nthe following points are addressed:\nMajor Concerns:\n1. According to figure S11, only 4% of sites identified in HeLa cells overlap with previous\nmapping datasets, and 8% overlap in 293T. The authors do not address this discrepancy beyond\nthis figure, except to suggest it is due to poor performance of CMC-based methods.\nMore evidence is needed to demonstrate that the mRNA data is capturing mostly genuine Ψ\nmodifications, including performing PUS enzyme depletions to demonstrate PUS-dependence of\nsignal, and far more targeted validations (ex. using CLAP), specifically of sites not recapitulated\nby other datasets.\nOur response:\nWe are grateful for the constructive feedback from reviewer #1. The reviewer has noted a 4%\nand 8% overlap between our data and previously reported sites for HeLa and 293T mRNA,\nrespectively. The lower percentage of overlap is attributed to a greater number of sites identified\nby our method due to its higher sensitivity, which we view as an advantageous feature.\nIn this revised version of manuscript, we performed a comparative analysis with the three most\nrecently published bisulfite-based methods (BID-seq, RBS-seq, and PRAISE). Our method\ndemonstrated substantial coverage, and successfully identified 42%, 49%, and 60% of the\nreported sites, respectively (Figure S10a). Overall, our data encompassed 1,904 of the 4,333\npreviously reported sites for 293T mRNA (combined data from BID-seq, PRAISE and CeU-seq)\n(Figure S10a). Our dataset not only validated these modification sites but also facilitated the\nprecise localization of Ψ modifications within multiple Us, a challenging task for the current\nbisulfite-based approaches.\nTo address discrepancies with previous methods, we further looked at signals of PRAISE-unique\nsites in our dataset. Our data uncovered 1,237 of the 2,070 sites reported by PRAISE (Figure S10\na-b). Among the 833 PRAISE-unique sites, 526 exhibited prominent stop signatures in at least\none replicate of our data, with 113 showing signals in two replicates. However, these sites were\nexcluded from our final list due to being present in only one replicate or having low read\ncoverage (Figure S10e). Only 288 sites (13.9%) lacked any signal in our dataset, and these sites\nwere found to have lower deletion ratios and read coverage than the rest Ψ sites in the PRAISE-\nseq dataset (Figure S10f), suggesting that they are intrinsically more difficult to identify and may\nhave lower confidence. Our examination also identified 11 sites determined by PRAISE that\nwere covered by multi-mapped reads, potentially leading to repeated counts, and 8 sites flagged\nas false positives due to identical stop signatures in both the pull-down and input in our data.\nSimilarly, for the 302 BID-unique sites, we identified 65 additional sites showing prominent stop\nsignals in at least one replicate of our dataset, with those lacking signals in our dataset mainly\nsituated in unfavorable sequence contexts and displaying low deletion levels in BID-seq dataset.\nFurthermore, our data covered 70% of TruB1 sites reported by PRAISE (Figure S10h) and all\nmitochondria mRNA sites reported by PRAISE and BID-seq (Supplementary table 3-4).\nNotably, by comparing the U sites and non-U sites identified by our method, we found that\nincreasing cutoff for stop ratio could increase the degree of confidence for the reported sites\n(Figure 4 g-h). Additionally, by performing a third replicate for each cell line, we found that the\ndegree of confidence increases with the number of replicates (Figure 4i). This suggests a means\nto better characterize the authenticity of each candidate site. Based on these observations, we\nmore stringently defined the degree of confidence for each candidate Ψ site based on the stop\nratio and the number of replicates the site can be repetitively identified in. Among the 10106\nsites identified from HeLa cells, 4260 are defined as “highest-confidence” sites as they fulfil\nboth criteria of being identified in all three biological replicates and having stop ratios higher\nthan 0.3, 3801 sites are defined as “higher-confidence” sites for fulfilling one of the criteria, and\n2044 are defined as “lower-confidence sites” for fulfilling none of the criteria (Figure 4j and\nLines 329-336).\nThe assignment of confidence, coupled with the enhanced accuracy in estimating the\nmodification levels as addressed in the response to question #2 from reviewer #1, provides\nvaluable insights for characterizing each individual candidate site.\nUpon inspecting raw signals in BID-seq for 6,333 ELAP-seq-only higher or highest-confidence\nsites, we found that 3,012 sites showed deletion ratio of higher than 0.05 in at least one replicate\nof BID-seq, thus these sites demonstrate modification signals in BID-seq as well (Figure S10j).\nAmong the rest sites, 2,139 sites have deletion ratio of between 0.0 and 0.05 in BID-seq with\n82% of them having average deletion ratios among three replicates of between 0.02 to 0.05,\nwhich is below the deletion ratio cutoff used by the method (Figure S10j). There are 1,182 sites\nhaving deletion ratio of 0 in BID-seq data, and these sites exhibited substantially lower read\ncoverage than sites detected by BID-seq, suggesting that these sites were missed by BID-seq due\nto low gene expression (Figure S10k). This further supports the authenticity of sites identified by\nELAP-seq and underscores the high sensitivity of our approach in validating weak signals in\nBID-seq.\nIn addition to more thorough data analysis, we also performed CLAP and CMC-RT-qPCR\nvalidation for more sites. In addition to the ten sites presented in the previous version, we\nevaluated 1 more site using CLAP (Figure S9c) and 7 more sites using CMC-RT-qPCR (Figure\nS7f). Overall, 18 sites were successfully validated, including 10 sites not previously identified by\nany method.\nWe have incorporated the above results in the revised manuscript (Lines 398-433), as well as\ndetailed discussions for the origins of false positives and false negatives to address concerns\nusers might have when approaching our dataset or using our method (Lines 573-607).\nIn response to reviewer’s suggestions, we employed siRNA knockdown (KD) of TruB1 and\nDKC1 and strived to assess the dependence of our reported sites on PUS to further prove their\nauthenticity. Unfortunately, we have encountered unforeseen technical issues which have\ndelayed our data acquisition. We believe that the significant improvements we have made would\nhave already sufficiently addressed the concerns raised by the reviewers. Therefore, we are\nsubmitting this revised manuscript without the inclusion of data from the PUS KD studies.\nHowever, we are still actively acquiring data from the PUS KD experiments and will be willing\nto provide the data if the reviewers see them remain necessary and critical.\nWe believe that this revised manuscript has been substantially improved by: 1) a more rigorous\ncharacterization of the confidence level for each candidate site; 2) a comprehensive comparison\nof our dataset with PRAISE and BID-seq to address the discrepancies; 3) improved accuracy in\nquantifying the modification levels; 4) thorough evaluation of the reproducibility of our method;\n5) site-specific validation of more sites not captured by other methods; 6) detailed discussions of\nthe origins of false negatives and false positives; as well as 7) streamlined pipeline with a single\ncommand for user-friendly analysis. We believe these revisions have greatly increased the\nusability of our dataset and the reliability of our method.\nIn summary, our revision includes additional experiments, meticulous data analysis, and detailed\ndiscussions, reinforcing the validity of the identified sites. Moreover, our method excels in\ndetecting sites overlooked by other methodologies albeit may be less effective in identifying sites\nin less favorable sequence contexts, especially those with low modification levels. We posit that\nour method represents a potent innovation and serves as a crucial complement to the existing\ntoolbox for RNA modification detection and investigation.\n2. The authors claim that 1/3 of sites have over a 50% modification rate, much higher than the\naverage mRNA modification rates estimated from other methods such as BID-Seq, and\ninconsistent with their claim that the massive number of new sites are from enrichment of lowly-\nmodified positions. The authors must more rigorously demonstrate that their estimation of\nmodification occupancy is accurate, by correlating their estimates with estimates from\nmutational-profiling methods such as BID-Seq, and by demonstrating correlation with many\nmore than 2 sites estimated by CLAP. Also, the CLAP vs ELAP-Seq comparison is difficult to\nsee in figure 4h and should have its own figure panel (with the additional comparisons).\nOur response:\nThank the reviewer for this constructive critique. In response to the reviewer's critique, we\nconducted a correlation analysis comparing the quantification results obtained from ELAP-seq\nwith those from CLAP and PRAISE (Lines 379-397). We performed CLAP and obtained\nabsolute quantification for one mRNA modification site (Figure S9c). In addition, there are 9\nother mRNA modification sites quantified by CALP previously. All these 10 mRNA sites are in\nnon-ΨU sequence context. Notably, we observed a strong correlation between ELAP-seq and\nCLAP quantifications for these sites (r = 0.83) (Figure 4m). Additionally, we correlated\nquantification by our method to that by PRAISE due to a high number of overlapping sies. While\nthe correlations for sites within non-ΨU contexts were moderate (r = 0.46-0.56), the correlation\nwas high for sites within ΨU contexts (r = 0.72) (Figure S9d). These results, combined with the\nobservation that highly modified rRNA sites are distributed at the upper end of the quantification\nspectrum, underscore the general accuracy of quantification by ELAP-seq.\nOur method relies on adjusted enrichment levels to estimate the relative modification levels,\nrather than determining exact modification ratios, a limitation of pull-down based strategies.\nNevertheless, by utilizing the absolute quantification provided by CLAP (for non-ΨU context\nsites) and PRAISE (for ΨU context sites) as benchmarks, we were able to estimate the\napproximate modification levels more confidently in this revised manuscript. This analysis\nrevealed that approximately 20-25% of sites exhibit moderate modification levels (between 30%-\n60%) or high modification levels (exceeding 60%), while over 70% of sites are characterized as\nlowly modified (below 30%) (Figure S9e). These results align more closely with conclusions\ndrawn from other studies (PMID: 36997645 and 36302989). Such semi-quantitative estimations\nof modification levels, coupled with the capacity to probe modification level changes across\nconditions, would render our method applicable to a wide range of studies.\nWith the improved quantification accuracy, we re-examined enriched GO terms for genes\ncontaining moderately and highly modified Ψ, and found a notable enrichment of categories\nrelated to mitochondrial electron transport and ATP synthesis (Figure 5e) (Lines 459-470). This\nfinding aligns with the report that PUS1 is associated with mitochondrial myopathy and calls for\nfuture investigation into the roles of Ψ in mitochondrial function and associated diseases.\n3. The language explaining the rational for cutoffs and filtering is vague. Please provide specific\ndetails with statistical justification for the following: “We optimized cutoffs for stop ratio and\nreads coverage and allow for gradients of coverage for different ranges of stop ratios” [275-276]\n“We required higher reads coverage for lower levels of stop ratios to decrease false positives”\n[685-686] “We allowed slightly lower reads coverage for sites in the out-group” [688-689]\n“…we screened cutoffs for stop ratios and coverage to find a balance between the number of\nsites identified and percentage of U” [689-691]\nOur response:\nWe appreciate the reviewers' feedback regarding the clarity of cutoffs and filtering rationale. The\ncutoffs were determined empirically. Initially, we set a stop ratio of 0.1, guided by our visual\nexamination of IGV, which indicated that some clustered sites tended to exhibit low stop ratios.\nAdditionally, we allowed a read coverage of 5, a practice used in several other studies (PMID:\n31548705 and 26075521). However, upon using the percentage of U as an indicator, we found\nthese cutoffs to be suboptimal (U% = 81%). Consequently, we further refined them with more\nstringent cutoffs. The rationale behind this optimization was requiring higher read coverage for\nsites with lower stop ratios. In addition, in-group sites situated inside of the enrichment peaks\ntypically had higher read coverage, thus necessitating overall higher coverage thresholds for\nthese sites. As the number of sites detected decreases with increasing U%, we settled on cutoffs\nwhich gave U% of 87.9%. The script used has been provided in GitHub and the cutoffs can be\nfurther user-defined for desired confidence. For example, users may raise the cutoffs for stop\nratio to filter for sites with higher confidence.\nWe have now provided the specific details in the Materials and Methods section (Lines 875-\n880):“…The optimized process was as follows: 1) we required higher reads coverage for lower\nlevels of stop ratios to decrease rate of false positives: In-group: at least 6-read coverage for sites\nwith stop ratios > 0.8, at least 8-read coverage for sites with stop ratios ranging between 0.2 and\n0.8 , and stopped reads of at least 11 for sites with stop ratios ranging between 0.1 and 0.2; Out-\ngroup: at least 5-read coverage for sites with stop ratios > 0.8, at least 7-read coverage for sites\nwith stop ratios ranging between 0.2 and 0.8 , and stopped reads of at least 2 for sites with stop\nratios ranging between 0.1 and 0.2; …”\n4. For the IFN dataset it is mentioned that “the interferon stimulated samples (IFN+) resulted in\n~40% less reads compared to untreated samples (IFN-), thus we allowed lower reads coverage\ncutoff for the IFN+ samples.”[715-719] How does this not increase the false positive rate for the\nINF+ samples – if the read cutoffs were set to minimize false positive rate? Intuitively this cutoff\nshould not be determined by read depth. This comparison would be more accurate with the INF-\nsamples randomly sub-sampled down to the same depth as INF+, using identical read cutoffs.\nOur response:\nThis is a great point. In the previous version of the manuscript, we implemented a stricter cutoff\nfor IFN- samples and a looser cutoff for IFN+ samples, based on the established thresholds used\nfor HeLa and HEK cell lines. As noted by the reviewer, this adjustment could lead to an\nincreased rate of false negatives and false positives for each sample group, respectively. In this\nrevised manuscript, we reverted to the original cutoffs, resulting in the identification of 9117 and\n3737 candidate sites from IFN- and IFN+ samples, respectively (Figure 6c and Lines 498-501).\nTo ensure more accurate comparisons, we conducted a reanalysis of the data by employing an\nequal number of sub-sampled reads for both IFN- and IFN+ samples. Utilizing the same pipeline\nwith consistent read cutoffs, we reanalyzed data subsampled from reads of libraries constructed\nwith SuperScript III, or combined libraries built with SuperScript III and SuperScript IV. In both\nscenarios, we observed a 40~50% reduction in identified sites from the IFN+ samples compared\nto the IFN- samples (Figure S12b and lines 501-503). This finding supports a consistent trend of\ndownregulation of Ψ following interferon treatment, which aligns with observations made using\nLC-MS/MS.\nMinor critiques:\n1. “a significant proportion of Ψ sites in mRNA appear in the 5’-ΨU-3’ context…” [203-204]\nwhat proportion exactly?\nOur response: Based on the paper of Li et al., ~30% of Ψ sites appear in the 5’ ΨU 3’ context.\nWe have changed our main text to provide a more precise description (Lines 209-210).\n2. “30 – 60% of Ψ in cellular mRNA can be labeled by this enzyme” [205-206] Do you expect\n30-60% of Ψ positions are being labeled – or each site is being labeled at a 30-60% rate?\nOur response:\nWe apologize for the confusion caused by this sentence. The ratio 30-60% is m1Ψ/Ψ percentage\ndetermined through LC-MS/MS analysis on digested RNAs. Due to the sequence preference of\nthe enzyme, the efficiency of enzymatic labeling should be higher for Ψ sites in the 5’-ΨU-3’\ncontext compared to sites in other sequence contexts. Ψ sites in the 5’-ΨU-3’ constitutes around\n30% of all Ψ sites in mRNA while Ψ sites in other sequence contexts constitutes about 70%.\nTherefore, more than 30% of Ψ sites would undergo labeling, with the degree of labeling at\nindividual positions contingent upon the specific sequence context, which together gave an\noverall labeling efficiency of 30-60% across all Ψ sites.\nWe have modified our main text to clarify this (Line 139): “The overall 33-60% m1Ψ/Ψ ratio\nsuggests that a considerable proportion of Ψ modification sites were successfully labeled.”\n3. Why is there a stop in untreated RNA at the expected psi position in figure 2e?\nOur response:\nWe appreciate the reviewer's observation, and upon reflection, we acknowledge that our\nmanuscript lacked clarity in describing the experiment's design.\nIn our assay, we used a synthetic 33-mer RNA as the template, which harbors a Ψ modification\nin the middle. A fluorescently labeled RT primer was designed to target the position immediately\ndownstream of the modification site. The RT primer and template were mixed at a 1:10 ratio to\ndrive most primers to hybridize with the template. In the absence of any treatment, the RT primer\nreads through the modification site, producing the full-length RT product (upper FL band).\nConversely, when encountering the propargyl-Ψ adduct (with or without subsequent\nbiotinylation at a low dNTP concentration), the elongation of the RT primer is impeded, resulting\nin a significant amount of un-extended RT primer (lower Primer band in Figure 2e).\nThe anticipated stop signatures in treated RNA align precisely with the expected Ψ position. The\nlimited presence of un-extended RT primer bands in untreated RNAs is likely attributed to a low\nlevel of unhybridized RT primer (Figure 2e, lane 2 and lane 6). Notably, the un-extended RT\nprimer bands in treated RNAs exhibited considerably stronger signals compared to those in\nuntreated RNAs (Figure 2e, lane 3-5 compared to lane 2, and lanes 7-9 compared to lane 6). This\ndisparity indicates that the introduction of bulkier groups during Ψ labeling significantly\nincreased the arrest rate of reverse transcription at the modification site, aligning with our\nexpectations.\nTo make the experimental design clearer, we have incorporated additional details in the legend\nof Figure 2 and in the Materials and Method section (Lines 714-724).\n4. How do the NNpsiNN pools perform in the standard library prep? E.g. are there any RT-stop\nor ligation biases due to sequence context?\nOur response:\nWe appreciate the reviewer's insightful question. In assessing the sequence preference of\nenzymatic labeling with the NNpsiNN pools, we employed a protocol designed to avoid\ngenerating RT stops so that we could comprehensively evaluate the sequence context spectrum\non both sides of the modification site. Specifically, following on-bead pull-down, we eluted\nRNAs from the beads by heat-denaturing the streptavidin-biotin interaction. Subsequently, we\nemployed Superscript IV with 0.5 mM final dNTP to facilitate full readthrough at the\nmodification site. Additional details about this methodology have been incorporated into the\nMaterials and Methods section for clarification (Lines 816-817).\nOur ELAP-seq results consistently identified uracil as the primary nucleobase, aligning with the\nin vitro characterization of Mj1640's specificity towards Ψ. This evidence supports the\nconclusion that ELAP-seq specifically detects Ψ sites in RNA. We want to emphasize that, in our\nlibrary preparation workflow, we do not perceive sequence context-dependent biases in RT-stop\nor ligation. While previous studies have reported sequence biases for 3’ RNA adapter ligation\nusing T4 RNA ligase 2 truncated KQ, we mitigated this bias by employing RNA adapters with\nmultiple Ns on the 5’ end, as recommended (Nucleic Acids Research, 2011, Nov; 39(21):e141).\nOur analysis of input libraries using the same pipeline revealed no enrichment for uracil,\ndemonstrating that ligation and reverse transcription were not influenced by sequence context.\n5. It is not described how the wobble pair analysis was performed, please add more detail in the\nmethods or main text.\nOur response:\nWe thank the reviewer for pointing this out. We have now provided a more detailed description\nof the analysis in the Materials and Methods section (Lines 934-946). We have also deposited the\ndata analysis pipelines and the in-house python transcripts to GitHub.\n“RNA sequences of 100 nt surrounding all Ψ sites located in different mRNA subregions (5’\nUTR, 3’ UTR and CDS) with 50 nt upstream and 50 nt downstream were used for RNA\nsecondary structure prediction. The program viennaRNA was used with the command of\n‘RNAfold -p < input.fa > output.res’. The resulting file was cleaned up in notepad++ to keep the\ninformation of the topmost prediction for each site. Then, we used an in-house python script\nmatch.py to extract the information of the base pairing status of the modification site and the\nlocation of its matched base if it is base paired. We then determined the identity of the matched\nbase for each base-paired modification site using the ‘MID’ function in excel. Next, we\ncalculated the frequency of the modification site residing in a mismatched/UG wobble status\nversus in UA matched status. Control U sites were chosen at both the -20 position and +20\nposition of Ψ sites and analysis were performed the same way as for Ψ sites. mfold were used to\ndraw the secondary structures of 100 nt surrounding the chosen Ψ sites.”\n6. Enrichment of Ψ at splice junctions is only informative with a comparison against the\nbackground distribution of U. Plot both distributions from your dataset and demonstrate that the\ndifference in Ψ distribution is significant.\nOur response:\nWe express gratitude to the reviewer for providing valuable insights into the analysis. In\nresponse to the suggestion, we randomly selected 10,000 U sites from the transcriptome and\nconducted a parallel analysis. To our disappointment, these randomly chosen sites were also\nfound to be situated within 100 nt of the nearest splice sites. To delve deeper into this\nobservation, we computed the average distance to splice sites within 800 nt, 400 nt, 200 nt, and\n100 nt windows, and still, no substantial differences were observed between the Ψ group and the\ncontrol U groups across any of these window sizes. We regret the oversight that led to an\nincorrect conclusion in the previous version of the manuscript. To rectify this, we have removed\nthe erroneous figure and texts that depicted the correlation of Ψ with splicing."}, {"round": 2, "reviewer_comments": {"reviewer_3": "During revision, the authors have incorporated a substantial number of bioinformatics analyses. These analyses encompass overlap studies with the PRAISE, BID-seq, and CeU-seq methodologies. Additionally, the authors have performed DKC1 knockdown experiments to identified PUS-dependent sites.The potential false positives and false negatives associated with the technology have been thoroughly discussed in the discussion section. However, two minor concerns remain to be addressed:\n1. The authors should provide data from NNUNN spike-in experiments to further validate the specificity of the technology.\n2. In Figure S8b, despite the authors having performed three replicates, the reproducibility of the modification sites is only approximately 20%. A detailed explanation and clarification of this discrepancy are required.\n(Remarks on code availability)", "reviewer_4": "In this manuscript, Wang et al. present a new transcriptome-wide pseudouridine (Ψ) site profiling method (ELAP-seq) using the Mj1640 enzyme to modify Ψ sites with a propargyl group and click-chemistry for biotin labeling. The authors achieved enrichment of Ψ by streptavidin pulldown and single nucleotide resolution by harnessing the RT-stop signature. For the most part, the authors have addressed concerns voiced by Reviewer 1 and Reviewer 3. The most appealing finding in this current manuscript is the detection of many more Ψ sites compared to previous methods, with some sites further validated by orthogonal methods. However, the potential strength of this enzyme-based method appears to be underdeveloped, and the inconsistency of cutoffs and filtering makes it difficult to follow. These concerns, combined with the relative simplicity, quantitativeness, and cost-effectiveness without need for commercially unavailable enzyme of other methods (particularly bisulfite-Ψ adduct methods PRAISE for transcriptome-wide Ψ profiling and BID-Seq for targeted Ψ profiling) severely limit the utility of ELAP-seq to the Ψ research community.\nAs such, the current version of this manuscript is not well suited for publication in this journal, in our opinion. The manuscript would be more useful as either a resource of novel Ψ sites (which would require a lot more validation, particularly for lower confidence sites) or a revelation of novel biology that existing methods cannot capture.\nMajor comments:\n1. The authors claim that ELAP-seq could be used in live cells for Ψ profiling, but in the ELAP-Seq protocol they denature fragmented isolated RNA prior to enzymatic labeling. This is a contradiction, because RNA secondary structure is influenced by conditions unique to cells not present in denatured fragmented isolated RNA, including long-range RNA-RNA interactions and RNA-binding proteins. For this reason, it is not obvious to the reader that in situ ELAP-Seq would work. Certainly in vivo ELAP-Seq, as hinted at in Supplementary Fig. S3, would not work.\n2. In the Introduction, the authors emphasize that many other methods for profiling Ψ sites require a high level of RNA input and/or deep sequencing. Although the authors showed the detection sensitivity of this assay on synthesized oligos in Figure 2b, the detection limit of this assay on real mRNA samples from cells is unclear (based on the Methods section, ELAP-seq appears to be carried out with 200-1000ng of RNA, but it is unclear if it is total RNA or isolated mRNA). In this aspect, it remains uncertain whether ELAP-seq can be used for rare cell populations that other methods could hardly approach.\n3. In lines 234-236 the authors decided to use data generated with Superscript III for detecting Ψ sites, while in lines 287-289 the authors decided to combine Superscript III and Superscript IV data to yield better signal. The inconsistency of data analysis makes it difficult to follow. It might be better to fix a filtering and cutoff strategy with rRNA benchmarking and then stick to the same filtering so that the level of false positives can be more easily estimated.\n4. In lines 278-280 and a few cases later in the manuscript, as well as in Fig. 4c-e, the authors detected Ψ sites at U, G, A, C with different proportions, while when describing their filter strategy (lines 246-247) they only include sites for which “the stop event occurs at uridine (U).” It is unclear why these G, A, C sites were still included with their pipeline, unless the authors really mean ΨU, ΨG, ΨA, and ΨC.\n5. In the paragraphs to study Ψ sites in response to IFN treatment, the phenotype was not well explained. By LS-MS/MS, there seemed to be a 20-25% reduction of global Ψ, while by ELAP-seq there seemed to be 41-52% fewer sites. The authors should explain this discrepancy. Additionally, since the protein level of PUS enzymes did not change, it is unclear how robust this assay was. The authors should either demonstrate Ψ change with a different method or identify other enzymes/regulators to explain the phenotype.\n6. Pull-down of both enzymatic and control reactions in ELAP-seq means that the methodology will likely be even less quantitative than other semi-quantitative methods for lower RNA inputs. It is hard to imagine that spike-ins mentioned in lines 564-565 would alleviate this concern.\nMinor comments:\n1. In the rebuttal, the authors write that, “Performed knockdown experiment for one of the pseudouridine writers, DKC1, which led to the identification of 1014 sites confidently downregulated due to DKC1 depletion, including 882 sites newly identified by ELAP-seq (Figure 4h and the new supplementary Figure S10).” The authors should analyze how much of this change is due to transcript expression changes leading to lower confidence in calling Ψ sites vs. bona fide Ψ downregulation.\n2. In line 114, the authors write, “data not shown.” They should show the data or remove/restate the sentence.\n3. In line 232-234, the authors do not show Superscript IV data.\n4. The authors should state the precise statistical test used to generate the p-value stated in line 241.\n5. In line 466, Fig. 5f should be Fig. 5e.\n6. In the Discussion in lines 604-605, the authors speculate that Mj1640 could be used to selectively introduce m1Ψ to mRNA vaccines. It is hard to imagine that such a strategy would be productive, given that such a change would inherently increase immunogenicity and reduce protein translational efficiency.\n(Remarks on code availability)"}, "author_response": "We thank the reviewers for their thoughtful comments that helped us substantially\nimprove this manuscript and strengthen the practical value of ELAP-seq. In this revision,\nwe have:\n Included additional biochemical characterization data (Figure S4) and control\nspike-in data (Figure S5), further validating the sensitivity and specificity of the\ntechnology.\n Improved reproducibility assessment with additional data analysis, addressing the\ndiscrepancies across replicates.\n Refined and benchmarked the data processing pipeline and applied it uniformly\nacross all datasets, making the analysis consistent and the interpretation of false\npositive rate more reliable.\n Provided further clarification of input requirements, quantification capabilities,\nand enzyme accessibility considerations for ELAP-seq.\n Refined the discussion of potential applications to avoid overstating in-cell\napplicability and to include additional applications suggested by new data.\n Responded to other comments raised by the reviewers, including clarifying\nspecific points of the experimental design and data analysis workflow.\nBelow we address each comment point-by-point.\nReviewer #3 (Remarks to the Author):\nDuring revision, the authors have incorporated a substantial number of bioinformatics\nanalyses. These analyses encompass overlap studies with the PRAISE, BID-seq, and\nCeU-seq methodologies. Additionally, the authors have performed DKC1 knockdown\nexperiments to identify PUS-dependent sites.The potential false positives and false\nnegatives associated with the technology have been thoroughly discussed in the\ndiscussion section. However, two minor concerns remain to be addressed:\n1. The authors should provide data from NNUNN spike-in experiments to further validate\nthe specificity of the technology.\nOur response: We have conducted ELAP-seq on the NNUNN spike-in, added a new\nfigure (Figure S5) to report the result and incorporated this result in the main text.\nImportantly, no specific sequence was enriched for the NNUNN spike-in, in clear\ncontrast to the enrichment pattern observed for NNΨNN spike-in. Moreover, compared\nwith the NNΨNN libraries, both the NNUNN and its click control (mock click chemistry)\nlibraries were built with overamplification. Altogether, these results suggest that there is\nno specific labeling or enrichment of the NNUNN spike-in oligos.\nThe high specificity of the enzymatic labeling assay was also demonstrated by the new\ndata added in Figure S4 comparing labeling on modified vs unmodified RNAs (in the\nNNΨ/UNN context) in dot blot assays.\n2. In Figure S8b, despite the authors having performed three replicates, the\nreproducibility of the modification sites is only approximately 20%. A detailed\nexplanation and clarification of this discrepancy are required.\nOur response: Thanks to the reviewers for this insightful question. We noticed that some\nsites with low read coverage may lack any coverage in other replicates simply due to\ntechnical variation, leading to low overlap. To remove the effect of technical variation at\nlow-coverage sites and better evaluate reproducibility in identifying modification-specific\nsignals, we focused on sites covered by at least five reads in at least one other replicate\n(within our filtering criterion of ≥5 stopped reads) and then applied our standard filter\ncriteria for calling RT stop sites. This updated analysis yielded markedly improved\nreproducibility: approximately 64% of stop sites in each replicate were detected in at\nleast one other replicate, and 32–39% of sites were shared across all three replicates\n(Figure S10c).\nTo further address the remaining discrepancies, we examined RT stop sites identified in\nonly one replicate. Because RT stop sites are defined mainly based on the number of\nstopped reads and the stop ratio, we compared the number of stopped reads between sites\nunique to a single replicate and those detected in multiple replicates. For example, sites\nidentified only in replicate 1 (Rep1-only sites) exhibited substantially fewer stopped\nreads than sites also detected in other replicates, and most fell below the detection\nthreshold in the other replicates (Figure S10d). This indicates that the discrepancies\nmainly arise from sites with near-threshold stopped reads counts. In this study, we define\ncandidate modification sites as those detected in at least two independent replicates,\nwhich effectively enriches for sites with robust signals.\nThese updated results are now shown in Figure S10b-d. We have added this explanation\nto the main text (page 16, lines 317-322). The associated code was also updated in Github\n(code availability section).\nReviewer #4 (Remarks to the Author):\nIn this manuscript, Wang et al. present a new transcriptome-wide pseudouridine (Ψ) site\nprofiling method (ELAP-seq) using the Mj1640 enzyme to modify Ψ sites with a\npropargyl group and click-chemistry for biotin labeling. The authors achieved\nenrichment of Ψ by streptavidin pulldown and single nucleotide resolution by harnessing\nthe RT-stop signature. For the most part, the authors have addressed concerns voiced by\nReviewer 1 and Reviewer 3. The most appealing finding in this current manuscript is the\ndetection of many more Ψ sites compared to previous methods, with some sites further\nvalidated by orthogonal methods. However, the potential strength of this enzyme-based\nmethod appears to be underdeveloped, and the inconsistency of cutoffs and filtering\nmakes it difficult to follow. These concerns, combined with the relative simplicity,\nquantitativeness, and cost-effectiveness without need for commercially unavailable\nenzyme of other methods (particularly bisulfite-Ψ adduct methods PRAISE for\ntranscriptome-wide Ψ profiling and BID-Seq for targeted Ψ profiling) severely limit the\nutility of ELAP-seq to the Ψ research community. As such, the current version of this\nmanuscript is not well suited for publication in this journal, in our opinion. The\nmanuscript would be more useful as either a resource of novel Ψ sites (which would\nrequire a lot more validation, particularly for lower confidence sites) or a revelation of\nnovel biology that existing methods cannot capture.\nWe thank the reviewer for the constructive comments. We would like to address the\nreviewers’ concerns as below:\nUtility of ELAP-seq: We have successfully applied the method to multiple cell lines and\nconditions, reproducing many previously reported Ψ sites and uncovering new ones,\nmany of which have been directly or indirectly validated. These newly detected sites will\nprovide a valuable resource and can guide future targeted investigations using low-\nthroughput assays (e.g., CMC-based primer extension assays). Notably, 10% of ELAP-\nseq identified sites are of low confidence, thus the majority of sites are reliable.\nWhile enrichment-based approaches are semi-quantitative, they have been highly\ninformative for RNA modification studies, such as the MeRIP-seq method for profiling\nN6-methyladenosine. Based on the existing data and new data presented in Figure S4 and\nS5, ELAP-seq shows improved specificity and sensitivity and will outperform previous\nchemical- or antibody-based enrichment strategies.\nMj1640 is highly stable and can be readily prepared in most molecular biology or\nbiochemistry laboratories in a few days without requiring repeated purification.\nAlternatively, institutional core facilities and commercial providers (e.g., GenScript,\nThermo Fisher Scientific, OriGene, Novus Biologicals) can supply purified protein at\nmodest cost. To support accessibility, we have deposited the Mj1640 expression plasmid\nto Addgene.\nMany widely used biochemical modification-mapping methods similarly make use of\npurified enzymes (e.g., SAC-m6A-seq (PMID:35288668) and eTAM-seq for m6A (PMID\n36593412)).\nAdvantages of ELAP-seq: 1) a high signal-to-noise ratio; 2) low sequencing depth and\ncomputational demand compared to non-enrichment-based methods, making it cost- and\ntime-efficient; and 3) requiring a lower RNA input compared to earlier enrichment-based\nΨ mapping strategies.\nBroad value of this study to the Ψ community: In addition to developing ELAP-seq, this\nstudy establishes an efficient and specific enzymatic labeling assay for pseudouridine.\nThe labeling reaction, combined with click chemistry, offers a versatile platform for\npotential future applications such as: 1) sensitive detection of Ψ for diagnosis purposes\n(based on new data in Figure S4 and page 7, lines 142-157); 2) enrichment of modified\nRNAs for interactome studies; and 3) compatibility with future low-input or in-cell\nadaptations.\nMajor comments:\n1. The authors claim that ELAP-seq could be used in live cells for Ψ profiling, but in the\nELAP-Seq protocol they denature fragmented isolated RNA prior to enzymatic labeling.\nThis is a contradiction, because RNA secondary structure is influenced by conditions\nunique to cells not present in denatured fragmented isolated RNA, including long-range\nRNA-RNA interactions and RNA-binding proteins. For this reason, it is not obvious to the\nreader that in situ ELAP-Seq would work. Certainly in vivo ELAP-Seq, as hinted at in\nSupplementary Fig. S3, would not work.\nOur response: We thank the reviewer for this insightful comment. We agree that the\ncurrent ELAP-seq protocol involves denaturation of fragmented, isolated RNA prior to\nenzymatic labeling, and that such conditions differ from those in live cells. In our study,\nwe observed that when the enzyme is expressed in mammalian cells, it can label\nendogenous modifications with approximately 10% efficiency. The apparently lower\nlabeling efficiency in cells compared to that in vitro (Figure 1) may be due to factors\naffecting accessibility of the enzyme to RNA such as protein protection, RNA tertiary\nstructure, and local enzyme versus substrate concentrations.\nWe agree that the current protocol is not yet optimized for in situ or in vivo mapping,\nand\nwe have now toned this down in the main text to avoid overclaiming. Future\noptimizations may enable these applications, such as fusing it with helicase or strand-\ndisplacement domains, to transiently unwind structured RNAs for better labeling.\nDespite this, our data show that the enzyme is catalytically active in cells and thus has\npotential for in-cell labeling of pseudouridine. We envision that this system could be\ndeveloped further as an in-cell site-directed RNA modification tool, where a guide RNA\nwill be used to disrupt local RNA structure, similar to other RNA-guided modification or\nediting systems. Since pseudouridine modification is irreversible, no current method\nexists to specifically remove it at defined sites to study its biological function. We\npropose that our approach could be used to perturb this modification and investigate its\nfunctional role. We highlight this as a potential future application in the discussion\nsection (page 29, lines 586-591).\n2. In the Introduction, the authors emphasize that many other methods for profiling Ψ\nsites require a high level of RNA input and/or deep sequencing. Although the authors\nshowed the detection sensitivity of this assay on synthesized oligos in Figure 2b, the\ndetection limit of this assay on real mRNA samples from cells is unclear (based on the\nMethods section, ELAP-seq appears to be carried out with 200-1000ng of RNA, but it is\nunclear if it is total RNA or isolated mRNA). In this aspect, it remains uncertain whether\nELAP-seq can be used for rare cell populations that other methods could hardly\napproach.\nOur response: In the Introduction, we referred specifically to CeU-seq, which functions\nsimilarly to ELAP-seq as an enrichment-based Ψ mapping method, that requires much\nhigher RNA input (10 µg of mRNA, approximately 50-fold higher) than our method.\nSimilarly, another enrichment-based method using antibody also requires 10 µg of\nmRNA (Campos et al., PMID: 34376564). For other existing methods, our discussion\nfocused primarily on their need for deep sequencing, rather than their RNA input\nrequirements. We have further clarified this in the Discussion, under the section of\nAdvantages and Limitations of ELAP-seq (page 27, line 532).\nAs stated in the materials and method section, 200–1000 ng of RNA (mRNA) was used\nfor the initial fragmentation step, with a typical yield of ~80-100%. From this fragmented\nRNA, we aliquoted ~200 ng for downstream library preparation. Therefore, ELAP-seq\ncurrently requires ~250 ng mRNA input and is comparable to PRAISE in input\nrequirements (500 ng mRNA). These details were included in the original version and are\nnow further clarified in the revised Materials and Methods section (page 34, lines 692-\n693).\nEnrichment-based methods generally require higher RNA input than non-enrichment-\nbased approaches, due to loss of materials during binding and washing. ELAP-seq uses\nmild labeling conditions, which preserves RNA, and the streptavidin-biotin pair, which\nprovides high binding affinity and can withstand stringent washes. These features make\nthe required RNA input lower than previous enrichment-based methods and comparable\nto that of recent non-enrichment-based methods.\nWe clarify that the current implementation is not intended for rare-cell applications. It is\nnoteworthy that the use of a low dNTP concentration in our library preparation pipeline\nnecessitates a starting material of 250 ng of mRNA. We are actively developing a next-\ngeneration, low-input version that does not require low dNTP or enrichment, enabling\nsubstantially lower input requirements and reduced background.\n3. In lines 234-236 the authors decided to use data generated with Superscript III for\ndetecting Ψ sites, while in lines 287-289 the authors decided to combine Superscript III\nand Superscript IV data to yield better signal. The inconsistency of data analysis makes it\ndifficult to follow. It might be better to fix a filtering and cutoff strategy with rRNA\nbenchmarking and then stick to the same filtering so that the level of false positives can\nbe more easily estimated.\nOur response: We appreciate the reviewer’s comments about the inconsistency in data\nanalysis. We have now explained the rationale for using Superscript III and IV data. We\nalso applied refined and fixed filtering criteria throughout the analysis. These criteria\nwere first benchmarked using rRNA to ensure robustness before being applied to mRNA\ndata. This new analysis changed the numbers of identified candidate modification sites to\n5,366 and 5,741 for HeLa and HEK293T poly(A) RNAs, respectively. We have updated\nthe manuscript throughout accordingly.\nWe have modified the main text to clarify the usage of Superscript III and IV (page 15-\n16, lines 291-306):\n“While Superscript III with low dNTP produces prominent RT stop signature, we noticed\nthat the coverage at some sites in mRNA is below the reads coverage cutoff although\nboth replicates share strong stop pattern. To achieve higher read coverage, we also used\nsuperscript IV with the commercially recommended concentration of dNTP (0.5 mM\nfinal) to generate libraries in parallel from the same pull-down product. Interestingly,\nsuperscript IV also generated stop signatures at Ψ, although to a lesser extent (Figure\nS7a), suggesting that the on-bead RT reaction used in our library construction protocol\nfurther facilitated the generation of the stop signature (Figure S7a). To rescue additional\nsites by increasing the reads coverage, we combined data from superscript IV and III for\nfurther analysis. We first tested this analysis strategy on rRNA and verified that, with\nslightly more stringent cutoff in stop ratio (See Materials and Methods section), this\nstrategy did not introduce additional false positives. Additionally, this analysis enabled\nthe identification of the two modification sites in 5.8S rRNA, which were previously\nmissed due to low coverage (Figure S7e).\nApplying the same analysis on mRNA identified an additional 1,628 and 2,410 candidate\nsites from HeLa and HEK293T poly(A) RNA, respectively, yielding a total of 5,366\ncandidate sites for HeLa and 5,741 for HEK (Figure 4d, Figure S9b). Of these, 2, 531\ncandidate Ψ sites are shared between HEK293T cells and HeLa cells, accounting for\n~46% of all candidate sites in each cell line (Figure 4f). We note that some sites passing\nthe filter in one cell line but found in only one replicate in the other cell line were not\nincluded. Therefore, the actual number of shared sites between the two cell lines may be\nhigher than reported here.”\nWe have updated and moved the description of the filter criteria to the Materials and\nmethods section:\nOn page 39-40, lines 813-823: “A candidate site in rRNA is defined when meeting the\nfollowing criteria: 1) the stop ratio at the site in the pull-down library is > 0.05 in both\nreplicates; 2) stop ratio in the pull-down library – stop ratio in the input library is > 0.05\nin both replicates; 3) the number of stopped reads in the pull-down library is >=5 in both\nreplicates; 4) the average value of stop ratio * stopped reads between the two pull-down\nreplicates is >=1.5 (this imposes a more stringent requirement on the number of stopped\nreads when the stop ratio in the pull-down sample is lower than 30%); and 5) the stop\nevent occurs at uridine (U). Finally, we further remove sites meeting the following\ncriteria: 1) the stop ratio in the input library is > 0.1 and > 1/3 of the stop ratio in the pull-\ndown library, and 2) number of stopped reads in input is > =3 in any replicate. When\nusing superscript III and IV data together, we follow the same criteria except for\nrequiring that the stop ratio in the pull-down library- stop ratio in the input library is >\n0.1.”\nWe have added more details about the library preparation:\nOn page 36, lines 730-733: “Beads after pull-down were divided to two parts for the\ndownstream reverse transcription step, with three-quarters (corresponding to 150 ng\nstarting RNA) used for SuperScript III RT and one-quarter (corresponding to 50 ng\nstarting RNA) for SuperScript IV RT.”\n4. In lines 278-280 and a few cases later in the manuscript, as well as in Fig. 4c-e, the\nauthors detected Ψ sites at U, G, A, C with different proportions, while when describing\ntheir filter strategy (lines 246-247) they only include sites for which “the stop event\noccurs at uridine (U).” It is unclear why these G, A, C sites were still included with their\npipeline, unless the authors really mean ΨU, ΨG, ΨA, and ΨC.\nOur response: We appreciate the opportunity to clarify this point. The figures in Fig. 4c–e\nshow the distribution of reverse transcription (RT) arrest sites at A, U, C, and G, rather\nthan the final calls of pseudouridylation. Only U positions exhibiting RT-stop signatures\nare subsequently defined as candidate Ψ sites. Our analysis pipeline records RT stops at\nall four nucleotides (A, U, C, G), but eventually filters for stop events occurring at U to\nidentify potential Ψ modification sites. To be clear, we do not infer pseudouridylation at\nG, A, or C; these positions are presented solely to illustrate background RT-arrest\nfrequencies and to assess specificity, as elaborated below.\nFor mRNA, where modification sites are largely unknown, evaluating stops across non-U\nnucleotides provides a measure of specificity and potential false positives. This is based\non the rationale that, if no specific labeling and enrichment of pseudouridine occurs, RT\nstop would occur randomly with similar frequencies at A, U, C and G. Instead, we\nobserved strong U-specific enrichment (>95% for SuperScript III and >90% for\nSuperScript IV) of RT stop signatures in the pull-down samples whereas input samples\nshowed no such enrichment, confirming that detected stops are specific to pseudouridine,\nrather than resulting from RNA structure or sequencing artifacts. Similar specificity\nanalyses have been used in prior studies (Zhang et al.,PMID 26075521 and Cui et al.,\nPMID 33313824).\nIn contrast, for rRNA, we omitted reporting results at non-U because the well-\ncharacterized Ψ modification sites on rRNA can be directly used for specificity\nvalidation, and we don’t need to use the data at non-U sites for this analysis. For rRNA,\nonly two non-U stops were detected (Supplementary table 5), consistent with the very\nlow false positive rate and supporting the reliability of using non-U stops to assess\nspecificity in mRNA.\nIn summary, non-U positions are included in the pipeline for specificity assessment in\nmRNA. Only the U sites are further defined as Ψ modification sites. We have now further\nclarified this in the main text (page 15, lines 277-290).\n5. In the paragraphs to study Ψ sites in response to IFN treatment, the phenotype was not\nwell explained. By LS-MS/MS, there seemed to be a 20-25% reduction of global Ψ, while\nby ELAP-seq there seemed to be 41-52% fewer sites. The authors should explain this\ndiscrepancy. Additionally, since the protein level of PUS enzymes did not change, it is\nunclear how robust this assay was. The authors should either demonstrate Ψ change with\na different method or identify other enzymes/regulators to explain the phenotype.\nOur response: We appreciate the reviewer’s insightful comments regarding the\ndiscrepancy between LC-MS/MS and ELAP-seq results. The discrepancy likely arises\nfrom the different nature of the two assays: LC-MS/MS quantifies the total Ψ content at\nthe nucleotide level, whereas ELAP-seq detects site-specific Ψ modifications, which\ndepend on detection thresholds. As shown in the figure below, the sites that become\nundetectable after IFN treatment exhibit substantially lower read coverage than sites that\nremain detectable. Thus, a modest reduction in modification level could cause many low-\ncoverage sites to fall below the ELAP-seq detection threshold, resulting in a larger\napparent decrease in the number of detectable sites (~41–52%) than in total Ψ abundance\n(~20–25%).\nWe agree with the reviewer that the study of the IFN treatment is incomplete due to the\nunresolved mechanism underlying Ψ downregulation phenotype. Fully resolving the\nmechanism would require further investigation beyond PUS protein abundance,\npotentially involving post-translational regulation, protein subcellular localization,\nsubstrate availability, or unknown RNA-binding regulators. As this mechanistic\nclarification is beyond the scope of the present study, we have removed the IFN treatment\ndata from this manuscript (previously Figure 6 and S17) and updated the text\naccordingly. We are actively continuing this line of research as part of a follow-up study.\nImportantly, this removal does not affect the robustness of ELAP-seq for differential\nmodification analysis, which has been demonstrated in the DKC1 knockdown\nexperiment. In the DKC1 knockdown experiment, the DKC1-dependent sites determined\nby PRAISE are also predominantly downregulated in our dataset, supporting the reliable\ndetection of site-specific differential modification using our method.\nIn addition, this removal does not affect the core findings of the paper, which center on\n(1) the biochemical characterization of Mj1640, (2) the development of ELAP-seq, and\n(3) the validation of Ψ detection across transcriptomes.\nIn the revised version, we have strengthened the manuscript by adding new comparisons\nwith antibody-based detection methods (Figure S4, page 7 lines 142-157) and control\nspike-in validation data (Figure S5), refining the ELAP-seq analysis pipeline, and\nexpanding discussion of future applications.\n6. Pull-down of both enzymatic and control reactions in ELAP-seq means that the\nmethodology will likely be even less quantitative than other semi-quantitative methods for\nlower RNA inputs. It is hard to imagine that spike-ins mentioned in lines 564-565 would\nalleviate this concern.\nOur response: We thank the reviewer for this comment.\nPull-down of the control group is used solely as quality control during library\npreparation. Under the same conditions, only modified RNAs in the experimental group\nare expected to be successfully pulled down and yield detectable libraries, whereas the\ncontrol group is expected to produce no DNA in the end (Figure S9a), unless being\noveramplified relative to the experimental group. Therefore, the control reaction does not\nproceed to sequencing and is not required for data analysis.\nIn our analysis pipeline, quantification starts with calculating the enrichment level at each\nsite, defined by normalizing the RPM at that site in the pull-down sample by that in the\ninput sample to account for gene expression difference. This enrichment value is further\nnormalized to account for sequence context difference, generating the normalized\nenrichment. The relative values of normalized enrichment reflect the relative\nmodification levels across different sites. Importantly, our measurements showed\nmoderate to high correlations with other quantitative methods at overlapping sites (Figure\n4j and S15e). Based on these correlations, we determined the correlation between the\nvalues of normalized enrichment and absolute modification levels, and then accordingly\nassigned remaining non-overlapping sites to various modification level categories (Figure\nS15f and the Materials and Methods section).\nWe envision that including barcoded spike-ins in various sequence contexts with defined\nmodification stoichiometries can improve the current quantification strategy. These\nspike-ins will undergo the same enzymatic labeling and pull-down treatments as\nbiological RNAs, providing calibration curves that allow direct estimation of\nmodification levels in various sequence contexts. Such spike-in calibration methods have\nbeen adopted by many studies (e.g., BACS (PMID:39349603) and PRAISE\n(PMID:36997645)).\nWe also note that reads coverage can influence quantification accuracy: mRNA sites in\nthe ΨU context, rRNA sites and mRNA sites in the non-ΨU context readily detectable by\nCLAP, which generally have high read coverage, correlate better with PRAISE\nquantification than mRNA sites in non-ΨU context, which generally have low read\ncoverage (Figure 4i and S15e). Therefore, optimizing enzyme sequence preference will\nalso help to address concerns about quantification accuracy. We have now updated the\ndiscussion section to further clarify this (page 27, lines 536-544).\nMinor comments:\n1. In the rebuttal, the authors write that, “Performed knockdown experiment for one of\nthe pseudouridine writers, DKC1, which led to the identification of 1014 sites confidently\ndownregulated due to DKC1 depletion, including 882 sites newly identified by ELAP-seq\n(Figure 4h and the new supplementary Figure S10).” The authors should analyze how\nmuch of this change is due to transcript expression changes leading to lower confidence\nin calling Ψ sites vs. bona fide Ψ downregulation.\nOur response: Our method to measure differential RNA modification uses changes in the\nenrichment level, defined as the RPM in the pull-down library normalized by that in the\ninput library (please see main text at page 20, lines 380-381 and Material and Method\nsection at page 42, lines 874-876), which already takes into account the change in the\ntranscriptional expression level. Therefore, these are bona fide downregulation, instead of\nfalse positives due to a decrease in the gene expression.\n2. In line 114, the authors write, “data not shown.” They should show the data or\nremove/restate the sentence.\nOur response: Thank the reviewer for pointing this out. Since these negative data were\nnot systematically documented for inclusion in the manuscript, we have now removed\nthis sentence.\n3. In line 232-234, the authors do not show Superscript IV data.\nOur response: This data has now been shown in figure S7a.\n4. The authors should state the precise statistical test used to generate the p-value stated\nin line 241.\nOur response: Please see the updated filtering and cutoff strategy in the Materials and\nMethod section (page 39-40, lines 813-823).\n5. In line 466, Fig. 5f should be Fig. 5e.\nOur response: This has been fixed.\n6. In the Discussion in lines 604-605, the authors speculate that Mj1640 could be used to\nselectively introduce m1Ψ to mRNA vaccines. It is hard to imagine that such a strategy\nwould be productive, given that such a change would inherently increase immunogenicity\nand reduce protein translational efficiency.\nOur response: We originally proposed such an application based on the hypothesis that\nthe modification may work through decorating a structurally important region in RNA\n(Nance and Meier, PMID 34075344). However, a recent publication revealed that the\nmodification inhibits TLR7 and TLR8 activation through blocking degradation by\nRNAse T2 (Bérouti et al., PMID 40580950), making our original hypothesis less likely.\nTherefore, we have deleted this from the discussion. We have added other potential\napplications for our method in the discussion section (page 29-30, lines 583-598).\nReviewer #3 (Remarks to the Author):\nRegarding the first point, while the inclusion of the spike-in results in Figure S5 is appreciated,\nthe color scheme requires adjustment as it is currently visually harsh. Additionally, the figure\nlegend needs to be expanded to provide more detailed descriptive information.\nOur response: We thank the reviewer for the helpful suggestion. Figure S5 has now been\nupdated with a different color scale, and the legend has been expanded with more details.\nReviewer #3 (Remarks on code availability):\nThe provided code is user-friendly and the documentation is sufficient. I was able to replicate\nthe installation and execution processes seamlessly.\nWe thank the reviewer for the positive feedback.\nReviewer #4 (Remarks to the Author):\nIn their updated manuscript, Wang et al. substantially improve their report on the\nenzymebased pseudouridine detection method ELAP-Seq. Beyond demonstrating enhanced\nsensitivity of the method, they have removed incomplete data and unsupported claims. With a\nfew minor changes, their work could be suitablel for publication in Nature Communications.\nWe thank the reviewer for the positive comments of our revised manuscript.\n1. The authors should add accessibility information on the Addgene-deposited Mj1640\nexpression plasmid to the Materials and Methods.\nOur response: We have added the deposit ID of the plasmid in the Materials and Methods\npaper (page 30). The plasmid will be released after the publication of this manuscript.\n2. The authors should remove references to the interferon treatment experiment from the\nMaterials and Methods.\nOur response: The section about interferon treatment has been removed from the Materials\nand Methods section.\n3. I am not convinced by the authors’ rebuttal to my previous minor comment 1.\na. For reference, I originally wrote: “In the rebuttal, the authors write that, “Performed\nknockdown experiment for one of the pseudouridine writers, DKC1, which led to the\nidentification of 1014 sites confidently downregulated due to DKC1 depletion, including 882\nsites newly identified by ELAP-seq (Figure 4h and the new supplementary Figure S10).” The\nauthors should analyze how much of this change is due to transcript expression changes\nleading to lower confidence in calling Ψ sites vs. bona fide Ψ downregulation.”\nb. The authors respond: “Our method to measure differential RNA modification uses\nchanges in the enrichment level, defined as the RPM in the pull-down library normalized by\nthat in the input library (please see main text at page 20, lines 380-381 and Material and\nMethod section at page 42, lines 874-876), which already takes into account the change in the\ntranscriptional expression level. Therefore, these are bona fide downregulation, instead of\nfalse positives due to a decrease in the gene expression.”\nc. The specific problem I have with this response is what the authors write in the previous\nparagraph of their rebuttal: “We also note that reads coverage can influence quantification\naccuracy: mRNA sites in the ΨU context, rRNA sites and mRNA sites in the non-ΨU context\nreadily detectable by CLAP, which generally have high read coverage, correlate better with\nPRAISE quantification than mRNA sites in non-ΨU context, which generally have low read\ncoverage (Figure 4i and S15e).”\nd. There appears to be a clear contradiction—if lower read coverage results in lower\nquantification accuracy, presumably it also results in lower confidence of called Ψ sites. If this\nis true, then decreased gene expression can indeed lead to an increase in sites confidently\ndownregulated. As is common in RNA enrichment method benchmarking, the authors should\nperform additional analysis of their RNA-seq data (Fig. S12a). Specifically, they should assess\nthe degree to which genes harboring 717 (previously 1014—I am not sure what caused the\nchange) Ψ sites confidently downregulated are themselves downregulated in gene expression.\nI propose that the authors first generate a differential gene expression analysis volcano plot of\nall RNA-Seq genes with the genes harboring these 717 Ψ sites plotted in a different color.\nAdditionally, they should quantify the distribution in fold-change of the Ψ site-harboring\ngenes relative to the fold-change of all genes. Such analysis will suitably resolve this open\nquestion.\nOur response: We thank the reviewer for the insightful comments and apologize for not\npreviously fully addressing the concern. We now understand that the key issue is whether\nreduced gene expression upon DKC1 knockdown could lead to insufficient starting material\nfor efficient pull-down and thereby cause apparent downregulation rather than true\nmodification level change. Our analyses below indicate that the degree of gene expression\nchange is insufficient to affect the confident detection of downregulated Ψ site, which are\nnow incorporated into Figure S12d–g.\nWhen evaluating modification level changes, we used enrichment level at the modification\nsite, which is defined as the local read coverage per million mapped reads (RPM) in the pull-\ndown library normalized by that in the input library, which accounts for changes in\nexpression levels (a strategy similarly applied in analyzing differential modification in\nMeRIP-seq data). To evaluate how much of this observed change is due to transcript\nexpression changes leading to lower confidence in calling Ψ sites versus bona fide Ψ\ndownregulation, we performed the analysis described below:\n1. We analyzed differential gene expressions for genes harboring the downregulated\nmodification sites and for all genes detected by RNA-seq using data from the input\nlibraries (prior to pull-down), using DESeq2. We first looked at overall trends of gene\nexpression changes without applying cutoffs for p-value or fold change, and the\nresults are presented in figure S12d-f. The results showed that there is an equal\ndistribution of upregulation and downregulation for all genes (Figure S12d-e),\nwhereas there are more genes being downregulated than being upregulated for\ngenes harboring the downregulated modification sites (Figure S12d-e). The median\nvalues of gene expression fold-changes are 1.03 for all genes and 0.83 for genes\nharboring the downregulated modification sites (Figure S12f). When applying a p-\nvalue cutoff of 0.05, only 262 genes harboring the downregulated modifications are\nconfidently downregulated, and when additionally applying a fold-change cutoff of\n2, only 2 genes cross the thresholds (Figure S12d). These overall suggest that\nalthough the expression of genes harboring the downregulated modification sites\ntend to be overall downregulated upon DKC1 knockdown, the extent of that\ndownregulation is mild. Given that modification level changes were assessed by\nnormalizing pull-down-RPM values by input-RPM values to account for gene\nexpression changes, and that stringent cutoffs (p-value < 0.05 and fold change > 2)\nwere applied to define downregulated sites, the mild gene expression changes are\nunlikely to affect the accuracy of the differential modification analysis.\n2. To further evaluate whether the reduced gene expression after DKC1 knockdown\ncould prevent efficient pull-down in the absence of true modification changes, we\nperformed an additional analysis. Genes harboring Ψ sites detected before\nknockdown represent transcripts with sufficient expression levels to support reliable\npull-down and modification detection. Under the assumption that Ψ modification\nlevels remain unchanged, transcripts whose post-knockdown expression levels\nremain within this range should still be efficiently pulled down. To test this, we\ncompared the distribution of normalized gene-level read counts for genes harboring\ndownregulated Ψ sites after DKC1 knockdown with the distribution for all genes\nharboring Ψ sites prior to knockdown. We found that these distributions are highly\ncomparable (Fig S12g), suggesting that despite modest expression decreases, the\npost-knockdown expression levels remain sufficient for efficient pull-down (Fig 12g).\nBased on these results, we updated the main text on page 17: “The changes in\nexpression of genes harboring these downregulated sites are mild and unlikely to\nsubstantially affect the efficiency of pull down of this modification or the reliability of\nthe detected Ψ downregulation events (Figure S12d-g).”\n3. To further increase the confidence in calling downregulated Ψ sites, we applied\nadditional conservative filtering in this revision. Specifically, from the list of\ndownregulated sites, we removed sites with average pull-down RPM values < 1.2\nbefore treatment (corresponding to at least 10× local read coverage at the\nmodification given the current sequencing depth of (6-10M mapped reads)), as such\nlow-coverage sites are more susceptible to technical variation or subtle decrease in\ngene expression. This filtering reduced the number of confidently downregulated Ψ\nsites to 654. The updated results are shown in Fig. 4h.\nAdditionally, to focus on genes with considerable expression levels and to enable\nmore reliable quantification of enrichment fold change, readers can apply a cutoff of\n1 for Input-RPM under both siCtrl and siDKC1 conditions (which corresponds to at\nleast 10x local read coverage at the modification under the current sequencing depth\n(10-20 M mapped reads) (PMID: 31870409)). Notably, 81% of the downregulated\nsites are above this cutoff, suggesting high reliability for most sites. This cutoff can\nbe user-defined based on their sequencing depth and research needs.\nThese changes have been incorporated into the Materials and Methods section: “For\ndownregulated sites, we further excluded sites with average pull-down RPM values <\n1.2 before treatment (corresponding to approximately 10× local read coverage given\nthe current total mapped reads of 6-10M), as their enrichment measurements are more\nsusceptible to technical variation. For more reliable fold-change quantification,\nreaders can apply a cutoff of 1 for Input-RPM under both siCtrl and siDKC1 conditions,\nwhich corresponds to at least 10x local read coverage given the current total mapped reads of\n10-20 M. Such cutoff can be user-defined based on their sequencing depth and\nresearch needs” (page 44)\n4. For clarity, the reduction in the number of downregulated Ψ sites across manuscript\nrevisions reflects the cumulative application of increasingly stringent analysis criteria,\nrather than changes in the underlying biological conclusions. During the last revision,\nwe applied an improved and more stringent analysis pipeline which has decreased\nthe total number of confident modification sites detected before treatment and\nconsequently reduced the number of sites identified as downregulated upon DKC1\nknockdown. The additional conservative filtering applied in the current revision\nfurther refined this list."}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have addressed the two minor comments raised during the revision. Regarding the first point, while the inclusion of the spike-in results in Figure S5 is appreciated, the color scheme requires adjustment as it is currently visually harsh. Additionally, the figure legend needs to be expanded to provide more detailed descriptive information.\n(Remarks on code availability)", "reviewer_4": "In their updated manuscript, Wang et al. substantially improve their report on the enzyme-based pseudouridine detection method ELAP-Seq. Beyond demonstrating enhanced sensitivity of the method, they have removed incomplete data and unsupported claims. With a few minor changes, their work could be suitablel for publication in Nature Communications.\n1. The authors should add accessibility information on the Addgene-deposited Mj1640 expression plasmid to the Materials and Methods.\n2. The authors should remove references to the interferon treatment experiment from the Materials and Methods.\n3. I am not convinced by the authors’ rebuttal to my previous minor comment 1.\na. For reference, I originally wrote: “In the rebuttal, the authors write that, “Performed knockdown experiment for one of the pseudouridine writers, DKC1, which led to the identification of 1014 sites confidently downregulated due to DKC1 depletion, including 882 sites newly identified by ELAP-seq (Figure 4h and the new supplementary Figure S10).” The authors should analyze how much of this change is due to transcript expression changes leading to lower confidence in calling Ψ sites vs. bona fide Ψ downregulation.”\nb. The authors respond: “Our method to measure differential RNA modification uses changes in the enrichment level, defined as the RPM in the pull-down library normalized by that in the input library (please see main text at page 20, lines 380-381 and Material and Method section at page 42, lines 874-876), which already takes into account the change in the transcriptional expression level. Therefore, these are bona fide downregulation, instead of false positives due to a decrease in the gene expression.”\nc. The specific problem I have with this response is what the authors write in the previous paragraph of their rebuttal: “We also note that reads coverage can influence quantification accuracy: mRNA sites in the ΨU context, rRNA sites and mRNA sites in the non-ΨU context readily detectable by CLAP, which generally have high read coverage, correlate better with PRAISE quantification than mRNA sites in non-ΨU context, which generally have low read coverage (Figure 4i and S15e).”\nd. There appears to be a clear contradiction—if lower read coverage results in lower quantification accuracy, presumably it also results in lower confidence of called Ψ sites. If this is true, then decreased gene expression can indeed lead to an increase in sites confidently downregulated. As is common in RNA enrichment method benchmarking, the authors should perform additional analysis of their RNA-seq data (Fig. S12a). Specifically, they should assess the degree to which genes harboring 717 (previously 1014—I am not sure what caused the change) Ψ sites confidently downregulated are themselves downregulated in gene expression. I propose that the authors first generate a differential gene expression analysis volcano plot of all RNA-Seq genes with the genes harboring these 717 Ψ sites plotted in a different color. Additionally, they should quantify the distribution in fold-change of the Ψ site-harboring genes relative to the fold-change of all genes. Such analysis will suitably resolve this open question.\n(Remarks on code availability)"}, "author_response": ""}]

Decision letter