Human colon cancer cell lines, RKO and LOVO, were cultured in RPMI-1640 medium (BioSharp Life Sciences) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) and 1% penicillin-streptomycin (Gibco; Thermo Fisher Scientific, Inc.) at 37°C in a humidified atmosphere containing 5% CO₂. Small interfering (si)RNAs targeting METTL3, YTHDC1 and STK25 were synthesized by Shanghai GenePharma Co., Ltd. To ensure efficient gene silencing and reduce potential off-target effects, different siRNA sequences targeting the same gene were evaluated in the two cell lines, and the siRNA exhibiting the highest knockdown efficiency in each cell line was used for subsequent experiments. Transient transfections were performed using siRNA-Mate Plus transfection reagent (Shanghai GenePharma Co., Ltd.) according to the manufacturer's instructions. The siRNA premix was prepared at a concentration of 1.5 pmol/µl (~1.5 µM). Cells were transfected at 37°C in a 5% CO₂ incubator, and transfection efficiency was assessed 48 h post-transfection. RNA extraction and subsequent RNA-based experiments were performed 48 h after transfection. The siRNA sequences used in the present study were as follows (5′-3′): RKO: siMETTL3-Homo-929 sense, GUGCAAGAAUUCUGUGACUTT; siMETTL3-Homo-929 antisense, AGUCACAGAAUUCUUGCACTT; siYTHDC1-Homo-1883 sense, GCGUCGACCAGAAGAUUAUTT; siYTHDC1-Homo-1883 antisense, GCGUCGACCAGAAGAUUAUTT; siSTK25-Homo-518 sense, GAGACAUACAUUGCCACGATT; and siSTK25-Homo-518 antisense, UCGUGGCAAUGUAUGUCUCTT. LOVO: siMETTL3-Homo-208 sense, GCACUUGGAUCUACGGAAUTT; siMETTL3-Homo-208 antisense, AUUCCGUAGAUCCAAGUGCTT; siYTHDC1-Homo-1595 sense, GCUGGGAGGUGUCUUUAAATT; siYTHDC1-Homo-1595 antisense, UUUAAAGACACCUCCCAGCTT; siSTK25-Homo-288 sense, GCAUCGAUAACCACACAAATT; and siSTK25-Homo-288 antisense, UUUGUGUGGUUAUCGAUGCTT.

Colon cancer tissues and matched adjacent normal tissues were collected from 10 patients who underwent radical surgical resection at Fujian Medical University Union Hospital (Fuzhou, China) between December 1, 2025, and December 31, 2025. Adjacent normal tissues were obtained from areas located at least 5 cm away from the tumor margin. The cohort included 5 male and 5 female patients, with a median age of 55.5 years (range, 35–74 years). The inclusion criterion was histologically confirmed colon adenocarcinoma in patients who had not received preoperative chemotherapy or radiotherapy. Patients who had received neoadjuvant therapy or had incomplete clinical information were excluded. Written informed consent was obtained from all participants, and the study was approved by the Institutional Ethics Committee of Fujian Medical University Union Hospital (approval no. 2025KY307).

Cell apoptosis was assessed by flow cytometry using an Annexin V-FITC/PI Apoptosis Detection Kit (Beyotime Biotechnology). Briefly, cells were resuspended in binding buffer, incubated with Annexin V-FITC and PI in the dark for 30 min at room temperature. Samples were analyzed using a BD FACSCanto II flow cytometer (BD Biosciences), and data were processed using FlowJo software (version 10.8.1; BD Biosciences).

Cell pellets derived from paired tumor and adjacent normal tissues were lysed in RIPA buffer (cat. no. P0013B; Beyotime Biotechnology), and protein concentrations were determined using a BCA Protein Assay Kit. Equal amounts of protein (40 µg) were separated on 10% gels using SDS-PAGE and transferred onto PVDF membranes. Membranes were blocked with 5% non-fat milk for 2 h at room temperature and then incubated overnight at 4°C with the primary STK25 polyclonal antibody (cat. no. 25821-1-AP, 1:1,000 dilution; Proteintech Group, Inc.). The membranes were also incubated with GAPDH antibody (cat. no. 60004-1-Ig; 1:5,000 dilution; Proteintech Group, Inc.) as a loading control. After three washes with TBST (TBS containing 0.05% Tween-20; 5 min each), membranes were incubated with horseradish peroxidase-conjugated secondary antibodies (cat. no. A0208; 1:5,000 dilution; Beyotime Biotechnology) for 1 h at room temperature. Membranes were washed again as aforementioned and protein bands were detected using an enhanced chemiluminescence detection kit (cat. no. P0018; Beyotime Biotechnology), with reagents A and B mixed at a 1:1 ratio in the dark. The chemiluminescent signals were captured using a ChemiDoc XRS+ imaging system (Bio-Rad Laboratories, Inc.).

Total RNA from the cells was extracted using QIAzol reagent (Qiagen, Inc.). Subsequently, 2 µg of total RNA was reverse-transcribed into cDNA using a commercial reverse transcription kit (PrimeScript™ RT Reagent Kit; cat. no. RR047A; Takara Bio, Inc.) according to the manufacturer's instruction. PCR amplification was performed using SYBR GREEN qPCR Master Mix TB Green^® Premix Ex Taq™ II; cat. no. RR820A; Takara Bio, Inc.) on an ABI 7500 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). The thermocycling conditions were as follows: Initial denaturation at 95°C for 30 sec, followed by 40 cycles of 95°C for 5 sec and 60°C for 34 sec. GAPDH was used as an internal control. Relative gene expression levels were calculated using the 2-^ΔΔCq method (18). The primer sequences used were as follows (5′-3′): STK25 forward, GCACAGCAAGCCCTTCAAGG; STK25 reverse, CGTCCCCTTGTGAAGCTTGC; METTL3 forward, CTACGGAATCCAGAGGCAGCATTG; METTL3 reverse, GAGATGGCAAGACAGATGGACACAG; YTHDC1 forward, TGGAGAGGAGGAGGAGGAAGAGG; YTHDC1 reverse, CGGACAGCACGAACGGAAGATG; GAPDH forward, CAAGGCTGTGGGCAAGGTCATC; and GAPDH reverse, GTGTCGCTGTTGAAGTCAGAGGAG.

RIP assays were performed using a RIP Kit (cat. no. 17-700; Merck KGaA) according to the manufacturer's instructions. Briefly, magnetic beads were incubated overnight at 4°C with 5 µg of antibodies against METTL3 (Cohesion Biosciences Limited; cat. no. CQA1783) or YTHDC1 (Abcam; cat. no. ab264375) or normal IgG as a negative control (Beyotime Biotechnology, cat. no. A7016). The antibody-conjugated beads were then incubated with 500 µg of cell lysates, followed by washing five times using the RIP washing buffer provided in the kit. RNA-protein complexes were digested with proteinase K at 55°C for 30 min, and the bound RNA was subsequently eluted using the RIP elution buffer supplied with the kit, purified and analyzed by RT-qPCR to determine enrichment levels. The enrichment of target transcripts (STK25 mRNA) was assessed. The aforementioned primers used in the RT-qPCR were used for the analysis.

MeRIP was performed using the MeRIP m6A Kit (cat. no. GS-ET-001A; Cloud-Seq Biotech Ltd. Co.) following the manufacturer's protocol. Briefly, 5 µg of chemically fragmented RNA was immunoprecipitated using 5 µg of anti-m6A antibody. The precipitated RNA and input RNA (control) were then analyzed by RT-qPCR using specific primers for STK25 (5′-3′): Forward, TTCTGAAGGGCCTGGATTATCTGC and reverse, TTCTTGGGAATCAGGAACAGGACGC.

Gene sets associated with different forms of PCD were curated from published literature (19,20) and databases (FerrDb, www.zhounan.org/ferrdb; HADb, http://autophagy.lu/clustering/index.html). For each PCD type, a corresponding gene list was compiled based on these sources. The number in parentheses indicates the number of genes included in each gene set. Specifically, the curated gene sets included: Apoptosis (n=580), autophagy (n=525), lysosome-dependent cell death (n=220), necroptosis (n=105), ferroptosis (n=104), pyroptosis (n=64), immunogenic cell death (n=40), cuproptosis (n=21), entotic cell death (n=19), NETotic cell death (n=10), alkaliptosis (n=8), parthanatos (n=7), oxeiptosis (n=5) and disulfidptosis (n=4). The gene expression dataset for COAD was obtained from TCGA database (https://www.cancer.gov/ccg/research/genome-sequencing/tcga). Data preprocessing was conducted using the R package (version 4.5.0), TCGAbiolinks (version 2.30.1). Raw expression matrices were subjected to quality control (QC) with the TCGAanalyze_Preprocessing function (cor.cut=0.6) to remove low-quality samples. Technical biases were corrected using the TCGAanalyze_Normalization function and low-expression genes were filtered out with the TCGAanalyze_Filtering function (qnt.cut=0.25). Differential expression analysis between tumor and adjacent normal tissues was performed using the edge R package (version 3.42.4). Genes with P<0.05 were considered significantly differentially expressed genes (DEGs). PCD-related DEGs were identified by intersecting the DEG list with predefined PCD-related gene sets. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of upregulated and downregulated DEGs were conducted separately using the enrichGO and enrichKEGG functions from the clusterProfiler (version 4.8.1).

Univariate Cox regression analysis of PCD-related genes was performed using the survival (version 3.5–5) R package. Spearman correlation analysis between PCD-related genes and m6A regulators was conducted using the Hmisc package (version 5.1–1). Genes that were significant in the univariate analysis and correlated with m6A regulators were combined and LASSO Cox regression was applied using the glmnet package (version 4.1–8) to select the final signature genes. The risk score for each patient, designated the MCDI, was calculated using the following formula: MCDI = (0.011510365 × ATG13 exp.) + (−0.013042804 × CARS1 exp.) + (−0.020456124 × DEF8 exp.) + (0.016717586 × FBH1 exp.) + (0.020419640 × FLCN exp.) + (0.016380593 × IDUA exp.) + (−0.031405568 × ING5 exp.) + (0.048748645 × MIR210 exp.) + (0.025716205 × PIDD1 exp.) + (0.009329190 × PIP4K2B exp.) + (−0.015689712 × PLA2G15 exp.) + (0.033027922 × STK25 exp.) + (0.099700295 × TGFB2 exp.) + (0.015789901 × TRADD exp.) + (0.001302880 × TRAF2 exp.) + (0.346873020 × TRIM6 exp.) + (0.068723205 × TRIM68 exp.) + (0.006337617 × TYK2 exp.) + (−0.014371960 × USP21 exp.) + (−0.012858944 × WDR6 exp.) + (0.003176967 × ZFYVE1 exp.). Within these formulae, ‘exp.’ designates the normalized expression level of each gene in the TCGA dataset.

A co-expression network among the 21 signature genes was constructed based on Spearman correlation analysis (|R|>0.15, P<0.05) and visualized using Cytoscape software (version 3.10.1; Cytoscape Consortium).

Differential expression analysis between the MCDI-defined high- and low-risk groups in the TCGA-COAD cohort was conducted using the edge R package (version 3.42.4). Patients were stratified into high- and low-risk groups according to the optimal cutoff value of the MCDI score determined by the surv_cutpoint function in the survminer (version 0.4.9) R package. Gene expression data were first normalized using the Trimmed Mean of M-values method. Genes with a false discovery rate (FDR) <0.05 and an absolute log2 fold change (|log2FC|)>1 were considered statistically significant DEGs. The resulting DEG list was subsequently used for functional enrichment analyses.

Clinical data from the TCGA-COAD cohort and external validation cohorts (GSE39582, GSE33113, GSE38832, GSE67501, GSE78220) (21–25), including age, sex, overall tumor stage (Stages I–IV) and T, N and M stages of the TNM system, were collected. Univariate and multivariate Cox regression analyses were performed, incorporating these clinical variables alongside the MCDI score, to assess the independent prognostic value of the MCDI in predicting patient survival. Additionally, a prognostic nomogram was constructed using multivariate Cox regression combined with stepwise regression, integrating the MCDI score and clinical characteristics (age, sex, overall stage, T stage and N stage). The nomogram was visualized using the regplot (version 1.1) R package.

The association between the MCDI score and the expression levels of immunomodulatory genes was assessed by calculating Spearman correlation coefficients using the rcorr function from the R Hmisc package.

The relative proportions of 22 immune cell types in each TCGA-COAD sample were estimated using the CIBERSORT tool (https://github.com/Moonerss/CIBERSORT). The Tumor Immune Dysfunction and Exclusion (TIDE) score for each sample was calculated via the online TIDE tool (http://tide.dfci.harvard.edu/) to predict potential responses to immune checkpoint blockade (ICB) therapy. Somatic mutation data (simple nucleotide variation) for the TCGA-COAD cohort were obtained from the Genomic Data Commons portal (https://portal.gdc.cancer.gov) using the TCGAbiolinks R package. Tumor mutational burden (TMB) was calculated with the maftools package (version 2.8.05), normalizing by an exonic capture region size of 38 Mb, and the resulting values were log-transformed. Microsatellite instability (MSI) status data were retrieved from the TCGA-COAD phenotype dataset available through the UCSC Xena database (https://xenabrowser.net/datapages/).

In total, 2 colon cancer-related scRNA-seq datasets (GSE132465 and GSE205506) were used for analysis (26,27). Data processing and analysis were conducted using the Seurat R package (v4.4.0). QC was applied to remove low-quality cells. After QC, data were normalized and integrated to correct for batch effects using the reciprocal principal component analysis (RPCA) method. Cells were clustered and annotated into six major cell types based on canonical marker genes. Epithelial cells were then extracted and divided into tumor and control groups according to the sample type. These epithelial cells were subjected to unsupervised re-clustering to identify distinct subpopulations. The copy number variation (CNV) profile of each epithelial subcluster was assessed using the inferCNV (version 1.18.1) R package to evaluate malignant potential. Malignant epithelial subpopulations were identified based on CNV inference for downstream analyses. Expression patterns of the gene of interest (STK25) and a known colorectal cancer malignant epithelial marker (carcinoembryonic antigen-related cell adhesion molecule 6; CEACAM6) were examined within these malignant cells. Differential expression analysis between tumor and control groups was performed using the FindMarkers function. The MCDI score was then calculated for each single cell and its distribution was visualized across different cell types, sample groups and epithelial subpopulations.

The association between STK25 expression and clinical parameters (MCDI, T stage, N stage, M stage and age) in the TCGA-COAD cohort was assessed. Kaplan-Meier survival curves were generated using the ggsurvplot function from the survminer (version 0.4.9) R package to visualize the impact of STK25 expression levels on patient overall survival (OS).

To investigate the expression correlation between STK25 and 23 known m6A regulatory enzymes, Spearman rank correlation analysis was performed to calculate correlation coefficients (ρ) and corresponding P-values. Results were visualized using the ggplot2 (version 3.5.1) R package.

The RM2Target database (http://rm2target.canceromics.org/), which provides information on potential associations between m6A regulators and their target genes, was subsequently utilized. STK25 was input as the target gene to search for potential m6A regulatory enzymes modulating its expression. During screening, a confidence score threshold >4 was applied to ensure high reliability of the results. Based on this criterion, two potential m6A regulators were identified: YTHDC1 and METTL3.

To explore potential functional interactions among the signature genes, a protein-protein interaction (PPI) network was constructed using the STRING database (version 12.0; http://string-db.org/). The list of the 21 signature genes was uploaded to STRING, and interactions with a minimum required interaction score >0.4 (medium confidence) were retained. The resulting interaction network was then imported into Cytoscape (version 3.10.1; http://cytoscape.org/) for visualization and network analysis.

Statistical analyses were conducted using R software (version 4.3.2) and GraphPad Prism (version 9; Dotmatics). Data are presented as the mean ± standard deviation. Comparisons between two groups were performed using both paired and unpaired two-tailed Student's t-test. For comparisons among multiple groups, one-way analysis of variance (ANOVA) followed by the Tukey post hoc test was applied. When the data did not follow a normal distribution, the Wilcoxon rank-sum test was used. The χ² test was applied for categorical variables. Spearman correlation analysis was performed to evaluate associations between continuous variables. Two-sided P<0.05 was considered to indicate a statistically significant difference. Multivariate Cox proportional hazards regression analysis was performed on the 21 PCD-related genes using the coxph function in the R survival package to assess whether their expression levels were independent prognostic factors.

To characterize the landscape of PCD in COAD, the TCGA-COAD transcriptomic dataset was analyzed, focusing on 1,379 genes representing 14 distinct PCD modalities curated from the literature and databases (Fig. 1A). Differential expression analysis revealed widespread dysregulation of PCD-related genes in tumor tissues compared with adjacent normal tissues, with 438 genes significantly upregulated and 393 downregulated (Fig. 1B).

GO enrichment analysis showed that upregulated PCD genes were primarily involved in apoptosis- and autophagy-related processes, including ‘regulation of apoptotic signaling pathway’, ‘intrinsic apoptotic signaling pathway’ and ‘regulation of autophagy’. Downregulated genes were enriched in pathways associated with ‘regulation of autophagy’, ‘lysosomal organization’ and ‘myeloid leukocyte mediated immunity’, indicating that multiple PCD-related processes are transcriptionally altered in COAD rather than uniformly activated or suppressed (Fig. 1C).

KEGG analysis revealed similar patterns: ‘Lysosome’ and ‘Necroptosis’ pathways were enriched among the upregulated genes, whereas ‘autophagy - animal’, ‘Lysosome’, ‘Apoptosis’ and ‘Necroptosis’ pathways were enriched among downregulated genes (Fig. 1D). The recurrence of apoptosis-, autophagy- and lysosome-related terms in both upregulated and downregulated gene sets suggests broad remodeling of these pathways in COAD, involving simultaneous activation and inhibition of distinct gene subsets within the same cell death networks.

Overall, these results indicate that PCD-related pathways undergo complex transcriptional alterations in COAD, providing a molecular foundation for subsequent integrative analysis of m6A modification and PCD.

To develop an integrated m6A-PCD prognostic model, a systematic analysis of 1,379 genes across 14 PCD pathways was first conducted from the TCGA-COAD cohort. Univariate Cox regression identified 109 PCD-related genes significantly associated with OS (P<0.05; Fig. 2A). Concurrently, Spearman correlation analysis showed that 957 PCD genes were significantly co-expressed with multiple m6A regulators (|R|>0.2, P<0.05; Fig. 2B), indicating extensive transcriptional interplay between m6A modification and PCD in COAD. By integrating differential expression, survival association and m6A correlation analyses, 31 overlapping hub genes that were dysregulated, prognostically relevant and co-expressed with m6A regulators were identified (Fig. 2C). Subsequent LASSO Cox regression refined this panel to a 21-gene core signature, constituting the MCDI (Fig. 2D), of which 15 genes were upregulated and 6 were downregulated in COAD tumors (Fig. 2E).

To investigate functional relationships, a protein-protein interaction network was constructed, revealing a central cluster of highly connected hub proteins including FBH1, STK25, TYK2, DEF8, PLA2G15, PIP4K2B, ATG13, ING5, WDR6 and USP21 surrounded by peripheral nodes (Fig. 2F). These core genes may serve as key integrators linking m6A modification with PCD regulation, potentially underlying the heterogeneous cell-death dynamics and immune responses observed in COAD.

Using the 21 key m6A-regulated PCD genes identified above, a prognostic model for COAD, termed the MCDI, was constructed. Patients in both the internal validation cohort (TCGA-COAD) and external validation cohorts (GSE39582, GSE33113 and GSE38832) were stratified into high- and low-risk groups based on the median MCDI score. Each patient received a composite MCDI score calculated from the expression levels of the 21 genes, with higher scores indicating a higher predicted mortality risk. Kaplan-Meier analysis demonstrated that patients in the high-risk group had a significantly shorter OS time than those in the low-risk group in both the internal (TCGA-COAD) and external validation cohorts (all log-rank P<0.001), confirming the robust and reproducible prognostic performance of the model (Fig. 3A and B). Time-dependent ROC curve analysis was performed to evaluate the predictive performance of the MCDI model. The area under the ROC curve (AUC) was calculated to assess prognostic accuracy. In the training cohort, the AUC values for predicting 1- and 3-year overall survival were 0.609 and 0.572, respectively. In the testing cohort, the AUC values increased to 0.783 and 0.711, respectively. Similar predictive performance was observed in the validation cohorts, with AUC values of 0.805 and 0.699, and 0.584 and 0.635 for 1- and 3-year survival (Fig. 3C), indicating a moderate but stable prognostic performance of the MCDI model across different cohorts.

To explore the biological basis of MCDI-based risk stratification, a comprehensive transcriptomic comparison between the high- and low-risk groups in the TCGA-COAD cohort was performed. Differential expression analysis using the edge R package, with thresholds of FDR <0.05 and |log2FC|>1, identified 2,520 significantly dysregulated genes, including 802 upregulated and 1,718 downregulated genes in high-risk patients compared with low-risk patients (Fig. 4a).

Functional enrichment analysis of these DEGs was subsequently performed to elucidate their biological significance. GO and KEGG pathway analyses consistently indicated that the DEGs were highly enriched in processes related to apoptosis, autophagy, lysosomal function and immune regulation (Fig. 4B-D). Further association analyses with clinical characteristics showed that the MCDI score increased progressively with tumor advancement. Patients with advanced-stage clinicopathological features tended to exhibit higher MCDI scores. Stage IV tumors showed significantly higher MCDI scores compared with Stage I and Stage II tumors, whereas differences between other stage groups were not statistically significant (Fig. 4E). For T stage, a significant difference was observed only between T2 and T3 tumors, while other comparisons were not significant. Regarding lymph node status, patients with N1 or N2 disease exhibited significantly higher MCDI scores than those with N0, although no significant difference was observed between N1 and N2 groups. In addition, patients with distant metastasis (M1) showed significantly higher MCDI scores than those without metastasis (M0). MCDI scores were also modestly higher in patients aged >55 years and were significantly elevated in deceased patients compared with those alive (Fig. 4E). Although differences across early-stage and age subgroups were relatively modest, patients who died exhibited significantly higher MCDI scores than those who remained alive, suggesting that elevated MCDI is associated with poorer survival outcomes.

To determine whether the MCDI serves as an independent prognostic indicator beyond conventional clinical parameters, clinicopathologic data from the TCGA-COAD cohort were analyzed. Cox regression analyses confirmed the independent prognostic value of the MCDI. In univariate analysis, tumor stage, age, T, N and M stages and MCDI were all significant predictors of OS (Fig. 5A). All clinically relevant variables, including sex, were simultaneously entered into the multivariate Cox model. Although sex did not reach statistical significance in univariate analysis, it was included as a pre-specified covariate based on its clinical relevance, to control for potential confounding on the primary variables of interest (age, T stage, and MCDI).(Fig. 5B). Based on these findings, a prognostic nomogram integrating MCDI score, age and T stage was constructed to provide individualized survival probability estimates (Fig. 5C). Patients were stratified into high- and low-score groups using the surv_cutpoint function. Kaplan-Meier analysis revealed a marked survival difference between these groups (log-rank P<0.001; Fig. 5D).

Collectively, these results establish MCDI as a reliable and independent prognostic indicator for COAD, linking m6A-regulated PCD gene expression with disease progression, patient survival and immune-related molecular alterations.

Given the enrichment of immune-related pathways in the MCDI-associated gene signature, its immunologic and therapeutic relevance were systematically evaluated. First, the relationship between the MCDI score and the expression of key immunomodulators in the TCGA-COAD cohort was assessed. Spearman correlation analysis revealed that MCDI score was significantly associated with a broad array of molecules, including PDCD1LG2 (PD-L2), CD276 (B7-H3), TNFRSF4 (OX40), IL12A, IL1B, VEGFA, TGFB1 and ADORA2A, indicating potential immunoregulatory significance of the MCDI signature (Fig. 6A).

To further evaluate its predictive value for ICB therapy, the MCDI model was applied to independent immunotherapy cohorts (GSE67501, renal cell carcinoma; GSE78220, melanoma). In both datasets, patients with high MCDI scores exhibited markedly lower response rates to anti-PD-L1 treatment compared with those with low scores, suggesting that MCDI-low tumors may be more responsive to ICB therapy (Fig. 6B and C).

To elucidate the cellular composition of the tumor immune microenvironment, the CIBERSORT tool was applied to estimate the relative abundances of 22 immune cell types in TCGA-COAD samples. This analysis revealed significant differences in immune cell composition between the MCDI-defined risk groups. Compared with the low-risk group, the high-risk group exhibited significantly increased infiltration of M0 macrophages, plasma cells, CD4+ memory resting T cells and regulatory T cells (Tregs), whereas monocytes were significantly decreased. Other immune cell types did not show significant differences between the two groups (Fig. 6D and E).

Consistent with these findings, TIDE analysis was performed to quantitatively assess immune evasion potential. The MCDI-high group showed significantly higher TIDE scores, indicating a microenvironment with enhanced immune evasion and impaired antitumor immune activity (Fig. 7A-C). The relationship between MCDI risk score and two key genomic biomarkers for immunotherapy (TMB and MSI status) were further examined. No significant association or correlation was observed between the MCDI and either TMB or MSI in the TCGA-COAD cohort (Fig. 7D-F), suggesting that the prognostic and immunological features captured by the MCDI may operate independently of these established biomarkers.

Together, these results indicate that MCDI may not only predict patient prognosis but also reflect distinct immunological states and potential responsiveness to immune checkpoint inhibitors in COAD.

To identify key genes driving MCDI-associated phenotypes, multivariate Cox regression analysis was first performed on the 21-gene signature. This analysis highlighted 5 genes (MIR210, STK25, TGFB2, TRIM6 and TRIM68) as independent prognostic factors for COAD (Fig. 2G). To examine the cellular context of these genes, scRNA-seq data from colon cancer tissues (GSE132465) were analyzed. After stringent QC and batch-effect correction via RPCA integration, six major cell populations were identified (epithelial, stromal, myeloid and T, B and mast cells) and annotated based on canonical marker gene expression (Fig. 8A-C). To systematically evaluate core prognostic gene expression across cell types, a dot plot was generated for STK25, TGFB2, TRIM68 and TRIM6 (Fig. 8D), revealing distinct cell-type-specific expression patterns.

STK25 was markedly upregulated in tumor epithelial cells compared with epithelial cells derived from adjacent normal tissues, highlighting its elevated expression in the malignant epithelial compartment, whereas TRIM6, TRIM68 and TGFB2 were expressed at low or undetectable levels in this cell type (Fig. 9A and B). Notably, STK25 expression was also enriched in mast cells within tumor tissues, suggesting potential roles in both epithelial and immune compartments (Fig. 9C). Consistent with this pattern, the MCDI score projected at the single-cell level was predominantly enriched in epithelial and mast cell populations from tumor tissues. Notably, enrichment was also observed in stromal cells, paralleling the STK25 distribution and further supporting its contribution to the high-risk phenotype (Fig. 9D and E).

To confirm the reproducibility of the single-cell observations, particularly the tumor-specific enrichment of STK25 in epithelial cells, an independent colon cancer scRNA-seq dataset (GSE205506) was analyzed. Consistently, STK25 was markedly upregulated in tumor-derived epithelial cells, whereas TRIM6, TRIM68 and TGFB2 remained minimally expressed across cell types (Fig. 8A-D).

To further validate the epithelial origin of STK25 enrichment and its association with malignancy, re-clustering of epithelial cells from the integrated GSE132465 dataset was performed, identifying 7 distinct subclusters (Cluster 0–6) (Fig. 10A). Subsequent inferCNV analysis, visualized via a heatmap of chromosomal CNVs (Fig. 10B) and violin plots comparing CNV aberration scores (Fig. 10C), revealed that Clusters 0, 1, 2 and 5 exhibited significantly elevated CNV scores relative to the normal control group, indicating that these subclusters represent malignant epithelial cells. Notably, within the GSE132465 cohort, high STK25 expression (Fig. 10D), along with the established malignant epithelial marker CEACAM6 (Fig. 10E and F), was predominantly concentrated in these CNV-defined malignant subclusters (Clusters 0, 1, 2 and 5) from tumor samples, providing cross-validation and confirming tumor-specific enrichment.

Collectively, these results consistently highlighted STK25 as a tumor-enriched, prognostically relevant MCDI core gene, supporting its selection for further functional investigation.

The expression profile and clinical significance of STK25 in COAD was next evaluated. The Wilcoxon rank-sum test revealed that STK25 expression was significantly higher in patients with elevated MCDI scores. STK25 expression did not show significant differences across tumor stages (T stage) or distant metastasis status (M stage). However, higher expression levels were observed in patients with lymph node metastasis, particularly in the N2 group (Fig. 11A). Although STK25 expression showed a modest upward trend in distant metastasis (M stage), this did not reach statistical significance. Notably, younger patients exhibited slightly higher STK25 expression than older patients (Fig. 11A). Patients were divided into high- and low-expression groups based on the median STK25 level. Kaplan-Meier survival analysis demonstrated that patients with high STK25 expression had a significantly poorer OS compared with those with low expression (Fig. 11B), indicating that elevated STK25 may be associated with aggressive tumor behavior and an unfavorable prognosis.

To validate the differential expression results, qPCR and western blot analyses were performed on paired tumor and adjacent normal tissues, which consistently showed markedly higher STK25 expression in tumor tissues (Fig. 11C and D). Functional assays further explored the biological role of STK25; STK25 knockdown in LOVO cells, confirmed by qPCR (Fig. 11E), significantly increased apoptosis, as measured by flow cytometry (Fig. 11F).

Collectively, these findings indicated that STK25 may promote colon cancer progression, at least in part by suppressing apoptosis, supporting its role as a potential oncogenic driver and therapeutic target.

The upstream mechanisms driving STK25 upregulation were next investigated, focusing on potential regulation via m6A-dependent post-transcriptional modification. Spearman correlation analysis between STK25 and 23 known m6A regulators in the TCGA-COAD cohort revealed a distinct pattern: STK25 expression was positively correlated with PCIF1, ZC3H7B, YTHDF1, HNRNPD, METTL3, GNL3 and YTHDC1, while negatively correlated with METTL14, RBM27, YTHDC2, RBM15, YTHDF3, WTAP and MSI2 (Fig. 12A).

To identify the most credible upstream regulators, the RM2Target database was queried using STK25 as the target gene. METTL3 and YTHDC1 emerged as the highest-confidence regulators predicted to target STK25 (Fig. 12B), as they showed the highest confidence scores in RM2Target and were supported by binding and perturbation evidence. Both regulators also exhibited positive associations with STK25 expression, suggesting a potential regulatory mechanism involving m6A installation by METTL3 and subsequent recognition by YTHDC1. Functionally, METTL3, as an m6A ‘writer’, deposits m6A markers on RNA, while YTHDC1, a nuclear m6A ‘reader’, binds methylated transcripts to enhance stability. Consistent with this mechanism, knocking down METTL3 or YTHDC1 expression in RKO and LOVO cells reduced STK25 mRNA levels (Fig. 12C).

To determine if this regulation is direct, RIP assays were performed using antibodies against METTL3 and YTHDC1. RIP-qPCR demonstrated significant enrichment of STK25 transcripts in both METTL3 and YTHDC1 IPs compared with IgG controls (Fig. 12D), confirming direct physical interaction with STK25 mRNA. Finally, MeRIP-qPCR was conducted to verify METTL3-mediated m6A modification of STK25. Following METTL3 knockdown, m6A enrichment on STK25 mRNA was significantly reduced (Fig. 12E).

Collectively, these results demonstrated that METTL3-catalyzed m6A methylation and YTHDC1 recognition cooperatively enhanced STK25 mRNA stability, establishing an m6A-YTHDC1-STK25 axis that promotes COAD progression, with METTL3 knockdown effectively suppressing STK25 expression.