14 publicly available CRC bulk transcriptomic cohorts comprising tumor and adjacent normal tissues with clinical and prognostic information were integrated (Supplementary Table S1). Cohorts with sample size <50, non-original data, or incomplete clinical annotation were excluded. The Cancer Genome Atlas Colon Adenocarcinoma (TCGA-COAD) (n = 491) was obtained from the GDC portal (https://gdc.cancer.gov/), and 13 GEO datasets (GSE103479, GSE14333, GSE17536, GSE18105, GSE21510, GSE25010, GSE28702, GSE29621, GSE29638, GSE39084, GSE71187, GSE72970, and GSE87211) were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/) (Cancer Genome Atlas Network, 2012; Allen et al., 2018; Jorissen et al., 2009; Smith et al., 2010; Matsuyama et al., 2010; Tsukamoto et al., 2011; Ågesen et al., 2011; Tsuji et al., 2012; Ågesen et al., 2012; Bruun et al., 2018; Kirzin et al., 2014; An et al., 2015; Cherradi et al., 2019; Hu et al., 2018). Among these, a prognostic meta-validation cohort comprising seven datasets with three types of clinical outcome records (total n = 1,550) was established. The remaining datasets were applied to analyses of clinical subgroup differences, pathway enrichment, and drug sensitivity. For cross-cohort comparison of target gene expression, raw data were normalized using log2(count + 1) transformation.

To elucidate the tumor biological characteristics of CRC, this study integrated 54 tumor and normal colon samples from four datasets (GSE166555, GSE299427, GSE232525, GSE200997) (Uhlitz et al., 2021; Wong et al., 2025; Yang et al., 2023; Khaliq et al., 2022) (Supplementary Table S2). Cell annotation was performed using a multi-tiered prioritization approach: (i) adopting original labels when available (GSE166555, GSE299427); (ii) manual annotation based on reference cell markers supplied within the original datasets (GSE232525, Supplementary Table S3); and (iii) confirming epithelial and stromal populations using established CRC epithelial markers (GSE200997, Supplementary Table S4). Data normalization, variable gene selection, and principal component analysis (PCA) were performed using the “Seurat” R package (Hao et al., 2024). Neighborhood graphs were constructed from PCA results, followed by Louvain clustering (resolution 0.3–0.5) and UMAP visualization. Batch effects were corrected using the “Harmony” algorithm.

The spatial dataset GSE225857 was obtained from GEO (Wang et al., 2023) (Supplementary Table S2). Data normalization was performed using the “Sctransform” method to reduce technical artifacts while preserving biological heterogeneity. Downstream analyses involved clustering procedures analogous to those used in single-cell omics. Finally, spot annotation was guided by pathologist-provided histological labels, and spatial gene expression patterns were visualized using “SpatialFeaturePlo”.

Differential expression analysis was performed using the “limma” algorithm. DEGs were defined as those with adjusted p < 0.05 and |log2FC| > 0.5. Ranked gene lists from each dataset were integrated using the ‘RobustRankAggreg’ (RRA) algorithm to identify genes consistently differentially expressed across cohorts (Kolde et al., 2012).

The “survival” package was utilized to calculate the hazard ratio (HR), 95% confidence interval (CI), log-rank p-value, and Kaplan–Meier (KM) curves. Meta-analysis was conducted using the “MIME” package, applying fixed- or random-effects models with the DerSimonian–Laird estimator to account for heterogeneity (Liu et al., 2024). Effect sizes were pooled via inverse-variance weighting. Heterogeneity was quantified using Q, I², and H statistics, and publication bias was assessed via Egger’s test.

Gene Set Enrichment Analysis (GSEA), Gene Set Variation Analysis (GSVA), and single-sample GSEA (ssGSEA) were performed using “clusterProfiler” and “gsva” R packages (Yu et al., 2012). For bulk data, non-significant genes (p > 0.05) were removed prior to GSEA. GSVA and ssGSEA were used to evaluate Hallmark gene sets and pathway activity at the single-cell level.

Drug sensitivity analysis was inferred using the “oncoPredict” package (bulk transcriptomics) and the “beyondcell” package (spatial transcriptomics) (Fustero-Torre et al., 2021). Specifically, for bulk samples, drug screening data were obtained from the Cancer Therapeutics Response Portal (CTRP) and the Genomics of Drug Sensitivity in Cancer (GDSC). Linear regression models were built for each drug and applied to CRC expression matrices to predict the half-maximal inhibitory concentration (IC50) values. For spatial transcriptomic samples, the “beyondcell” algorithm was applied to identify cells or spots with distinct drug response profiles. A Beyondcell Score (BCS) ranging from 0 to 1 was calculated for each drug–cell/spot pair, with higher scores indicating greater sensitivity to the specific drug.

Additionally, the Subclass Mapping (SubMap) algorithm was used to compare patient expression profiles with reference datasets of immune therapy responses (Hoshida et al., 2007). Samples were stratified based on either target gene expression patterns or immune response phenotypes to assess commonalities between subgroups. This analysis supported its utility in predicting differential responses to immunotherapy across distinct subgroups.

All statistical analyses and data visualization were performed using R (version 4.3.2) or Python (version 3.13). Continuous variable correlations were evaluated using Pearson’s correlation. Group comparisons were performed using the Wilcoxon rank-sum test. Classifier performance was assessed by Receiver Operating Characteristic (ROC) analysis, with area under the curve (AUC) indicating discriminatory ability. All statistical tests were two-sided, with p < 0.05 considered significant.

Human CRC cell lines SW480 and HCT116 and human embryonic kidney HEK293T cells were obtained from the American Type Culture Collection (ATCC, United States). Cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; Gibco, Thermo Fisher Scientific, Carlsbad, CA) supplemented with 10% fetal bovine serum (FBS; Gibco) at 37 °C in a humidified atmosphere with 5% CO₂. Cells were routinely passaged at 70%–80% confluence using trypsin-EDTA and were used for experiments during the logarithmic growth phase. Unless otherwise specified, cells were seeded 12–24 h before treatment/transfection to achieve ∼50–70% confluence at the start of experiments. Oxaliplatin (L-OHP), erastin, Ferrostatin-1 (Fer-1), RSL3, Brequinar (BQR), and Actinomycin D were purchased from MedChemExpress (MCE, NJ, United States). Stock solutions were prepared according to the manufacturers’ instructions and stored aliquoted at −20 °C or −80 °C to avoid repeated freeze–thaw cycles. The working concentrations and treatment durations for each assay are provided in the corresponding figure legends.

Primary CRC tissues and paired adjacent normal tissues were collected from patients at the Fourth Affiliated Hospital, Zhejiang University School of Medicine (2024–2025). Diagnosis was confirmed independently by two pathologists. Written informed consent was obtained from all participants. The study was approved by the Institutional Ethics Committee (K2025210) and conducted in accordance with the Declaration of Helsinki.

Lentiviral vectors encoding two independent NSUN2-targeting shRNAs (shNSUN2-1 and shNSUN2-2) or a scrambled control (shCtrl) were purchased from GeneChem (Shanghai, China). The target sequences were: shNSUN2-1: 5′-GAT​CGG​AGC​GAT​GCC​TTA​GGA​TAT​TAC​TCG​AGT​AAT​ATC​CTA​AGG​CAT​CGC​TCT​TTT​TT-3′; shNSUN2-2: 5′-GAT​CGT​GCA​GTG​TCC​CAT​CGT​CTT​ATC​TCG​AGA​TAA​GAC​GAT​GGG​ACA​CTG​CAT​TTT​TT-3′. Lentiviruses were generated in HEK293T cells using psPAX2 and pMD2.G packaging plasmids, and viral supernatants were harvested 48 h post-transfection. Target cells were infected at a multiplicity of infection (MOI) of 10–20, followed by puromycin selection (2 μg/mL, 7 days). NSUN2 siRNAs and negative control siRNA were obtained from GenePharma (Shanghai, China), with the following sequences: siNSUN2-1: forward 5′-GAG​CGA​UGC​CUU​AGG​AUA​UUA​tt-3′ and reverse 5′-UAA​UAU​CCU​AAG​GCA​UCG​CUC​tt-3′; siNSUN2-2: forward 5′-UGA​GAA​GAU​GAA​GGU​UAU​UAA​tt-3′ and reverse 5′-UUA​AUA​ACC​UUC​AUC​UUC​UCA​tt-3′. Constructs for NSUN2 wild-type (WT), enzymatic-dead mutants (C271A, C321A), and DHODH-WT were purchased from GenePharma. Transfections were performed using jetPRIME (Polyplus-transfection, Illkirch, France) according to the manufacturer’s protocol. Unless otherwise specified, cells were harvested for downstream analyses 24–72 h after transfection depending on the assay.

For wound healing assays, cells were seeded in 12-well plates and cultured to a confluent monolayer. A linear wound was created using a sterile 200-µL pipette tip, followed by gentle PBS washing to remove detached cells. The medium was then replaced with serum-free DMEM to minimize the confounding effect of cell proliferation, and images were acquired at 0 h and 24 h or 30 h at pre-marked positions using an inverted microscope. Wound closure was quantified by measuring the wound area in ImageJ, and migration was expressed as the percentage of wound closure relative to baseline (0 h).

For transwell invasion assays, Matrigel-coated transwell inserts (8-µm pore size; Corning, United States) were used. Briefly, 2 × 10⁴ cells were seeded into the upper chamber in serum-free medium, whereas complete medium containing serum was added to the lower chamber as a chemoattractant. After 24 h, non-invaded cells were removed from the upper surface with a cotton swab. Invaded cells on the lower surface were fixed, stained with crystal violet, imaged, and counted under a microscope. Quantification was performed using representative fields per insert, and values were averaged across biological replicates.

Viability was determined using the Cell Counting Kit-8 (CCK-8, Beyotime, China). Cells were seeded into 96-well plates and allowed to attach overnight. After the indicated treatments, CCK-8 reagent was added to each well and incubated for 2 h at 37 °C. Absorbance was measured at 450 nm using a microplate reader. Background values from medium-only wells were subtracted, and viability was normalized to the corresponding control group.

Cells were seeded into 6-well plates at low density and cultured for 10–14 days until visible colonies formed, with medium refreshed every 2–3 days. Colonies were fixed and stained with crystal violet. Colonies containing more than 50 cells were counted using ImageJ, and colony numbers were normalized to the control group.

Apoptosis and lipid peroxidation were quantified by flow cytometry. For apoptosis analysis, cells were stained using an Annexin V-FITC/PI apoptosis detection kit (Meilunbio, China) and analyzed on a CytoFLEX S flow cytometer (Beckman Coulter, United States). Apoptotic cells were defined as Annexin V–positive populations. Compensation and gating were set using unstained and single-stained controls and applied uniformly across samples. For lipid peroxidation, lipid reactive oxygen species (lipid ROS) were assessed using BODIPY 581/591 C11 (Beyotime, China). Cells were incubated with the dye as recommended, washed, and analyzed.

Malondialdehyde (MDA) levels were determined using a Lipid Peroxidation MDA Assay Kit (Beyotime, China). Cells or frozen tumor tissues were lysed or homogenized in Beyotime Western and IP cell lysis buffer (Beyotime, China), followed by centrifugation to collect supernatants for MDA determination. Total protein concentrations were measured using a BCA assay (Beyotime, China) for normalization, and MDA levels were expressed as nmol/mg protein.

TEM was used to examine mitochondrial morphology in cells treated with DMSO or erastin. After ethanol dehydration and resin embedding, samples were sectioned, stained and observed.

RNA was extracted using TRIzol (Invitrogen, United States) and reverse-transcribed with HiScript III (Vazyme, China). qRT-PCR was performed with SYBR Green (Vazyme) on a CFX-96 amplifier (Bio-Rad, United States). GAPDH was used as an internal control, and expression was calculated by the 2^−ΔΔCT method. NSUN2 primers were 5′-ACA​GCC​ACT​GAG​TTG​GTA​TCC-3′ (sense) and 5′-TTA​TGA​TGA​GGC​CGC​ACG​TT-3′ (antisense). DHODH primers were 5′-CCA​CGG​GAG​ATG​AGC​GTT​TC-3′ (sense) and 5′-CAG​GGA​GGT​GAA​GCG​AAC​A-3′ (antisense). GPX4 primers were 5′-GAG​CCA​GGG​AGT​AAC​GAA​GAG-3′ (sense) and 5′-TGG​TGA​AGT​TCC​ACT​TGA​TGG-3′ (antisense). GAPDH was used as an internal control for normalization of the results. GAPDH primers were 5′-AGG​TCG​GTG​TGA​ACG​GAT​TTG-3′ (sense) and 5′-GTA​GAC​CAT​GTA​GTT​GAG​GTC​A-3′ (antisense).

Proteins were extracted using RIPA lysis buffer supplemented with PMSF (Beyotime, China), separated by SDS-PAGE, and transferred onto PVDF membranes (Millipore, United States). Membranes were blocked with 5% non-fat milk dissolved in TBST at room temperature for 1 h, and then incubated with primary antibodies against NSUN2 (20854-1-AP, Proteintech), GPX4 (ab125066, Abcam), DHODH (ab24689, Abcam), and β-actin (66009-1-Ig, Proteintech) at 4 °C for 12 h. After washing with TBST, membranes were incubated with HRP-conjugated secondary antibodies (EpiZyme) at room temperature for 1 h. Bands were visualized using an ECL detection system (Bio-Rad) and quantified by integrated density using ImageJ, with β-actin used as a loading control.

For dot blot, total RNA was denatured at 95 °C for 10 min in a metal bath and immediately chilled on ice. Denatured RNA was then spotted onto nylon membranes, followed by UV crosslinking at 254 nm for 1 h. Membranes were blocked with 5% non-fat milk in TBST at room temperature for 1 h, and incubated with an anti-m5C antibody (68301-1-Ig, Proteintech) at 4 °C for 12 h. After washing with TBST, membranes were incubated with HRP-conjugated secondary antibodies at room temperature for 1 h. Methylene blue staining was used to verify equal RNA loading, and signals were detected using an ECL detection system (Bio-Rad).

For RNA stability, cells were treated with 5 μg/mL actinomycin D (MCE) and harvested at 0–3 h for qRT-PCR analysis of DHODH mRNA decay.

All animal procedures were approved by the Animal Experimental Ethics Committee of Zhejiang University (ZJU20250668). Four-week-old male BALB/c nude mice (GemPharmatech, Jiangsu, China) were maintained under specific-pathogen-free (SPF) conditions. For the subcutaneous xenograft model, SW480 cells (5 × 10⁶) in 0.1 mL PBS were injected subcutaneously into the right flank. Tumors were measured every 2 days, and volume was calculated as L × S²/2. For growth validation, mice were randomized into shCtrl and shNSUN2 groups (n = 5 per group); mice were euthanized at a predefined endpoint (tumor volume of 2,000 mm³), and tumors were excised, photographed, and weighed. For treatment evaluation, mice were randomly assigned to five groups (n = 5 per group) before inoculation and injected with shCtrl- or shNSUN2-transduced SW480 cells according to group assignment. Oxaliplatin (OXA) was dissolved in sterile 5% glucose solution, and imidazole ketone erastin (IKE; MCE, NJ, United States) was prepared in 10% DMSO vehicle. Treatment was initiated for each mouse when its tumor volume reached approximately 100 mm³, with the first dosing day defined as treatment day 0. The groups were as follows: vehicle control (shCtrl tumors; i.p. vehicle), oxaliplatin (shCtrl tumors; OXA 5 mg/kg i.p. every 3 days), NSUN2 knockdown (shNSUN2 tumors; i.p. vehicle), NSUN2 knockdown in combination with oxaliplatin (shNSUN2 tumors; OXA 5 mg/kg i.p. every 3 days), and oxaliplatin in combination with IKE (shCtrl tumors; OXA 5 mg/kg i.p. every 3 days and IKE 30 mg/kg i.p. daily). Vehicle injections were matched to the corresponding treatment schedule. After 14 days of treatment, mice were euthanized and tumors were collected for imaging and weighing; MDA was measured using a Lipid Peroxidation MDA Assay Kit (Beyotime, China) according to the manufacturer’s instructions.

Data are presented as mean ± SD. Comparisons were performed using Student’s t-test, one-way ANOVA or two-way ANOVA with multiple-comparisons correction where applicable (GraphPad Prism 9). A p-value <0.05 was considered statistically significant (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).

To explore the role of m5C modification in CRC, the expression of major m5C-related genes was examined across multi-omics datasets (Supplementary Table S5). Most genes exhibited differential expression between CRC tumor tissues and matched normal tissues. Robust rank aggregation (RRA) analysis identified NSUN2 as the top-ranked gene, displaying consistent overexpression with marked fold changes across multiple cohorts (Figures 1A,B; Supplementary Table S6). Single-cell transcriptomics further demonstrated preferential expression of NSUN2 within epithelial cell clusters of tumor tissues, with higher expression frequency and abundance compared with other m5C-related genes (Supplementary Figure S1A; Figure 1C). Based on pathologist-annotated spatial transcriptomic data, NSUN2 showed higher expression in tumor regions than in adjacent peritumoral regions (Figure 1D; Supplementary Figure S1B), and was further upregulated in liver metastatic tumor regions compared with primary tumor regions (Wang et al., 2023) (Supplementary Figure S1C). Collectively, multi-omics analyses indicate that NSUN2 exhibits both elevated expression and notable specificity in CRC.

To assess its prognostic value, a meta-validation cohort of seven studies (n = 1,550) was constructed. A meta-analysis of Cox regression results demonstrated that NSUN2 serves as a risk factor for poor prognosis in CRC patients, with consistent effects observed across all included studies (meta-HR [95% CI]: 1.72 [1.46–2.03]) (Figure 1E; Supplementary Table S7). Importantly, after adjustment for established clinicopathological variables, including AJCC/Dukes stage, KRAS mutation status, and chemotherapy response where available, NSUN2 remained an adverse prognostic factor across cohorts (Supplementary Figure S1D). Kaplan-Meier analysis further confirmed shorter overall survival (OS) and disease-free survival (DFS) in NSUN2-high groups (Figure 1F; Supplementary Figure S1E). Notably, NSUN2 was enriched in poorly differentiated tumors, associated with PI3K mutations, and exhibited a stage-dependent gradient across the normal–adenoma–carcinoma sequence (Figure 1G). Consistent with a more aggressive clinical phenotype, NSUN2 expression was higher in metastatic lesions than in primary tumors (Supplementary Figure S1F) and was significantly associated with greater local invasion depth (T), lymph node involvement (N), distant metastatic status (M), and overall stage in TCGA (Supplementary Figures S1G–J), collectively supporting a clinical association between NSUN2 expression and tumor invasion and metastasis. Elevated NSUN2 protein expression was also validated in freshly collected CRC tissues (Figure 1H). These results suggest that NSUN2 is frequently overexpressed and serves as an unfavorable prognostic biomarker in CRC.

To examine the functional role of NSUN2, SW480 and HCT116 cells were transduced with lentiviral shRNAs, resulting in >90% knockdown efficiency as confirmed by Western blotting (Figure 2A). NSUN2 depletion significantly impaired colony formation (Figure 2B), migration (Figure 2C), and invasion (Figure 2D). In xenograft models, tumors derived from NSUN2-deficient SW480 cells showed reduced growth and weight compared with controls (Figure 2E). Collectively, these findings indicate that targeting NSUN2 impairs CRC progression. Transcriptomic profiling of NSUN2-high versus NSUN2-low tumors revealed significant enrichment of proliferation-associated programs, including MYC targets, E2F targets, and the G2M checkpoint, together with invasion-related signatures such as epithelial–mesenchymal transition, angiogenesis, and TGF-β signaling (Figure 2F; Supplementary Figure S2; Supplementary Table S8). In single-cell transcriptomic analyses, cancer cell clusters were stratified into NSUN2-positive and NSUN2-negative populations (Figure 2G). GSVA revealed that cancer cell clusters stratified by NSUN2 expression and chemotherapy exposure exhibited distinct enrichment patterns of tumor-promoting pathways. Notably, NSUN2-positive cells were preferentially enriched in oncogenic pathways, including cell cycle regulation and DNA damage repair, consistent with the bulk transcriptomic GSEA results (Figure 2H).

Given that NSUN2 is primarily expressed in malignant epithelial cells and is positively correlated with poor prognosis, we further investigated its potential association with oxaliplatin resistance. In GSE72970, progressive disease (PD) cases displayed significantly higher NSUN2 levels than responders (Figure 3A). ROC analysis demonstrated robust performance of NSUN2 in predicting PD, with an AUC of 0.734 (Figure 3B). Consistently, across both the GSE72970 and GSE28702 cohorts, non-responders to oxaliplatin-based chemotherapy showed elevated NSUN2 expression relative to responders (Figure 3C). The results of the “SubMap” algorithm indicate that there are statistically significant commonalities between the groups based on the expression level of NSUN2 and the groups based on chemotherapy response (Figure 3D). We next assessed the association between NSUN2 expression and oxaliplatin sensitivity across multiple CRC bulk transcriptomic cohorts. NSUN2 expression was positively correlated with oxaliplatin IC50 values, and tumors with high NSUN2 expression consistently exhibited reduced chemosensitivity (Figure 3E; Supplementary Tables S9, S10). In the GSE14333 cohort, NSUN2 expression was positively correlated with oxaliplatin IC50 values (Oxaliplatin_1806: R = 0.33, p = 1.1 × 10⁻⁶; Oxaliplatin_1089: R = 0.42, p = 1.7 × 10⁻¹¹) (Figure 3F).

To experimentally validate these findings, we assessed apoptosis in CRC cells following oxaliplatin treatment. Due to GFP interference from shRNA vectors in FITC-based assays, siRNA-mediated NSUN2 knockdown was employed, resulting in over 70% reduction in protein expression (Figure 3G). Flow cytometry revealed that siRNA-mediated knockdown of NSUN2 enhanced oxaliplatin-induced apoptosis (Figures 3H,I), and CCK-8 assays showed that NSUN2 depletion markedly increased cellular sensitivity to oxaliplatin (Figure 3J). Collectively, these findings demonstrate that high NSUN2 expression is associated with CRC progression and oxaliplatin resistance, whereas its inhibition suppresses malignant phenotypes and sensitizes tumor cells to treatment.

Building on the KEGG enrichment results (Figure 4A; Supplementary Table S8) and the observed association with oxaliplatin sensitivity, we next investigated whether NSUN2 promotes CRC progression and therapy resistance through ferroptosis-related mechanisms. In oxaliplatin-treated single-cell datasets, NSUN2-positive cells were enriched in resistant clusters (Figure 4B). Concurrently, ssGSEA revealed ferroptosis-suppressive activity in NSUN2-positive cell clusters, whereas NSUN2-negative clusters exhibited ferroptosis-promoting activity. These differences were more pronounced after oxaliplatin treatment (Figure 4C; Supplementary Tables S11, S12). Spatial transcriptomics further demonstrated a significant negative spatial correlation between NSUN2 expression levels and sensitivity to classical ferroptosis inducers (such as erastin and ML162) in independent CRC samples (Figures 4D,E; Supplementary Table S13). This finding was corroborated in multiple bulk transcriptomic CRC datasets, where NSUN2 expression showed a significant positive correlation with the IC50 values of various ferroptosis inducers (Figures 4F,G; Supplementary Tables S9, S10). Collectively, these multi-omics analyses indicate that elevated NSUN2 expression is associated with reduced ferroptosis sensitivity.

As NSUN2 functions as an RNA m5C methyltransferase, we next evaluated whether ferroptosis sensitivity was influenced by RNA methylation status. Dot blot analysis showed increased global RNA m5C levels after erastin treatment (Figure 4H). Pharmacological inhibition with Ferrostatin-1 (Fer-1) significantly attenuated both oxaliplatin- and erastin-induced cytotoxicity in CRC cells, indicating that ferroptosis constitutes a major cell death pathway in this setting (Figures 4I,J). NSUN2 knockdown markedly sensitized CRC cells to ferroptosis, evidenced by reduced viability (Figure 4K), elevated lipid ROS (Figure 4L), and increased malondialdehyde (MDA) production (Figure 4M). Transmission electron microscopy (TEM) also revealed classical ultrastructural changes observed during ferroptosis—including condensed mitochondrial membranes, cristae disruption, and outer membrane rupture—in NSUN2-deficient CRC cells, both at baseline and following erastin treatment (Figure 4N). To further validate the in vivo relevance of NSUN2-mediated ferroptosis resistance and oxaliplatin responsiveness, we performed nude mouse xenograft experiments evaluating oxaliplatin alone and oxaliplatin combined with the ferroptosis inducer imidazole ketone erastin (IKE). Tumor growth was more effectively suppressed by oxaliplatin in NSUN2-depleted xenografts than by oxaliplatin alone (Figures 4O,P), supporting in vivo chemosensitization upon NSUN2 depletion. In parallel, oxaliplatin combined with IKE further reduced tumor burden compared with oxaliplatin monotherapy (Figures 4O,P), highlighting the potential translational value of ferroptosis-promoting strategies to enhance oxaliplatin response. Mechanistically, tumor malondialdehyde (MDA) levels were increased in oxaliplatin-treated NSUN2-depleted tumors and in tumors treated with oxaliplatin combined with IKE, relative to tumors treated with oxaliplatin alone (Figure 4Q). These findings support enhanced lipid peroxidation and elevated ferroptotic stress in vivo.

Collectively, these data demonstrate that NSUN2 inhibition enhances ferroptosis sensitivity in CRC cells, whereas elevated NSUN2 expression is linked to ferroptosis resistance and reduced oxaliplatin efficacy.

To further dissect the molecular basis underlying NSUN2-mediated ferroptosis resistance, we systematically interrogated the expression correlation between NSUN2 and major ferroptosis-related genes across multiple CRC bulk transcriptomic cohorts. Meta-cohort analyses identified dihydroorotate dehydrogenase (DHODH) as the ferroptosis suppressor most consistently associated with NSUN2 expression (Figure 5A; Supplementary Tables S14, S15). While NSUN2 expression showed little association with GPX4 or other ferroptosis-related genes, DHODH expression was markedly elevated in NSUN2-high CRC samples, suggesting a unique regulatory relationship. At the single-cell level, DHODH expression was preferentially enriched in malignant epithelial cells with high NSUN2 expression, whereas low-NSUN2 populations exhibited weaker DHODH expression (Figure 5B). Furthermore, joint visualization of NSUN2 and DHODH co-expression patterns in spatial transcriptomic samples revealed significant spatial co-localization within specific tumor subclusters, implying overlapping enrichment in therapy-resistant epithelial regions (Figure 5C). These findings pinpoint DHODH as a putative downstream effector selectively regulated by NSUN2.

Functionally, NSUN2 knockdown sensitized cells to RSL3, a GPX4 inhibitor, but conferred relative resistance to the DHODH inhibitor Brequinar (BQR) (Figures 5D,E). This indicates that NSUN2 differentially modulates GPX4-and DHODH-mediated ferroptosis. qRT-PCR and WB further showed that NSUN2 knockdown markedly suppressed DHODH expression at both mRNA and protein levels, while GPX4 expression remained unaffected (Figures 5F,G). Mechanistically, only wild-type NSUN2, but not catalytically inactive mutants (C271A, C321A), restored DHODH expression in NSUN2-deficient CRC cells (Figure 5H). Global RNA m5C levels were reduced by NSUN2 knockdown and rescued by wild-type NSUN2 (Figure 5I). RNA decay assays further showed accelerated DHODH mRNA degradation in the absence of NSUN2 (Figure 5J). Collectively, these results indicate that NSUN2 stabilizes DHODH mRNA through m5C modification.

To determine whether DHODH mediates the ferroptosis resistance and oxaliplatin tolerance conferred by NSUN2, we ectopically expressed DHODH in NSUN2-silenced SW480 and HCT116 cells. Western blotting confirmed efficient restoration of DHODH expression following its overexpression (Figure 6A). Functionally, colony formation and CCK-8 assays demonstrated that NSUN2 knockdown markedly enhanced erastin-induced growth inhibition, whereas DHODH overexpression largely reversed these effects (Figures 6B,C). Consistent with this, flow cytometry analysis revealed that DHODH overexpression attenuated oxaliplatin-induced apoptosis in NSUN2-deficient cells (Figure 6D). Drug sensitivity assays further indicated that DHODH overexpression counteracted the enhanced susceptibility of NSUN2-silenced cells to both erastin and oxaliplatin (Figures 6E,F). Mechanistically, DHODH overexpression markedly attenuated MDA production and lipid ROS accumulation induced by NSUN2 knockdown (Figures 6G,H), suggesting that the NSUN2–DHODH axis promotes ferroptosis resistance by suppressing lipid peroxidation.

To assess the clinical relevance, we next interrogated public CRC transcriptomic datasets. High DHODH expression was associated with poor chemotherapy response in GSE72970 (Figure 7A). Meta-analysis of Cox regression across multiple CRC bulk transcriptomic datasets confirmed its prognostic significance (Figure 7B; Supplementary Table S16), and Kaplan-Meier analyses of GSE14333 and GSE87211 demonstrated worse survival in DHODH-high patients (Figures 7C,D). Importantly, after adjustment for tumor stage, KRAS mutation status, and chemotherapy response where available, DHODH remained an independent adverse prognostic factor in multivariable Cox analyses (Figure 7E). Together, these results establish the NSUN2–DHODH axis as a mediator of ferroptosis suppression and oxaliplatin resistance in CRC.