Decoding Senescence-Driven Heterogeneity in Early-Onset Colorectal Cancer for Prognostic and Therapeutic Stratification.

1
Introduction
Colorectal cancer (CRC) remains a leading cause of cancer‐related morbidity and mortality worldwide, posing a significant global health burden [1, 2, 3]. While CRC has historically been associated with aging populations, a striking epidemiological shift has emerged over the past two decades: The incidence of early‐onset colorectal cancer (EOCRC, diagnosed in individuals under 50 years) has risen by approximately 2%–4% annually worldwide, contrasting with declining rates in older cohorts [4, 5]. This demographic transition presents a pressing clinical paradox. Despite comparable overall survival outcomes between EOCRC and late‐onset CRC (LOCRC), EOCRC tumors exhibit distinct biological aggression, including higher frequencies of advanced‐stage diagnosis (58.4%–89.0% vs. 30.3%–62.5%) and lymph node metastasis (LNM, 20.96% vs. 11.90%) [6, 7]. Clinically, the absence of age‐specific therapeutic guidelines often leads to overtreatment without survival benefits, compounded by a critical lack of biomarkers to guide precision strategies [8].
Cellular senescence—a durable cell cycle arrest coupled with metabolic reprogramming and secretion of pro‐inflammatory factors (senescence‐associated secretory phenotype, SASP)—plays a dual role in cancer [9, 10, 11]. Initially, a tumor‐suppressive barrier by halting proliferation, senescent cells can paradoxically fuel tumor progression through SASP‐mediated remodeling of the tumor microenvironment (TME), fostering immune evasion, angiogenesis, and metastasis [12, 13, 14]. In CRC, the functional implications of senescence remain context‐dependent and poorly defined, particularly in EOCRC [15, 16]. Intriguingly, preliminary data suggest that compared to LOCRC, EOCRC exhibits pronounced senescence features, where senescence decouples from chronological aging and acts as an intrinsic tumor cell state that drives heightened malignancy. This challenges the view of senescence as a passive age‐related barrier, positioning it as a dynamic therapeutic target. Furthermore, EOCRC's significantly higher clinical aggressiveness (advanced stages, metastasis) underscores the prominent role of this senescence state. However, EOCRC lacks effective age‐specific guidelines and precision biomarkers, leaving the molecular heterogeneity and translational potential of these senescence pathways unexplored.
To address these critical knowledge gaps, we employed a comprehensive approach by integrating single‐cell sequencing (scRNA‐seq) and bulk RNA sequencing (RNA‐seq) data to analyze the impact of cellular senescence on EOCRC. Through this multi‐omics analysis, we identified two distinct senescence subtypes of EOCRC, characterized by significantly different biological features and prognoses. Furthermore, we constructed a senescence scoring system based on EOCRC patient data and assessed its utility as a prognostic biomarker to guide clinical treatment decisions. In summary, this study significantly contributes to our understanding of the role of senescence in EOCRC, providing new insights for the development of precision medicine strategies for EOCRC patients.

2
Material and Methods
2.1
Data Acquisition and Processing
The workflow of this study is outlined in Figure 1. A total of seven eligible bulk RNA‐seq cohorts and one single‐cell RNA‐seq (scRNA‐seq) cohort (GSE236581 [17]) were curated, all of which contained matched clinical annotations and high‐quality gene expression profiles. The seven bulk RNA‐seq cohorts comprised 2961 patients, including 473 with EOCRC and 2488 with LOCRC.
Our in‐house ICGC‐ARGO cohort (https://www.icgc‐argo.org/page/89/project‐list), generated at the Sixth Affiliated Hospital of Sun Yat‐sen University, included RNA‐seq data from 975 CRC tissue samples, among which 234 were EOCRC cases; this cohort was used as the training set. TPM (Transcripts Per Million) data for CRC patients from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) were downloaded from the TCGA portal, with 74 EOCRC cases selected as an external validation cohort.
In addition, five publicly available CRC cohorts (GSE39582, [18] GSE17538, [19] GSE39084, [20] GSE14333, [21] and GSE75315 [22]) were retrieved from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) using the “GEOquery” R package [23]. Three GEO datasets containing EOCRC survival information (GSE39582, GSE17538, and GSE39084) were merged to form the Meta‐GEO validation cohort (Table S1A,B). To correct for batch effects and ensure data comparability across platforms, the ComBat algorithm implemented in the “sva” R package was applied (Figure S1F) [24].
A predefined senescence‐related gene set (SenMayo), comprising 125 genes, was obtained from the literature by Saul et al. [25]. To assess the heterogeneity of senescence gene expression, we calculated the Euclidean distances of SenMayo gene profiles between EOCRC and LOCRC cohorts. Euclidean distance is a standard metric used to quantify differences in high‐dimensional data, measuring the straight‐line distance between samples in a multidimensional gene expression space [26]. A larger average Euclidean distance indicates greater dispersion and heterogeneity in gene expression.

2.2
Histological and Immunohistochemical Staining
Formalin‐fixed paraffin‐embedded (FFPE) tumor tissue blocks were sectioned (4‐μm). Hematoxylin and eosin (H&E) staining confirmed tumor regions. For immunohistochemical (IHC) analysis, slices were dried and processed using an automated system (BenchMark XT, Roche) for Ki67 detection (Cat No. MIBI, ZSGB‐BIO). Ki67‐positive signals appeared as brown nuclear staining. Whole‐slide images (WSIs) were acquired at 20× magnification using an Aperio scanner.

2.3
Identification of the Core Senescence‐Associated Genes
Based on transcriptomic data from the GEO, TCGA, and ICGC‐ARGO cohorts, univariate Cox regression was performed in EOCRC patients to identify genes significantly associated with disease‐free survival (DFS). Genes with p < 0.01 were considered prognostic. The intersection of these genes with the SenMayo gene set yielded a 13‐gene core senescence‐associated gene (CSG) signature for EOCRC.

2.4
Establishment of Senescence Subtypes
NMF was performed using the “NMF” package to identify distinct senescence subtypes in EOCRC patients [27], based on the expression profiles of the 13 CSGs.

2.5
Assessment of Biological Characteristics Between Clusters
To quantify the activity of specific biological pathways, we performed single‐sample gene set enrichment analysis (ssGSEA) using the “GSVA” R package [28]. Gene identifiers were converted to Entrez IDs to match the HALLMARK gene sets (h.all.v2023.1.Hs.entrez) from the Broad Institute [29], representing canonical biological pathways. Senescence‐related gene sets were retrieved from the Molecular Signatures Database (MSigDB) using keywords like “SENESCENCE,” “AGING,” and “SENESCENT” (Table S4), [30] and metastasis‐related gene sets were collected using terms such as “METASTASIS,” “EMT,” and “INVASION” (Table S5). Additionally, we assessed the activity of 10 oncogenic pathways [31].
Senescence‐associated scores, including the feature score by Lv et al. [32], Core SASP Score, and the Secretory Senescence‐Associated Secretory Phenotype (sSASP) Score by Basisty et al. [33], were used to evaluate senescence differences between EO‐Senscore subtypes. The apoptosis gene set was obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [34]. Genomic variation was analyzed using the “maftools” R package [35]. To explore tumor–immune interactions and predict immunotherapy response, immune‐related genes and signatures were analyzed using the “IOBR” R package [36].

2.6
Construction of EO‐Senscore
To quantify cellular senescence in EOCRC, we applied the LASSO algorithm using the “glmnet” R package for feature selection and model construction [37]. From a predefined senescence‐related gene set, 28 DFS‐associated genes (p < 0.05) were identified by univariate Cox regression. Ten‐fold cross‐validation was used to select the optimal penalty parameter (lambda), yielding 17 core genes. An EO‐Senscore was then calculated for each patient as a weighted sum of gene expression levels using LASSO‐derived coefficients (Figure S4B).
Patients were stratified into high‐ and low‐score groups based on an optimal cutoff determined by the “survminer” R package [38] and this classification was used in downstream analyses.

2.7
scRNA‐Seq Data Processing and Analysis
scRNA‐seq data from eight EOCRC patients receiving neoadjuvant anti‐PD‐1 therapy were retrieved from the GEO dataset GSE236581 and processed using the Seurat package (v5) [39]. Quality control (QC) was performed according to Chen et al. [17], filtering out cells with: (i) UMI counts < 600 or > 25,000, (ii) < 600 detected genes, or (iii) mitochondrial gene percentage > 5% (or > 70% in CD45− cells). Only cells passing QC were retained. Data were normalized using NormalizeData, and the top 2000 highly variable genes were identified using FindVariableFeatures. These genes were scaled with ScaleData to reduce noise. Principal component analysis (PCA) was performed for dimensionality reduction.
Patients were grouped into high‐ and low‐score categories based on EO‐Senscore pseudo‐bulk values, and comparisons were made regarding cell‐type composition and senescence‐related marker expression in tumor epithelial tissues.

2.8
Immunotherapy Response Analysis
Immune checkpoint‐related genes were obtained from Qin et al. [40] Established immunotherapy response predictors, including tumor mutation burden (TMB), Tumor Immune Dysfunction and Exclusion (TIDE) score [41], and MIRACLE score [42], were evaluated.
Additionally, three immunotherapy cohorts were used as supplementary datasets: IMvigor210 [43], comprising metastatic urothelial cancer samples treated with anti‐PD‐L1 agents, was sourced from the “IMvigor210CoreBiologies” package; GSE91061 [44], containing advanced melanoma samples treated with anti‐PD‐1 immune checkpoint blockade (ICB), was downloaded from the GEO database; and PRJEB23709, a cohort of metastatic melanoma samples treated with combination anti‐PD‐1 and anti‐CTLA‐4 therapies, was accessed from the European Molecular Biology Laboratory (EMBL, https://www.ebi.ac.uk/).

2.9
Cell Lines and Culture Conditions
Transcriptomic data of 79 colorectal adenocarcinoma cell lines were obtained from the Cancer Cell Line Encyclopedia (CCLE, https://sites.broadinstitute.org/ccle/). Human CRC cell lines (SW480, HCT116, LoVo, HCT15) were purchased from ATCC (Rockville, USA) and cultured in RPMI‐1640 or DMEM medium (Gibco, USA) supplemented with 10% FBS and 1% penicillin–streptomycin in a humidified incubator at 37°C with 5% CO2.
We investigated the sensitivity of cell lines from different groups to chemotherapeutic drugs, as predicted by the Genomics of Drug Sensitivity in Cancer (GDSC) database [45].

2.10
Cell Counting Kit‐8 (CCK8) Assays and Apoptosis Analysis
Cells (5,000/well) were plated in opaque‐walled 96‐well plates and incubated overnight in 100 μL culture medium. After 24 h, ABT‐263 (Selleck, USA) was added and incubated for 72 h (Table S6). Then, 10 μL of CCK8 solution (Abbkine, China) was added, and after 2 h, absorbance at 450 nm was measured using a Varioskan Flash reader (Thermo Scientific, USA). Drug dose–response curves were generated by nonlinear regression in GraphPad Prism 7.05. Apoptosis was assessed using the Annexin V‐APC/7‐AAD Apoptosis Kit (MultiSciences, China). Cells (3 × 105 per well) were treated with 5 μM ABT‐263 (or DMSO) for 48 h. Adherent and suspension cells were harvested, washed with cold PBS, and stained with Annexin V‐APC and 7‐AAD. Flow cytometry was performed on a CytoFLEX S flow cytometer (Beckman, USA), and data were analyzed with FlowJo V10 software. Graphing and statistical analysis were performed using GraphPad. Student's t‐test was used to compare differences between two groups. The p‐values < 0.05 were considered statistically significant.

2.11
Statistical Analysis
All data processing and statistical analyses were performed using R software (version 4.1.3). The Wilcoxon rank‐sum test was performed to assess differences between two groups, while the Kruskal‐Wallis test was used to compare variations among multiple groups. The Kaplan–Meier survival curves with log‐rank test were employed with the “survminer” R package. Correlation analyses were conducted using Spearman's rank correlation coefficient. All in vitro experiments were repeated at least three times, and the data were expressed as the means ± SEM. Statistical significance was determined by two‐tailed Student's t‐tests. Statistical significance was defined as p‐values < 0.05.

3
Results
3.1
Clinicopathological and Senescence Distinctions Between EOCRC and LOCRC
We analyzed transcriptomic profiles from a total of 2961 CRC patients after batch‐corrected integration of seven RNA‐seq datasets. Subsequent clinicopathological analyses revealed significant differences between EOCRC and LOCRC. Compared to LOCRC, EOCRC patients presented with more advanced TNM stages (59.7% vs. 51.4%, p = 0.001), a higher incidence of LNM (52.6% vs. 42.8%, p = 0.001), a higher proportion of poorly differentiated tumors (7.2% vs. 6.2%, p = 0.001), and a higher frequency of microsatellite instability high (MSI‐H) status (20.1% vs. 15.5%, p = 0.021). Additionally, EOCRC tumors were more frequently located in the left colon (64.3% vs. 57.9%, p = 0.021), while no significant differences in gender distribution were observed between the groups (Figures 2A, S1A and Table S2). However, no significant differences in overall survival (OS) and DFS were observed between EOCRC and LOCRC patients (Figure S1E).
To further explore the molecular distinctions between EOCRC and LOCRC, we performed ssGSEA on hallmark pathways. The analysis revealed that pathways associated with tumor metastasis, such as EMT, TGF‐β, and NOTCH signaling, were more enriched in EOCRC patients (Figure 2B). Gene set enrichment analysis (GSEA) further confirmed that, compared to LOCRC, EOCRC tumors exhibited significant enrichment of pathways including EMT (ES = 2.13, p < 0.001), MAPK signaling (ES = 1.26, p = 0.01), NOTCH signaling (ES = 1.55, p = 0.01), and WNT signaling (ES = 1.28, p = 0.02; Figure 2C). Additionally, other metastasis‐ and invasion‐related pathways were significantly more enriched in EOCRC than in LOCRC (Figure S1C), supporting the notion that EOCRC tumors possess a more aggressive and metastatic phenotype. Supporting these transcriptomic findings, immunohistochemical analysis of the ICGC‐ARGO cohort (n = 524) demonstrated significantly higher Ki67 expression in EOCRC tumors (Figure 2D). EOCRC showed a more dispersed Ki‐67 distribution (IQR = 20.00 vs. 11.25 for LOCRC) and a higher upper quartile (Q3 = 50 vs. 40) (Figure 2E). This suggests EOCRC tumors have a higher proportion of highly proliferative tumors. Furthermore, using a 40% threshold, EOCRC exhibited a higher prevalence of Ki‐67 high tumors (28.0% vs. 19.9% in LOCRC) (Figure S1G), reinforcing the trend toward stronger proliferative activity in EOCRC.
Interestingly, we found that the cellular senescence scores did not show a positive correlation with the patients’ actual age, but rather a trend of negative correlation (R = −0.05, p = 0.003) (Figure 2F). This minimal association suggests that the elevated senescence phenotype observed in EOCRC is largely independent of the patients’ chronological age. Notably, EOCRC patients exhibited significantly higher senescence scores (Figure S1B) and elevated expression of key biomarkers linked to cellular senescence mechanisms (Figure 2G). Senescence‐ and aging‐related pathways were also more enriched in EOCRC tissues (Figure S1D), highlighting a pronounced senescence phenotype. Additionally, Euclidean distance analysis of senescence gene expression revealed greater heterogeneity in EOCRC tumors compared to LOCRC (Figure 2H).

3.2
Identification of Senescence Subtypes in EOCRC According to CSGs
These findings suggest that cellular senescence may contribute to the increased malignancy and metastatic potential of EOCRC tumors. We hypothesized that EOCRC harbors distinct senescence‐related subtypes with different outcomes and features. To test this, we further investigated the heterogeneity of senescence in EOCRC to refine patient stratification and explore the prognostic and molecular impact of these subtypes.
Using the SenMayo gene set, we performed univariate Cox regression on DFS in the EOCRC cohort to identify CSGs. A total of 423 EOCRC samples with DFS information were included in this analysis. We applied NMF to the expression profiles of these CSGs, identifying two distinct senescence subtypes: Cluster 1 (n = 151) and Cluster 2 (n = 272) (Figure 3A). 52.6% of Cluster 1 and 58.9% of Cluster 2 patients received chemotherapy. All patients underwent surgical treatment, and no patients received immunotherapy (Table S3). Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction based on CSGs confirmed that these genes could effectively distinguish the two subtypes (Figure 3B). Survival analysis revealed that Cluster 2 was significantly associated with poorer DFS (HR = 1.95, 95% CI = 1.31–2.91, p = 7.73 × 10−4; Figure 3C) and OS (HR = 1.89, 95% CI = 1.15–3.10, p = 0.01; Figure S2A). Compared to Cluster 1, Cluster 2 exhibited more aggressive clinical features, including a higher incidence of LNM (63.1% vs. 34.0%, p = 0.001), distant metastasis (31.1% vs. 19.9%, p = 0.018), and a trend toward a higher proportion of poorly differentiated tumors (10.4% vs. 2.4%, p = 0.054), as well as more advanced TNM staging (68.0% vs. 39.3%, p = 0.001). No significant differences were observed in primary tumor location or patient gender distribution (Figures 3D, S2B and Table S3).
To characterize the senescence‐related differences between subtypes, we compared the expression of senescence biomarkers, which were consistently elevated in Cluster 2 (Figure 3E). Enrichment scores based on multiple published senescence signatures further supported this observation (Figure 3F). Additionally, senescence‐associated pathways were significantly enriched in Cluster 2 (Figure 3G). In summary, our analysis identified two distinct senescence subtypes within EOCRC, with Cluster 2 characterized by a more pronounced senescence phenotype and poorer clinical outcomes.

3.3
Biological Features of Senescence‐Associated EOCRC Subtypes
To uncover the molecular characteristics of EOCRC senescence subtypes, pathway analysis was performed (Figure 4A). Cluster 2 was enriched in pathways related to invasion and metastasis (EMT, angiogenesis), and oncogenic signaling (Wnt/β‐catenin, TGF‐β, and NOTCH). In contrast, Cluster 1 was enriched in cell cycle (MYC targets and G2M checkpoint) and apoptosis pathways (Figure S3A). GSEA confirmed metastasis enrichment in Cluster 2 (Figure 4B). Genomically, Cluster 2 showed a higher TP53 mutation frequency, and co‐occurrence and mutual exclusivity analyses of common mutations revealed potential biological interactions (Figure S2C–E). Cluster 1 exhibited significantly higher TMB (p = 0.0051) and MSI‐H cases (31.7% vs. 15.0%, p = 0.001) (Figures 3D, S2F). These findings suggest Cluster 1 patients may be more responsive to ICB therapies, providing a potential rationale for precision immunotherapy in this subgroup.

3.4
Immune Microenvironment Characterization Between EOCRC Senescence Subtypes
Cellular senescence is a cancer hallmark characterized by SASP factors modulating the TME and promoting immune evasion/immunosuppression [12, 14]. To investigate the immune landscape across EOCRC senescence subtypes, we examined their association with CMS (a classification system reflecting intrinsic tumor features) [46]. Notably, Cluster 2 was predominantly linked to CMS4 (angiogenesis, stromal infiltration) (Figure 4C), displaying a pro‐tumorigenic, immunosuppressive phenotype (immune exclusion and immunosuppression) (Figures 4D, S3D). Specifically, Cluster 2 showed increased M2 macrophages and fibroblast infiltration and elevated stromal scores (Figures 4E, S3B), suggesting a stromal‐rich, immunosuppressive TME [46]. It was also enriched in oncogenic pathways (Figure S3C).
In contrast, Cluster 1 was associated with CMS1 and CMS2 (high mutation burden, strong immune activation) (Figure 4C). It exhibited higher expression of immune checkpoint genes and increased dendritic cell (DC) infiltration (Figure 4E,F), suggesting an enhanced, immunologically active tumor microenvironment.

3.5
Development and Validation of EO‐Senscore
To overcome the challenges of clinical translation and high development costs associated with EOCRC senescence subtypes, we developed a senescence scoring model, EO‐Senscore, to guide personalized therapy for EOCRC patients. Based on senescence‐related genes from the SenMayo set, we constructed the EO‐Senscore using 17 genes selected via univariate Cox and LASSO regression in the ICGC‐ARGO training cohort (Figure S4A). The EO‐Senscore was calculated using the following formula: (0.5232 × ACVR1B) + (−0.2864 × CCL1) + (0.0024 × CCL13) + (0.0692 × CD9) + (0.0328 × CTNNB1) + (−0.2481 × CXCL10) + (−0.2638 × CXCL2) + (0.1677 × EDN1) + (0.8196 × IL2) + (0.1734 × ITPKA) + (0.0018 × JUN) + (0.2947 × PTGES) + (−0.3476 × SCAMP4) + (0.2055 × SERPINE1) + (0.1495 × SERPINE2) + (−0.1097 × TNFRSF1A) + (0.3328 × VGF). Based on this model, patients were stratified into high‐ and low‐score groups according to their EO‐Senscore (Figure S4C).
Subsequent analyses in both the training and validation cohorts revealed that patients with a high EO‐Senscore had significantly shorter DFS (Figure 5A) and exhibited stronger senescence‐related expression patterns (Figures 5B, S5A). Given the multifactorial determinants of survival in EOCRC, we performed rigorous univariate (HR = 3.61, 95% CI = 2.68–4.86, p < 0.001) and multivariate Cox regression analyses (HR = 2.40, 95% CI = 1.71–3.36, p < 0.001; Figure S4D). These results demonstrated that EO‐Senscore is an independent prognostic factor for DFS in EOCRC patients.

3.6
Molecular Characterization of EO‐Senscore Groups
To investigate biological differences, Hallmark pathway analysis showed the high‐score group was enriched in tumor invasion/metastasis (EMT, angiogenesis, and TGF‐β signaling), while the low‐score group was enriched in cell cycle regulation (MYC targets and G2M checkpoint) (Figures 5C, S5B). In both the training and validation cohorts, EO‐Senscore showed positive correlations with stromal activation and immune exclusion pathways, including EMT, TGF‐β signaling, and CAF‐mediated extracellular matrix (ECM) remodeling. Conversely, EO‐Senscore was negatively associated with immune‐activation pathways, including CD8+ T effector function, cell cycle regulation, and DNA damage repair (DDR) (Figure 5D). The Sankey diagram showed the high‐score group predominantly associated with Cluster 2 and CMS4 (active EMT) (Figures 5E, S4E). Collectively, the high‐score group corresponds to Cluster 2, sharing features of enhanced senescence, aggressive invasion, and poor prognosis.

3.7
Exploring the Immune Microenvironment of EO‐Senscore Groups
To further elucidate the immune landscape associated with EO‐Senscore, we analyzed the scRNA‐seq dataset GSE236581, which includes pre‐treatment transcriptomic data from EOCRC patients receiving neoadjuvant immunotherapy. Patients were stratified into high‐ and low‐score groups based on EO‐Senscore calculated from pseudo‐bulk expression profiles.
Analysis of tumor epithelial cells revealed distinct TME compositions. The low‐score group was enriched in immune‐activating cells, including T cells, DCs, and plasma cells, whereas the high‐score group showed elevated levels of fibroblasts, regulatory T cells (Tregs), and B cells (Figures 6A,B, S6A). Notably, the majority of patients in the low‐score group achieved favorable responses to immunotherapy, with a higher rate of complete remission (CR) (Figure 6A). Additionally, senescence‐associated markers, such as CDKN1A and CDKN2A, were significantly upregulated in the high‐score group (Figures 6A, S6B). To further characterize immune activity, ssGSEA demonstrated a negative correlation between EO‐Senscore and infiltration of activated immune cells (Figure S6C). The low‐score group also exhibited higher TMB and MIRACLE scores, alongside lower TIDE scores, suggesting stronger anti‐tumor immunity and reduced immune evasion potential (Figure 6C). These patterns were reproducible across multiple independent validation cohorts (Figure S6G).
Consistently, MSI status was more frequent in the low‐score group (Figure S6E), and the expression of immune checkpoint genes was broadly elevated (Figure S6D), further supporting enhanced immunogenicity. To further validate this prediction, in vitro experiments revealed that the low‐score CRC cell lines exhibited a significantly higher percentage of tumor cell death and elevated T‐cell IFN‐γ secretion after anti‐PD‐1 (Nivolumab) treatment, supporting that low‐senescence phenotypes are more responsive to immunotherapy (Figure S6F). Collectively, these findings demonstrate that EO‐Senscore is a potential predictive biomarker for ICB therapy, effectively stratifying low‐score tumors (immunologically active and likely responders) from high‐score tumors.

3.8
Differential Immunotherapy Responses Between EO‐Senscore Groups
The performance of EO‐Senscore in predicting the response to immune checkpoint inhibitors (ICIs) was validated in several immunotherapy cohorts. In the IMvigor210 cohort, which included 298 patients with urothelial carcinoma treated with the PD‐L1 inhibitor atezolizumab, the response rate was higher in the low‐score group compared to the high‐score group (28.4% vs. 13.0%), along with better survival outcomes (HR = 1.52, p = 4.49 × 10−3; Figure 6D). Similarly, Makarov et al.'s cohort, with 101 advanced melanoma patients treated with anti‐PD‐1 ICBs, showed a higher response rate (31.8% vs. 11.1%) and better survival outcomes in the low‐score group (HR = 1.95, p = 1.09 × 10−2; Figure 6E). Analysis of 91 metastatic melanoma patients from the PRJEB23709 cohort, who received anti‐PD‐1 and anti‐CTLA‐4 immunotherapy, revealed that the high‐score group had a lower immune response rate (44.1% vs. 71.9%) and worse overall survival (HR = 2.25, p = 3.87 × 10−2) compared to the low‐score group (Figure 6F).

3.9
In vitro Validation of EO‐Senscore in Predicting Senolytic Drug Sensitivity
To validate EO‐Senscore's predictive value in drug response, we analyzed RNA profiles of 79 CRC cell lines from CCLE and stratified them into high‐ and low‐score groups (Figure 7A). Pathway analysis confirmed the high‐score group was positively associated with metastasis‐related pathways (EMT, angiogenesis) (Figure 7B). Furthermore, GDSC prediction revealed that the low‐score cell lines were significantly more sensitive to conventional chemotherapeutic agents (including 5‐fluorouracil, oxaliplatin, and irinotecan) (Figure 7C).
In contrast, high‐score cell lines were predicted to be more responsive to the senolytic agent Navitoclax. In current clinical oncology, senolytic therapy—targeting and eliminating senescent cells based on their molecular characteristics—has emerged as a promising strategy [47, 48, 49]. Navitoclax (ABT‐263), a Bcl‐xL and Bcl‐2 inhibitor, has demonstrated potent senolytic activity by selectively inducing apoptosis in senescent tumor cells [50, 51]. Supporting this, pathway enrichment analysis revealed that the apoptosis pathway was significantly enriched in the low‐score group and negatively correlated with EO‐Senscore levels (Figure 5F,G), suggesting that tumors in the high‐score group may exhibit apoptosis resistance and could respond more favorably to senolytic therapy.
To experimentally validate this prediction, four representative cell lines—SW480 and HCT116 (high‐score), LoVo and HCT15 (low‐score)—were selected for in vitro analysis (Figure 7D). Flow cytometry analysis revealed that the total apoptotic rate was significantly higher in the high‐score group compared to the low‐score group after 5 μM Navitoclax (ABT‐263) treatment (Figure 7E,F; p < 0.01). These results were further confirmed by CCK‐8 assays, which showed enhanced sensitivity in the high‐score group, as reflected by lower half‐maximal inhibitory concentration (IC50) values compared to the low‐score group (Figure 7G,H; p < 0.01). Together, these findings support EO‐Senscore as a potential predictive biomarker for both chemotherapy response and senolytic drug sensitivity.

4
Discussion
The incidence of EOCRC has risen markedly in recent years and is frequently associated with aggressive clinical features and high tumor malignancy [5]. However, individualized therapeutic strategies tailored to the biological characteristics of EOCRC remain lacking [52, 53, 54]. In this study, we conducted a comprehensive multi‐omics analysis to investigate the molecular underpinnings of EOCRC and proposed a novel senescence‐based scoring model, EO‐Senscore. Our key findings are as follows: (i) EOCRC tumors exhibit higher levels of cellular senescence than LOCRC, independent of chronological age; (ii) EOCRC displays greater senescence‐related heterogeneity, identifying two distinct subtypes with divergent biological behaviors and outcomes; (iii) the EO‐Senscore model robustly stratifies patients by prognosis and predicts response to both immunotherapy and senolytic therapy, offering a potential tool for precision oncology in EOCRC.
Our findings challenge the traditional view of senescence as a passive, age‐driven process, revealing that tumor senescence is elevated in younger EOCRC patients and shows a trend of inverse correlation with age. EOCRC tumors with higher senescence also show increased proliferative and metastatic potential, marked by the activation of pathways like EMT, angiogenesis, and TGF‐β signaling. Despite similar clinical outcomes between EOCRC and LOCRC, EOCRC displays greater senescence‐related heterogeneity, suggesting that distinct senescence programs may influence tumor aggressiveness and treatment response. Molecular subtyping identified two EOCRC subgroups: Cluster 1 (low‐senescence) with a favorable prognosis and enriched in proliferative pathways, and Cluster 2 (high‐senescence) with a poor prognosis, enriched in pro‐metastatic pathways, high TP53 mutation frequency, and co‐activation of Wnt and Notch signaling, indicating that genomic instability and senescence drive tumor aggressiveness.
To facilitate clinical application and reduce the cost of senescence subtyping, we developed and validated the EO‐Senscore model. This score demonstrated strong and independent prognostic value across multiple cohorts. A high EO‐Senscore was significantly associated with poor outcomes, advanced TNM stage, and pro‐metastatic transcriptional programs. Notably, high‐score tumors were enriched in stromal‐ and immune‐exclusion–related features, including TGF‐β signaling, CAF‐mediated ECM remodeling, and tumor‐associated macrophages (TAM) infiltration [55, 56, 57, 58], aligning with a senescence‐driven, immunosuppressive tumor phenotype. In contrast, low‐score tumors displayed active immune characteristics, including CD8+ T cell infiltration, active cell cycle, and DNA repair pathways [59, 60], scRNA‐seq analysis further confirmed that low‐score tumors were enriched in activated immune cell subsets, such as T cells and plasma cells [60, 61], which play essential roles in antitumor immunity.
Moreover, low EO‐Senscore tumors showed elevated expression of immune checkpoint genes (e.g., CTLA‐4), higher TMB, increased MIRACLE scores, and reduced TIDE scores—features predictive of enhanced sensitivity to ICB therapy [40, 62]. Consistent with this, patients with low EO‐Senscore had better clinical responses across multiple ICI‐treated cohorts. To validate this prediction, we evaluated the sensitivity of high‐ and low‐score cell lines to anti‐PD‐1 treatment. These findings support the potential of EO‐Senscore as a predictive biomarker for ICB therapy, with low‐senescence phenotypes being more responsive to immunotherapy. Importantly, the EO‐Senscore subgroups corresponded well with the aforementioned senescence clusters (Cluster 1: Low‐senescence, Cluster 2: High‐senescence), further supporting its biological relevance and clinical utility.
Given that EO‐Senscore was negatively correlated with apoptosis‐related pathways, we hypothesized that high‐score tumors may be resistant to programmed cell death. These tumors also showed activation of the hypoxia pathway, which may promote anti‐apoptotic signaling via HIF‐1α and BCL‐2 [63]. This finding provides a mechanistic rationale for the application of senolytic therapies targeting anti‐apoptotic pathways in senescent cells. Supporting this, Navitoclax (ABT‐263)—a BCL‐2 family inhibitor—has been reported to selectively eliminate senescent cells in vivo [63, 64]. In our functional experiments, high‐score EOCRC cell lines were more susceptible to ABT‐263–induced apoptosis, suggesting that EO‐Senscore may also serve as a predictive biomarker for senolytic treatment.
To explore EO‐Senscore's specificity, we validated it in the LOCRC cohort (Figure S4F). The prognostic discrimination trend in LOCRC was considerably weaker (lower statistical significance and HR) than in EOCRC. This aligns with the design goal of our model: EO‐Senscore was specifically built using molecular features from EOCRC patients. Given the significant differences in pathogenesis, genetic mutations, and clinical characteristics between EOCRC and LOCRC [65, 66], this performance disparity underscores the importance of distinguishing between these two subtypes and highlights the unique clinical value of EO‐Senscore as a prognostic and therapeutic tool specifically for EOCRC.
Despite these promising insights, several limitations should be acknowledged. First, the predictive performance and clinical utility of EO‐Senscore require further validation in prospective clinical trials. Second, the precise molecular mechanisms underlying senescence in EOCRC remain to be fully elucidated. Finally, the potential of EO‐Senscore to guide senolytic therapy warrants additional investigation in preclinical and clinical settings. Additionally, as a retrospective study, the available cohorts may exhibit some heterogeneity and lack colorectal cancer immunotherapy datasets with annotated ICI response information, which could influence the generalizability and validation of our findings. Therefore, future studies utilizing well‐characterized CRC ICI cohorts will be necessary to further confirm the predictive potential of EO‐Senscore in immunotherapy response. In summary, this work redefines cellular senescence as a linchpin of EOCRC progression, challenging its historical perception as a passive aging marker. EO‐Senscore, the first quantitative tool linking senescence to treatment response, provides an actionable roadmap for precision oncology.

Author Contributions

Du Cai: formal analysis, funding acquisition, writing – original draft, writing – review and editing. Mingru Mai: formal analysis, visualization, writing – original draft. Rende Huang: data curation, formal analysis, writing – original draft. Haoning Qi: software, visualization. Xingzhi Feng: validation. Qianling Gao: validation. Yinmeng Zhang: investigation, visualization. Chenghang Li: data curation, software. Xiaojian Wu: funding acquisition, project administration, resources. Yize Mao: conceptualization, writing – review and editing. Zihuan Yang: project administration, writing – review and editing. Feng Gao: conceptualization, funding acquisition, project administration, resources, writing – review and editing.

Funding
This study was supported by Noncommunicable Chronic Diseases‐National Science and Technology Major Project (No. 2023ZD0501600), National Natural Science Foundation of China (No. 82272422, FG; 82573725, FG), Guangzhou Basic and Applied Basic Research Fund (No. 2024A04J9983, FG), Shenzhen Medical Research Special Project (No. D2401005, FG), Guangdong Basic and Applied Basic Research Foundation (No. 2024A1515012269, FG), Guangdong Key Research and Development Project (No. 2023B1111040003, XW), Guangdong Basic and Applied Basic Research Foundation (No. 2023B1515130008, XW), “Tianshan Talents · Leading Medical Talents in Guangdong Province” Cooperative Expert Studio (No. KSYJ2022001, XW), Guangdong Basic and Applied Basic Research Foundation (No. 2023A1515110196, DC), Postdoctoral Fellowship Program of CPSF (No. GZC20233214, DC), the program of Guangdong Provincial Clinical Research Center for Digestive Diseases (No. 2020B1111170004), National Key Clinical Discipline, Youth S&T Talent Support Programme of Guangdong Provincial Association for Science and Technology (SKXRC2025127), Guangzhou Key r&d Project (2024B01J1211).

Ethics Statement
Approval of the research protocol by an Institutional Review Board: The research was endorsed by the Ethical Committee of the Sixth Affiliated Hospital of Sun Yat‐sen University, Guangzhou, China, with ethical approval code 2025ZSLYEC‐154 on 24 March 2025.

Consent
The authors have nothing to report.

Conflicts of Interest
The authors declare that they have no competing interests. None of the authors is an Editor or Editorial Board Member of Cancer Science.

Supporting information

FIGURE S1: Differences in clinical and molecular characteristics between EOCRC and LOCRC.

FIGURE S2: Comparison of molecular characteristics between EOCRC senescence subtypes.

FIGURE S3: Immune microenvironment characterization between EOCRC senescence subtypes.

FIGURE S4: Development and validation of EO‐Senscore.

FIGURE S5: Senescence and biological characteristics of EO‐Senscore groups.

FIGURE S6: Validation of immune and pathway differences between EO‐Senscore subgroups across multiple cohorts.

TABLE S1: Summary of clinical characteristics for the study cohorts.

Table S1A: Clinical characteristics of the training and validation cohorts.

Table S1B: Clinical characteristics among the meta‐GEO cohorts.

TABLE S2: Clinical characteristics of the EOCRC and LOCRC cohorts.

TABLE S3: Clinical characteristics of the cluster 1 and cluster 2 cohorts.

TABLE S4: The summary of senescence pathways.

TABLE S5: The summary of metastasis‐related pathways.

TABLE S6: Parameters for cell proliferation.