A union list of m⁵C sites in human transcriptomes

BsRNA-Seq data from seven human tissues (two replicates each) and five human cell lines (typically replicated; several conditions of altered NSUN gene expression) were collected from six published studies (Yang et al, 2017; Chen et al, 2019; Huang et al, 2019; Janin et al, 2019; Schumann et al, 2020; Liu et al, 2021a) (Table S1). These were re-analysed with a pipeline for m⁵C site calling based on (Schumann et al, 2020) with minor modifications (see the Materials and Methods section). Briefly, after data pre-processing and mapping to the C-to-T and G-to-A converted genome (plus spike-in sequences as appropriate), non-conversion sites were called as variations by a custom script, for example, T-to-C in read2 and A-to-G in read1 (Fig 1A). Read distribution and overall cytosine conversion rate were checked to assess quality of the 50 datasets. For seven of the primary tissue samples, one of the two replicate datasets was excluded as metagene analysis showed strongly biased coverage toward mRNA 3′ ends, suggesting excessive degradation (Fig S1). Consistent with the original reports (e.g., [Schumann et al, 2020]), the remaining 43 datasets each showed good library complexity (Fig S2A) and high cytosine conversion rates, for example, ≥99.6% conversion rate for all annotated genes and ≥99.7% for protein coding genes (Table S1, Fig S2B, top panel). A set of filters was then applied to identify high-confidence m⁵C candidate sites. Reads with more than three non-converted cytosines were discarded via a custom script to remove clustered non-conversion likely caused by RNA secondary structure (3C filter) (Fig S2B, bottom panel). Furthermore, we also excluded those non-converted cytosine positions for which less than 90% of reads passed the 3C filter (S/N90 filter). We then retained “C” positions with ≥20 reads coverage (20RC filter), ≥3 reads supporting the non-conversion, and with a non-conversion ratio ≥10% (10 MM). Sites that passed all filters in at least two biological replicates were called as high-confidence m⁵C sites. Sites passing the filtering steps in only one of the replicates had to satisfy a higher threshold of ≥5 reads supporting non-conversion. After merging called sites from all datasets of non-manipulated NSUN gene expression, this yielded a union set of 6,393 m⁵C sites (Table S2).

Concordance in site calls between samples will be limited by several factors. Expression levels of m⁵C methyltransferase enzymes differ substantially across different tissues and cell lines (e.g., NSUN2/5/6; Fig S3A). Read coverage varies strongly between datasets for many sites (Fig S3B) because of differences in gene expression profile between samples and depth of library sequencing, which strongly reflects on the number of m⁵C sites called in each samples of our union set (Fig S3C). Moreover, most sites show a non-conversion rate just above the empirically chosen ≥10% cut-off (Fig S3E), implying considerable random fluctuation as to whether a site is called in a given sample. Indeed, overlap in union sites called between different samples is modest (Fig S4A). In any given two sample comparison, most m⁵C sites that are not called in common simply either did not reach the ≥20 read coverage threshold or exhibited stoichiometry below threshold (e.g., between 1–10%) in one of the samples (Fig S4B–D). Reassuringly, most union set sites are found in protein coding genes (Fig 1B) and they distribute nonrandomly along mature mRNAs (3,965 sites), showing enrichment around the start codon (Fig 1C), a pattern found repeatedly with single source material in published work (Yang et al, 2017; Chen et al, 2019; Huang et al, 2019; Schumann et al, 2020) and again here (Fig S4E). This establishes utility of our union set as a larger yet still prototypical list of human transcriptomic sites across multiple cell types for further exploration of m⁵C enzymology and function.

The spectrum of mRNA m⁵C writer enzymes

We extended our analyses to several datasets derived from cell lines depleted of NSUN2 (HeLa human cervical cancer cells; [Yang et al, 2017; Huang et al, 2019]) or NSUN6 (HEK293T human embryonic kidney cells; [Liu et al, 2021a]) (Table S3). By re-analysing NSUN2 KO or knockdown data from HeLa cells, we found that 798 m⁵C sites (91.3%) in Huang et al (2019) and 859 (52.8%) in Yang et al (2017) had a reduced non-conversion level upon NUSN2 depletion compared with wild-type control and were therefore regarded as NSUN2-dependent. Similarly we found that 275 (86.7%) sites had reduced methylation level in HEK293T cells upon NSUN6 KO (Liu et al, 2021a). Lastly, upon re-analysis of NSUN2 and NSUN6 double KO data (Liu et al, 2021a) we observed the m⁵C stoichiometry of 44 (8.19%) remained unaffected. In this way, we identified 1,368 NSUN2-mediated, 275 NSUN6-mediated and 44 NSUN2/6 independent m⁵C sites using experimental data. As expected, we found a predominance of sequence and structural features of type I and type II m⁵C sites, respectively (Liu et al, 2021a). Type I sites are characterized by a G-rich triplet downstream of the modified position at the 5′ end of a hairpin (Fig 2A and B—top) (Huang et al, 2019; Schumann et al, 2020), resembling the context of m⁵C sites (C48-50) modified by NSUN2 in the variable loop of multiple tRNAs (Auxilien et al, 2012). Type II sites are adjacent to a 3′ UCCA motif and located in the loop of a hairpin (Fig 2A and B–middle) (Liu et al, 2021a; Selmi et al, 2021). In this way they comply with a mix of sequence and shape-selective substrate recognition by NSUN6, which otherwise targets C72 next to a UCCA motif at the 3′ end of some tRNAs (Long et al, 2016; Liu et al, 2017). Based on the sequence characteristics at position 1 to 5 downstream of these experimentally identified NSUN2 and -6 targets, we generated a position weight matrix (PWM). We then predicted the enzymes responsible for m⁵C sites in our union set (Fig 2C) and individual samples (Fig S3D) using the PWM via FIMO software (Grant et al, 2011). This indicated that in the wider human transcriptome, around two-third of m⁵C sites are likely deposited by NSUN2, whereas around one-sixth are deposited by NSUN6 (Fig 2C) and these sites still comply with type I and type II structural characteristics, respectively (Fig S4F and G). The remaining one-sixth of sites that could not be attributed to either NSUN2 or NSUN6 showed no strong discernible commonalities (Fig 2A and B–bottom; Fig S4F and G–bottom). Given that MTase expression levels will vary between cell and tissue sources (e.g., Fig S3A), we tested for any co-variation of site methylation levels. Indeed, we found a selective positive correlation of the methylation level of NSUN2- or NSUN6-dependent sites with NSUN2 or NSUN6 expression levels, respectively. NSUN2/6-independent sites (“other”) showed no correlation with expression levels of any of the NSUN MTases tested (Fig 2D). This confirms and extends expectations that NSUN2 and -6 are the main m⁵C writers in mRNA but also leaves scope for other enzymes that might deposit m⁵C sites on mRNA.

We next looked at bsRNA-seq data from LN229 human glioblastoma cells, where NSUN5 was either epigenetically silenced (control) or overexpressed by lentiviral transduction (OE) (Janin et al, 2019). NSUN5 is a rRNA:m⁵C MTase responsible for modifying C:3782 in human 28S rRNA and equivalent positions in other organisms (Sharma et al, 2013; Schosserer et al, 2015; Heissenberger et al, 2019; Janin et al, 2019). We found 446 mRNA m⁵C sites in the control but 3,921 sites in NSUN5 OE (Fig 3A), in the context of similar read coverage (Fig S5). Sites in the control condition complied with type I sequence features and were enriched around start codons (Fig 3B). By contrast, in the NSUN5 OE condition, a GUNGCCANNUG motif was found to be prevalent, with broader m⁵C site enrichment along the mRNA coding sequence instead (Fig 3B). This sequence motif clearly resembles the sequence context on human 28S rRNA (C:3782) (Fig 3C), which is highly conserved from yeast to human (Heissenberger et al, 2019). Using the established PWM of NSUN2 and -6 targets (details in the Materials and Methods section), we predicted sites in LN229 cells to be either NSUN2 or -6 dependent (Fig S5). Looking specifically at sites determined to be NSUN5-dependent, based on their increased non-conversion levels upon NSUN5 overexpression (Fig S5, left - “NSUN5-dependent”), we found a narrow but specific pattern of predicted secondary structure (Fig 3D), which is consistent with the predicted structure in Liu et al (2023), and reminiscent of the situation of 28S rRNA C:3782, which is in a relatively unstructured region but immediately downstream of two base-paired positions (Fig 3E).

The NSUN5-targeted sequence motif described above has also just been identified in bsRNA-seq data from early developmental stage samples as well as Nocodazole-treated HeLa cells (Noc-HeLa), and termed “type III”, using a novel computational framework for epitranscriptomic motif discovery (Liu et al, 2023). Nocodazole arrests cells in mitosis when most nuclear content is released into the cytoplasm, providing more opportunity for predominantly nuclear/nucleolar MTases such as NSUN5 to modify mRNA (Liu et al, 2022). Comparison of sites between the two cell lines/conditions was severely limited by differences in gene expression. 42 m⁵C sites were called for NSUN5 in both conditions (Fig 3F). Remarkably, these were all the sites among 2,939 in LN229 OE that had sufficient coverage (≥20RC) also in data from Noc-HeLa. Conversely, whereas around half of the 531 NSUN5 sites called in Noc-HeLa also had coverage in LN229 OE, 42 still represented a substantial overlap. Altogether, this further substantiates that NSUN5 is another MTase capable of depositing m⁵C on mRNA.

RBP binding sites are enriched around m⁵C sites

To investigate potential spatial proximity between m⁵C sites and footprints of known RBPs we obtained enhanced UV cross-linking and immunoprecipitation (eCLIP) sequencing datasets from ENCODE, available for 122 RBPs from the human lymphoblast K562 cell line and for 104 RBPs from human liver cancer HepG2 cells (Van Nostrand et al, 2020). We further used iCLIP data for the known m⁵C reader ALYREF in HeLa cells from CLIPdb as positive control (Yang et al, 2015). It is reported that genes with comparable expression levels in different cell lines show similar RBP enrichment peaks (Van Nostrand et al, 2020), justifying our approach to compare eCLIP and m⁵C data from different sources. Subsequently, we performed enrichment analysis of RBPs binding near m⁵C site deposition in mature mRNAs. After identifying the RBP binding sites as footprint centre positions, we established for each RBP, sets of mRNAs that featured both, m⁵C sites (from the union set or selected subsets) and footprints of said RBP (from both, K562 and HepG2 cells, or individually; see the Materials and Methods section for details). RBP-adjacent regions around each RBP’s footprint centre (±20/30/50/70 nucleotides) were defined and the distribution of m⁵C and unmodified cytosine positions assessed within and outside these regions. This approach found statistically significant m⁵C enrichment in several RBP-adjacent regions that often persisted irrespective of interval width and dataset cross-comparison (Fig S6A–D). We settled on a ±50 nucleotide interval around RBP footprint centres for further analyses. Using the same method, ALYREF was also seen as enriched near m⁵C sites in HeLa cells (Fig S6E), indicating the reliability of this method. To test the robustness of the association between RBP and m⁵C, we compared the enrichment of HepG2 RBP footprints with m⁵C sites in mRNA from different cell contexts, including HepG2-derived (Liu et al, 2021a) and HeLa-derived sites (Yang et al, 2017; Huang et al, 2019; Schumann et al, 2020) (Fig S6E) and the full union set (Fig 4B). When using 708 HepG2-derived m⁵C sites, we identified m⁵C co-enrichment with five RBPs: DDX3X (DEAD-box helicase 3 X-linked), UPF1 (Up-frameshift protein 1), NCBP2 (nuclear cap-binding protein subunit 2), RPS3 (ribosomal protein S3), and DHX30 (DExH-box helicase 30) (Fig 4A). Using the 2,342 m⁵C sites derived from HeLa cells (Fig S6E) or the 3,965 mRNA sites from our union set (Fig 4B) yielded seven RBPs enriched in both comparisons, the five seen before plus FTO (Fat mass and obesity associated) and AKAP1 (A-kinase anchoring protein 1). We also repeated the analysis using eCLIP data from K562 cells and found six of the seven RBPs to co-enrich with m⁵C sites from the union set (lack of data for AKAP1 in K562 cells precluded its assessment; Fig S6F).

We performed several additional analyses to further characterize the m⁵C-RBP associations. First, we plotted footprint counts of RBPs (from Fig 4A; HepG2 eCLIP data) inside and outside a ±50 nucleotide interval around m⁵C sites (HeLa cell data). This further confirmed a positive association of footprints with m⁵C for five RBPs, DDX3X, UPF1, RPS3, NCBP2, and AKAP1 (Fig 4C). Second, we plotted the density of m⁵C and unmodified cytosines relative to the centre of footprints for the five most consistently enriched RBPs (Fig 4D), again showing marked co-enrichment.

Through their recognised functions, the seven enriched RBPs associate m⁵C sites with several mRNA-related processes. For example, NCBP2 promotes mRNA nuclear export together with ALYREF (Kataoka, 2023), a previously known m⁵C reader (Yang et al, 2017). RPS3, AKAP1, DHX30, and DDX3X all have roles in mRNA translation (Dong et al, 2017; Bosco et al, 2021; Ryan & Schröder, 2022; Bohnsack et al, 2023; Cohen et al, 2024). FTO is a demethylase involved in m⁶A metabolism (Mauer et al, 2017). Finally, UPF1 is involved in several mRNA decay pathways (Lavysh & Neu-Yilik, 2020).

UPF1 function is affected by the lack of NSUN2

We chose UPF1 for experimental follow-up to investigate a functional association with m⁵C. siRNA-mediated UPF1 knockdown was performed in both WT and NSUN2 KO HeLa cells (Fig 5A), to compare the effects of the lack of UPF1 with and without NSUN2-mediated mRNA methylation. Details on the generation of NSUN2 KO cell lines are described in Acera Mateos et al (2024). We performed RNA sequencing for three biological replicates of the RNA samples derived from UPF1 knockdown experiment (Fig S7A and B). Compared with the scrambled siRNA (siSCR) control condition, we found that the knockdown of UPF1 in WT and NSUN2 KO cells resulted in a similar number of down-regulated genes (∼2,000). However, the number of up-regulated genes was lower in NSUN2 KO (1,908) compared with WT (2,365) cells (Fig S7C). To investigate whether this differential regulation could be attributed to the different m⁵C content in mRNA in WT versus NSUN2 KO cells, we used bisulfite sequencing data of WT untreated cells and NSUN2 KO/KD HeLa cells (Yang et al, 2017; Huang et al, 2019) to compile a list of 1,633 NSUN2-dependent m⁵C sites in HeLa cells (Table S1). We then selected those m⁵C sites located on mRNA that are also enriched in UPF1 binding and looked at their expression levels upon UPF1 knockdown, in WT versus NSUN2 KO cells. Those mRNA with an m⁵C modification overlapping with UPF1 binding site (“in region”—within ±50 nucleotides) were found to significantly accumulate in WT cells upon UPF1 knockdown; this up-regulation is significantly reduced when UPF1 is knocked down in NSUN2 KO cells (Figs 5B and S7D). Interestingly, this differential regulation between WT and NSUN2 KO cells was either lost, or strongly reduced, for those mRNAs whose m⁵C sites do not overlap with UPF1 binding site (“not in region”—outside ±50 nucleotides) (Fig 5C). As an additional control, we looked at the RNA levels of those mRNAs with UPF1 binding overlapping with NSUN2-independent m⁵C sites. In particular, these sites have the m⁵CUCCA consensus motif associated with NSUN6 methyltransferase and showed no reduction in methylation upon NSUN2 KO/KD (Liu et al, 2021a). As expected, the down-regulation of UPF1 protein had the same regulatory effect on these mRNAs in both WT and NSUN2 KO (Fig S7E). RT-qPCR validation of four mRNAs selected from Fig 5B confirmed a significant accumulation of RNA upon UPF1 knockdown in WT cells and the lack of such regulation upon UPF1 knockdown in NSUN2 KO HeLa cells (Fig 5D).

We then checked the effects of the lack of NSUN2 on the steady state of mRNAs. We divided the transcriptome into five groups (Fig 5E): non-targets (I); NSUN2 targets only (II); UPF1 targets only (III); UPF1 and NSUN2 targets (IV)—same as Fig 5C; UPF1 binding around NSUN2-depentent m⁵C (V)—same as Fig 5B. When comparing the RNA levels of each group upon NSUN2 KO, with (siSCR) and without (siUPF1) UPF1 protein, we found that group V had the highest difference in RNA level relative to group I (Fig 5F, right panel), which recapitulates the effects of UPF1 knockdown and links the effects of UPF1 protein with the methyltransferase activity of NSUN2.

We then decided to assess whether the reduced methylation content in mRNA caused by the KO of NSUN2 also affects the binding ability of UPF1 protein. To test this, we performed UPF1 CLIP experiments in five biological replicates of both WT and NSUN2 KO HeLa cells (Fig 6A). RT-qPCR analysis on the RNA immunoprecipitated with UPF1 protein revealed a significant increase in binding activity for UPF1 upon NSUN2 KO. In particular, the same genes tested in UPF1 KD experiment (Fig 5D) were found to be significantly more enriched in NSUN2 KO samples compared with their WT control (Fig 6B and C).

Collectively, these results indicate that removing nearby m⁵C sites negatively affected UPF1’s capacity to promote target mRNA degradation, when improving its ability to persistently bind to those targets.

Quality control and mapping of RNA BS-seq reads

For each public dataset, raw reads were subjected to FastQC (v0.11.9). Low-quality bases and adaptor sequences were removed using Trimmomatic (v0.38) with options ILLUMINCLIP:Adapter.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW: 4:20 MINLEN:50. We checked the conversion level for each dataset using spike-in sequence (Fig S2B–top, Table S1). Clean reads were mapped on the 45S rRNA, pre-tRNA, and predicted mature tRNA using meRanT tool (align bsRNA-seq reads to reference using Bowtie2 v2.3.5) in MeRanTK (v1.2.1b) with options (‒k 10). Then, paired reads were obtained from unmapped reads using seqkit. Finally, unmapped reads pairs were mapped to the hg38 genome using the meRanGh tool (align bsRNA-seq reads to reference using HISAT2 v2.1.0) in MeRanTK with options (‒fmo). Only uniquely mapped reads were retained and used for m⁵C sites detection.

The library complexity was analysed for each dataset and represented by PCR Bottlenecking Coefficient 1, calculated as the number of genomic locations where exactly one read maps uniquely divided by the number of distinct genomic locations to which read maps uniquely. For low complexity datasets (T24 cells, [Chen et al, 2019]), PCR duplicates were marked by Picard MarkDuplicates (v2.7.1) with the default parameters, and removed via samtools (v 1.9) by filtering 1024 flag value.

m⁵C site calling and annotation for each dataset

The discovery of m⁵C sites from each dataset was performed similarly to Schumann et al (2020). Read coverage at each cytosine position in the genome was obtained using the “mpileup” function in samtools (v1.9). The C-to-T mismatches on forward reads and G-to-A mismatches on reverse reads were regarded as unconverted cytosines, and called using a custom script with parameters “-minBQ 30 ‐‐overhang 6,” where reads were filtered by minimum base quality score 30 and removed terminal 6 nt to avoid overestimation of non-conversion. For ribosome profiling data, data from individual fraction libraries were combined with their respective biological replicate.

Reads containing more than three unconverted cytosines were considered as conversion failure and removed from the bam files (“3C filter”). In the meantime, the RNA icSHAPE values determined by Sun et al (2021) in HeLa, HEK293 and HepG2 cells were downloaded from GEO (GSE145805) to confirm that m⁵C sites flagged by the “3C filter” were in regions with low icSHAPE values, that is, in highly structured regions. Then candidate sites with a signal-to-noise ratio <0.9 (3C/raw; “S/N ≥ 0.9”) were dropped to reduce false positives. To retain high-confidence non-conversion sites, the following criteria were applied: (1) Minimum total read coverage was set as 20 (“20RC”); (2) non-converted C ≥ 3 (“3C”); (3) C + T coverage ≥ 80% (“80CT”); (4) non-conversion of ≥ 10% (“10 MM”); (5). For replicates integration, we select m⁵C candidates detected in at least two replicates, and set non-converted C ≥ 5 for sites detected in only one replicate. Candidate m⁵C sites were annotated to gene locus and transcript types according to the hg38 GENCODE v32 annotation (UCSC) using “intersectBed” in bedtools and mapped to six features simultaneously: 5′UTR, CDS, 3′UTR, ncRNA_exonic, intronic, and intergenic.

Distribution and enrichment of m⁵C along mRNA

Only m⁵C sites on exonic protein-coding transcripts were used in this analysis. mRNAs with m⁵C sites were divided into three segments (5′UTR, CDS, and 3′UTR) and normalised according to their average length. All C positions on the transcripts with m⁵C candidates were selected as background control. The relative position of each candidate site and background C in the corresponding segment were noted. We identified the relative position of m⁵C sites and background C around the start codon and stop codon separately to calculate the spatial enrichment for each bin.

NSUN-enzyme dependence of m⁵C sites

We collected NSUN2 (Yang et al, 2017; Huang et al, 2019) and NSUN6 depletion datasets (Liu et al, 2021a) for NSUN2/6-dependent m⁵C sites analysis as described in Chen et al (2019). In detail, m⁵C sites with more than 0.05 methylation ratio decrease and methylation level reduced to less than 0.1 in NSUN2 or NSUN6 depleted cells were considered to be catalysed by NSUN2 or NSUN6, respectively. In particular, for NSUN2-dependent m⁵C we used the union set of Yang et al (2017); Huang et al (2019). The sequence of position 0 to +5 of NSUN-dependent m⁵C sites were then used to calculate position weight matrix (PWM) of consensus motifs. NSUN-dependence of m⁵C sites in union set can be further predicted by PWM through FIMO (Find Individual Motif Occurrences) in the MEME suite (v 5.4.1). For potential NSUN5-dependent sites, we used NSUN5-overexpression and epigenetically silenced NSUN5 data in LN229 cells (Janin et al, 2019), and regarded sites with 0.05 methylation ratio increase and methylation level higher than 0.1 in NSUN5-OE group as potential NSUN5-dependent sites. All the consensus motifs are plotted by ggseqlogo (Wagih, 2017).

Proximity to RBP binding sites

RBP footprints reported by Van Nostrand et al (2020) were downloaded from ENCODE; ALYREF footprints in HeLa cells line were obtained from CLIPdb database (Yang et al, 2015) as a positive control. For ENCODE datasets, the intersections of two biological replicates with fold-enrichment ≥4 and P ≤ 10⁻³, were filtered as significant peaks and the middle sites of the peaks were regarded as RBP binding sites. Genes with both RBP binding sites and m⁵C sites were considered for enrichment analysis. Bins were divided around RBP binding sites with the number of m⁵C sites and total C within the bins. Then Fisher’s exact test was applied to calculate the enrichment significance. We set up a gradient region as ±20, ±30, ±50 nt, and ±70 nt, to test the influence of bin size. Finally, ±50 nt was used for further analysis. The relative distance of m⁵C sites to RBP binding sites was identified and the background C within the same region to get the distribution of m⁵C around RBP binding sites.

Cell culture and transfection

HeLa cells (human cervical cancer) were obtained from ATTC and confirmed with CellBank Australia. Cells were grown in DMEM medium (Gibco) supplemented with 10% FBS and 1% antibiotic-antimycotic solution (Sigma-Aldrich) and passaged when 70–90% confluent.

For knockdown experiments, 150 × 10³ cells were plated in 35 mm plates and transfected 6–12 h later with the siRNA against the target selected (FlexiTube siRNA UPF-1, Cat. No. GS5976; QIAGEN) or the negative control (non-silencing siRNA, Cat. No. 1022076; QIAGEN) with a final concentration 30 nM, using 5 μl of Lipofectamine RNAiMAX Reagent (Thermo Fisher Scientific) and 300 μl of Opti-MEM (Thermo Fisher Scientific). The medium was replaced 12 h later and cells were harvested 48 h later.

Protein analysis

Cells were harvested with 200–500 μl of Protein Extraction Buffer (50 mM Tris pH 7.5, 5 mM EDTA, 150 mM NaCl, 21.5 mM MgCl₂, 10% glycerol, 1% Triton, 1X PIC [cOmplete, EDTA-free Protease Inhibitor Cocktail; Sigma-Aldrich]) and incubated 10 min on ice, then incubated on a rotator for 30 min at 4°C and centrifuged at 17,000g for 10 min at 4°C.

The supernatant was transferred to a clean tube, used, or stored at −80°C. Total protein concentration was measured through the Qubit Protein Assay Kit (Thermo Fisher Scientific) following manufacturer’s instructions. 30 μg of proteins were loaded on NuPage 4–12% Bis-Tris Protein Gels (Invitrogen) followed by transfer onto PVDF membrane. The membrane was blocked in Odyssey Blocking Buffer (for IR-Dye detection; LI-COR 927-40000) and probed with primary antibodies: (anti-NSUN2 (1:1,000, 20854-1-AP; Proteintech), anti-NSUN5 (1:1,000, 15449-1-AP; Proteintech), anti-UPF1 (1:1,000, ab86057; Abcam), anti-ACTB (1:1,000, sc-4778 AF790; SantaCruz). The membranes were probed with a secondary antibody, either anti-mouse-IR-Dye800 (1:10,000, 926-32210; LI-COR) or anti-rabbit-IR-Dye680 (1:10,000, 925-68071; LI-COR), and imaged using the Odyssey CLx Imaging System (LI-COR).

RNA analysis

Extraction of total RNA was performed using the Direct-zol RNA Miniprep (Zymo Research) kit according to the manufacturer’s instructions. Reverse transcription reactions were performed with PrimeScript RT Master Mix (Takara Bio) on 0.5–1 μg of total RNA in a 10 μl reaction. RT-qPCR analyses were performed with 20 ng equivalent of cDNA, 5 μl of 2X SYBR Mastermix (QIAGEN), 1 μl of 5 μM primers, and water to a final volume of 10 μl. DNA amplification (followed by melting curve analysis) was monitored with a QuantStudio 12k Flex Realtime PCR instrument. Primer sequences used for RT-qPCR are listed in Table S4.

Relative RNA quantity was calculated as the fold change (ΔΔCt) with respect to the experimental control sample set as 1 and normalised over ACTB. The ratio of each sample versus its experimental control was tested by two-tailed t test.

UPF1 knockdown RNA-seq and differential expression analysis

The raw reads were subjected to FastQC (v0.11.9). Low-quality bases and adaptor sequences were removed using Trimmomatic (v0.38) with options (ILLUMINCLIP:Adapter.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW: 4:20 MINLEN:50). Clean reads were mapped to 4S rRNA using bowtie2 (v2.3.5) with parameters (“-q ‐‐sensitive ‐‐reorder ‐‐no-unal ‐‐un-conc-gz”). Unmapped reads were mapped to hg38 genome using HISAT2 (v2.1.0) with parameters (“‐‐no-softclip ‐‐score-min L,-16,0 ‐‐mp 7,7 ‐‐rfg 0,7 ‐‐rdg 0,7 ‐‐max-seeds 20 -k5 ‐‐dta”). The mapped reads were processed to featureCounts (v2.0.1) for feature counting.

To quantify RNA, fragment abundance in genes with ≥50 mapped reads were selected and normalised using the size factors estimated by the median of all genes implemented in the DESeq2 (v1.30.1) Bioconductor package. Differential expression analysis was performed by DESeq2. Genes with log₂(fold change) ≥ 1.2 and adjusted P-value ≤ 0.05 were regarded as differentially expressed genes.

Cross-linking immunoprecipitation

One million HeLa cells were washed twice with 1X PBS and irradiated once at 150 mJ·cm⁻² at 254 nm using the Stratalinker UV cross-linker. Cells were then scraped in CLIP Lysis Buffer (50 mM Tris–HCl, pH 7.4, 100 mM NaCl, 1% Igepal CA-630, 0.1% SDS, 0.5% sodium deoxycholate, 1X PIC [cOmplete, EDTA-free Protease Inhibitor Cocktail; Sigma-Aldrich], and RNasin Plus Ribonuclease Inhibitor [Promega]), incubated on ice for 5 min and passed through a 21G needle. Lysates were treated with 3 μl of TURBO DNase (Thermo Fisher Scientific) for 3 min at 37°C while shaking at 1,100 rpm using a Thermomixer R (Eppendorf). Samples were then centrifuged for 10 min at 17,000g at 4°C to clear the lysate. The supernatant was collected and quantified using Qubit Protein Assay Kit (Thermo Fisher Scientific) following manufacturer’s instructions.

50 μl of Dynabeads Protein G (Thermo Fisher Scientific) were previously washed twice in CLIP Lysis Buffer and then resuspended in the same buffer with either 9 μg of anti-UPF1 goat polyclonal antibody (Cat. No. A300-038A), or 9 μg of Abcam mouse IgG2a Isotype Control antibody (Cat. No. ab18413) in a final volume of 100 μl and incubated in rotation at room temperature for 1 h for antibody conjugation. The conjugated beads were then incubated overnight on a rotator at 4°C with 500 μg of lysate. 50 μg of lysates was also rotated to be used as inputs. The next day beads were washed four times in High-Salt Buffer (50 mM Tris–HCl, pH 7.4, 1 M NaCl, 1 mM EDTA, 1% Igepal CA-630, 0.1% SDS, 0.5% sodium deoxycholate) and twice in PK buffer (100 mM Tris–HCl pH 7.4, 50 mM NaCl, 10 mM EDTA). Beads were resuspended in 50 μl of PK buffer, of which 10 μl were used for Western Blotting analysis and 40 μl were digested with Proteinase K Solution (Cat. No. 25530049; Thermo Fisher Scientific) while shaking at 1,100 rpm in a Thermomixer R (Eppendorf) for 30 min at 37°C for RNA analysis. RNA was purified by phenol/chloroform extraction and precipitated for RT-qPCR analysis.