Integrative analysis reveals the role of SUMOylation-related patterns in shaping the tumor microenvironment and predicting treatment sensitivity of colorectal cancer

Introduction

Colorectal cancer (CRC) ranks among the most malignant tumors globally, holding the third position in incidence among all cancers (1). Additionally, CRC ranks second in incidence and fourth in mortality among malignant tumors in China (2). The primary treatment strategy for CRC still revolves around surgical resection, supplemented by radiotherapy, chemotherapy, targeted therapy, and immunotherapy in a comprehensive “total management” approach (3). Despite advancements in surgical techniques and adjuvant therapies, the prognosis for CRC patients remains suboptimal, with a 5-year survival rate below 60%.

In recent years, immunotherapy has emerged as a significant adjunctive treatment for CRC patients. Immune checkpoint blockage (ICB) therapies, particularly those targeting the PD-1/PD-L1 pathway, aim to modulate and activate CD8⁺ T cells by inhibiting this signaling pathway, thereby facilitating the destruction of tumor cells. This approach has shown encouraging results in several cancers, including non-small cell lung cancer, gastric cancer, melanoma, as well as CRC (4-8). However, the response to anti-PD-1/PD-L1 immunotherapy in CRC is predominantly observed in patients with microsatellite instability-high (MSI-H)/deficient mismatch repair (dMMR), who comprise only about 5–15% of all CRC cases (9-12). For patients with microsatellite stable (MSS) or proficient mismatch repair (pMMR) CRC, who typically present with fewer neoantigens and an immune-excluded or desert tumor microenvironment (TME), the efficacy of ICB remains uncertain due to the lack of active tumor-infiltrating lymphocytes (13,14). Consequently, there is an urgent need to identify new TME classifications in CRC to better guide and predict the selection of ICB.

Post-translational modifications (PTMs) are pivotal regulatory mechanisms that modulate protein interactions, hydrophobicity, charge states, conformations, and stability, thereby influencing their functions across various biological processes (15). Identified in the mid-1990s, SUMOylation is a specific PTM involving the addition of small ubiquitin-like modifiers (SUMO). This dynamic, reversible modification is executed by a conserved enzymatic cascade comprising SUMO paralogs (SUMO1–4), the E1-activating heterodimer SAE1/SAE2, the sole E2 conjugase Ubc9, a set of SUMO E3 ligases, and SENP proteases that catalyze deconjugation (16,17). It plays a critical role in cellular growth, migration, and oncogenesis (18). Moreover, SUMOylation regulates the production and activity of interferons (IFNs) by modulating transcription factors either through inhibition or stimulation, leading to attenuated immune responses (19). It also influences IFN generation by regulating the cyclic GMP-AMP synthase (cGAS) and the stimulator of interferon genes (STING) pathway, which may be linked to the activation of CD8⁺ T cells (20,21). Given its contributions to tumorigenesis, disease progression, and immune regulation, SUMOylation may serve as a promising target for precision intervention. We present this article in accordance with the TRIPOD and MDAR reporting checklists (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2194/rc).

Methods

Data acquisition and integration

We retrospectively obtained 1,226 CRC gene-expression profiles and clinical annotations from Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA). Cases lacking survival data were excluded. Four cohorts were analyzed: GSE103479, GSE39582, TCGA-Colon Adenocarcinoma (COAD), and TCGA-Rectal Adenocarcinoma (READ). For TCGA, RNA sequencing (RNA-seq) expression [fragments per kilobase of transcript per million mapped reads (FPKM)] files were downloaded from the GDC via the TCGABiolinks R package and converted to transcripts per million (TPM). Cross-dataset batch effects were corrected using ComBat (sva). Somatic mutation data for TCGA-COAD/READ were also retrieved from the Genomic Data Commons (GDC).

Unsupervised consensus clustering for SUMOylation-related genes (SRGs)

Using MSigDB [gene set enrichment analysis (GSEA)], we assembled 189 SRGs and extracted their expression from a meta-cohort. To map SUMOylation patterns, patients were grouped via unsupervised clustering based on SRG expression. The final partition and stability metrics were established with ConsensusClusterPlus after 1,000 resampling runs.

Functional annotation and gene set variation analysis (GSVA)

We applied GSVA to profile pathway-level heterogeneity across copper-toxicity patterns. GSVA is a non-parametric, unsupervised approach that infers relative activity of pathways and biological processes from expression data. Hallmark gene sets from MSigDB (h.all.v7.4.symbols.gmt) were used.

Construction of SUMOylation-related risk model

In this study, samples from the meta-cohort were randomly divided into a control group and a training group. Within the training group, univariate Cox regression analysis was employed to identify differentially expressed genes (DEGs) associated with prognosis. Hazard ratios (HRs) for each DEG were calculated, with a P value of less than 0.05 considered prognostically significant. Based on these prognostic DEGs, least absolute shrinkage and selection operator (LASSO) Cox regression analysis was conducted using the glmnet R package. The optimal lambda value was determined using 10-fold cross-validation, selecting the minimal model (minimizing observed error) for developing a SUMOylation-related risk model. The risk score for all CRC samples was calculated using the following formula: riskscore=∑li(Coefi∗Expi), where ‘Coefi’ represents the regression coefficients of genes determined by LASSO Cox regression analysis, and ‘Expi’ denotes the expression values of genes in the SUMOylation-related risk model.

Analysis of immune cell infiltration and prediction for immune response

Algorithms, including single-sample GSEA (ssGSEA) (22), xCELL (23), EPIC (24), and the MCP-counter (25), were employed to quantify the relative abundance of immune cell infiltration in the TME and assess the immune function in CRC patients. The ESTIMATE algorithm was utilized to estimate the levels of immune and stromal cells in CRC by calculating immune and stromal scores, aiding in the prediction of immune cell and stromal cell infiltration (26). The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm, encompassing the two principal mechanisms of tumor immune escape-T cell dysfunction and exclusion, was applied to evaluate the TME and predict responses to ICB (27). A higher TIDE score indicates a lower likelihood of response to ICI therapy due to the tumor’s tendency to induce immune escape. The immune phenotype score (IPS) serves as a robust predictor of the efficacy of anti-CTLA-4 and anti-PD-1 regimens, quantifying determinants of tumor immunogenicity and characterizing the immune landscape and cancer immunogenomics within tumors (22). A higher IPS suggests potential sensitivity to ICB.

Chemotherapeutic response prediction

We estimated 5-fluorouracil (5-FU) sensitivity for each CRC sample using the GDSC resource and the R package pRRophetic, which applies ridge regression to infer half-maximal inhibitory concentration (IC₅₀) values. Prediction performance was evaluated by 10-fold cross-validation on the GDSC training set.

For the external pharmacogenomic context, expression and mutation profiles were taken from the Broad Institute’s Cancer Cell Line Encyclopedia (CCLE), and drug response data from Cancer Therapeutics Response Portal (CTRP) (481 compounds tested across 835 cell lines) (28,29). Lower area under the curve (AUC) values indicate greater sensitivity. Following Yang et al., hematopoietic/lymphoid lineages were excluded, and compounds with >20% missing AUCs were removed; remaining missing values were imputed by k-nearest neighbors (kNN), and CCLE molecular features were used for downstream analyses of CTRP data (30).

Genomic and clinical data collection and integration for the ICB cohort

In our study, four cohorts with expression and clinical data from immunotherapy treatments were included. Patients with metastatic melanoma were treated with pembrolizumab or nivolumab (31), patients with non-small cell lung cancer received nivolumab or pembrolizumab (32), and patients with urothelial cancer were treated with anti-PD-1/PD-L1 therapies (33), specifically with the anti-PD-L1 antibody atezolizumab in the IMvigor 210 cohort (34). Transcription expression profiles were compiled and converted to TPM format for further analysis.

Single-cell RNA-seq (scRNA-seq) analysis

We started by getting four CRC tissue samples from Ruijin Hospital. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Ruijin Hospital, Shanghai Jiao Tong University School of Medicine [No. 2022 (138)]. Informed consent was taken from all the patients. We cleaned these samples to remove fat and blood vessels, then broke them down into individual cells using both enzymes and mechanical force. For the sequencing step, we used the 10X Chromium system and ran the samples on a NovaSeq 6000 machine. We converted the raw data into FASTQ format and checked the quality before mapping everything to the human genome. We used Seurat software to filter the data and group the cells. To handle differences between batches, we used the Harmony tool. Finally, we identified specific marker genes and looked at the activity of metabolic pathways.

Cell culture and quantitative real-time polymerase chain reaction (qRT-PCR)

The CRC lines RKO (SCSP-5236), HCT116 (SCSP-5076), HT-29 (SCSP-5032), SW480 (SCSP-5033), SW620 (SCSP-5234), and the normal colonic line NCM460 (GNHu74) were obtained from the National Collection of Authenticated Cell Cultures (NCACC) (Shanghai, China). Cells were maintained in RPMI-1640 (C11875500BT) with 10% fetal bovine serum and 1% penicillin-streptomycin (15140122) at 37 °C in a humidified 5% CO₂ atmosphere. Total RNA was isolated using the FastPure^® Cell/Tissue Total RNA Isolation Kit V2 (RC112-01) and quantified on a NanoDrop 2000. Complementary DNA (cDNA) synthesis employed HiScript^® RT SuperMix for quantitative polymerase chain reaction (qPCR) with genomic DNA (gDNA) wiper (R323-01), and qRT-PCR was performed with ChamQ Universal SYBR qPCR Master Mix (Q711-02). Relative messenger RNA (mRNA) levels of SUMOylation-based scoring system (SUMOscore) genes were calculated by the 2^−ΔΔCt method. Primer sequences are listed in Table S1.

Immunohistochemistry

Fixed, paraffin-embedded CRC tissues were deparaffinized, rehydrated, and subjected to antigen retrieval. Sections were incubated with primary antibodies overnight at 4 °C (Cell Signaling Technology, Danvers, MA, USA; 2808, O15392), followed by a biotinylated secondary antibody for 30 min at 37 °C. Signal was developed with 3,3'-diaminobenzidine (DAB), counterstained with hematoxylin, and imaged by light microscopy.

Statistical analysis

All analyses were performed in R (v4.3.2). Two-sided P values <0.05 were deemed significant. HRs for SRGs and SUMOylation subtypes were estimated using univariable and multivariable Cox models; the latter was fitted in patients with complete clinical data to identify independent prognostic factors. Results from Cox analyses were visualized with the forestplot package. Correlations were assessed by Spearman and distance correlation. Between-group differences were tested with the Wilcoxon rank-sum test. Somatic mutation landscapes in TCGA-CRC were displayed using the maftools waterfall function.

Results

SUMOylation subtypes based on SRGs

We merged four datasets with survival information and clinical annotations (GSE103479, GSE39582, TCGA-COAD, and TCGA-READ) into a single meta-cohort. Using the ConsensusClusterPlus R package, we stratified samples exhibiting different SUMOylation patterns based on the expression levels of SRGs. This unsupervised clustering approach identified two distinct SUMOylation patterns, encompassing 793 cases in cluster C1 and 433 in cluster C2 (Figure 1A-1D). Prognostic analysis of the two clusters revealed that C2 exhibited a significant survival advantage, whereas C1 had the lowest survival probability within the meta-cohort (Figure 1E). GSVA enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets was conducted to explore the molecular mechanisms differentiating the two SUMOylation patterns (Figure 1F). Notably, immune signaling pathways were significantly enriched in C2, including B cell receptor signaling, Toll-like receptor signaling, cytokine-cytokine receptor interaction, natural killer (NK) cell-mediated cytotoxicity, and chemokine signaling pathways, suggesting that the superior prognosis in C2 may be related to higher immune cell infiltration and activation. Box plots constructed via ssGSEA clarified and compared the infiltration of 23 immune cell subgroups in each cluster (Figure 1G). As hypothesized, higher infiltration of CD8⁺ T cells, B cells, NK cells, and dendritic cells was observed in C2. However, cells associated with immune suppression and escape, such as myeloid-derived suppressor cells (MDSCs) and regulatory T cells (Tregs), were also more prevalent in C2. To further assess the differences in the TME between the two clusters, we evaluated their TME scores, a novel concept from previous studies representing the characteristics of the tumor immune microenvironment as TME score A – TME score B (34). Contrary to expectations, we found a relatively higher TME score B and significantly lower TME score in C2, indicating higher stromal cell infiltration, characteristic of an immune-excluded phenotype, which contradicts the favorable prognosis observed in this cluster (Figure 1H). Thus, we contend that stratification of CRC patients based solely on SRGs expression levels is insufficient and requires further subtyping.

Figure 1 SUMOylation pattern in CRC patients. (A) Unsupervised consensus clustering for SRGs in meta-cohort and the consensus matrices for k=2. (B) Consensus values range from 0 to 1. (C) Area under CDF curve. (D) The transcriptome profiles of two SUMOylation patterns were analyzed by principal component analysis, revealing a striking difference in transcriptome profiles between different patterns. (E) Kaplan-Meier curves for the two SUMOylation patterns based on 1,226 CRC patients from meta-cohort (log-rank test). The cluster C2 showed a significant better prognosis than C1. (F) Barplot depicting the GSVA score of representative KEGG pathways curated from MSigDB in two SUMOylation patterns. (G) The fraction of TME cell infiltration of two SUMOylation patterns using ssGSEA algorithm. The top end, median line, and bottom end of the box represent the 25%, 50%, and 75% value, respectively. the black dots showed outliers. (H) The fraction of different signatures (immune-relevant signature, mismatch-relevant signature, and stromal-relevant signature) and TME score. The line in the box represented the median value. The bottom and top of the boxes are the 25^th and 75^th percentiles (interquartile range). The whiskers encompassed 1.5 times the interquartile range. The asterisks illustrated the statistical P value [*, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001; ns, not significant (P>0.05)]. aDC, activated dendritic cell; APM, antigen processing machinery; CDF, cumulative distribution function; CRC, colorectal cancer; EMT, epithelial-mesenchymal transition; GSVA, gene set variation analysis; iDC, immature dendritic cell; KEGG, Kyoto Encyclopedia of Genes and Genomes; MDSC, myeloid-derived suppressor cell; NER, nucleotide excision repair; NK, natural killer; OS, overall survival; PC, principal component; pDC, plasmacytoid dendritic cell; ssGSEA, single-sample gene set enrichment analysis; SUMO, small ubiquitin-like modifiers; TBR, transforming growth factor-β response; Th, T helper; TME, tumor microenvironment; Treg, regulatory T cell.

Construction of SUMOylation-related risk model

By comparing the expression of SRGs in clusters C1 and C2, we identified 97 DEGs with false discovery rate (FDR) <0.05 and |log fold change (logFC)| >1 in the meta-cohort, comprising 77 upregulated genes and 20 downregulated genes (Figure 2A). A SUMOylation-associated risk model was constructed through LASSO Cox regression analysis, enabling comprehensive assessment of each CRC sample (Figure 2B). We identified five DEGs associated with CRC prognosis, forming the basis of the risk model defined as: SUMOscore = −0.1651 × BIRC5 + 0.0193 × NCOA2 − 0.0543 × NR3C2 − 0.0154 × TOP1 + 0.1800 × UHRF2.

Figure 2 Construction of SUMOscore risk model the correlation between SUMOscore and CRC clinical characteristics. (A) Volcano plot of differential expression genes between cluster C1 and cluster C2. (B) Construction of SUMOylation-related risk model by LASSO Cox regression analysis. (C-F) Kaplan-Meier curves for the high and low SUMOscore groups from meta-cohort (C), TCGA-CRC-cohort (D), GSE39582-cohort (E), and GSE103479-cohort (F) (log-rank test). (G-I) The ROC curves for OS prediction by SUMOscore, gender, age, T stage, N stage, and MSI status at 1- (G), 3- (H), and 5-year (I). (J-L) Differences in clinical characteristics such as MSI status (J), T stage (K), and N stage (L) of CRC cohort between high and low SUMOscore groups. (M) The difference of mRNA expression level of 5 key SRGs between normal and tumor CRC samples (TCGA-cohort). The asterisks represented the statistical P value [*, P<0.05; ***, P<0.001; ****, P<0.0001; ns, not significant (P>0.05)]. AUC, area under the curve; CRC, colorectal cancer; dMMR, deficient mismatch repair; LASSO, least absolute shrinkage and selection operator; mRNA, messenger RNA; MSI, microsatellite instability; MSS, microsatellite stable; N, node; OS, overall survival; pMMR, proficient mismatch repair; ROC, receiver operating characteristic; SRG, SUMOylation-related gene; SUMO, small ubiquitin-like modifiers; SUMOscore, SUMOylation-based scoring system; T, tumor; TCGA, The Cancer Genome Atlas.

The association between SUMOscore and clinical features

We subsequently assessed the prognostic capability of the SUMOscore by stratifying patients into high or low score groups using an optimal cut-off value. Patients with lower SUMOscore were significantly associated with better prognoses (Figure 2C). As a further validation, we explored the relationship between SUMOscore and patient outcomes across three CRC cohorts (GSE103479, GSE39582, and TCGA-CRC), finding a consistent association between lower SUMOscore and favorable survival outcomes (Figure 2D-2F). Additionally, time-dependent receiver operating characteristic (ROC) curves were constructed to assess the accuracy of SUMOscore in predicting overall survival (OS), with 1-, 3-, and 5-year AUCs of 0.621, 0.648, and 0.627 respectively, indicating the predictive relevance of SUMOscore for OS (Figure 2G-2I). Moreover, SUMOscore demonstrated superior prognostic performance compared to conventional clinical prognostic markers such as tumor stage (T stage) and node stage (N stage). Further analysis revealed a significant correlation between higher SUMOscore and advanced CRC stages (T3–T4, N1–N3) and a notably lower SUMOscore in MSI-H/dMMR patients (Figure 2J-2L). Lastly, the expression levels of the five key model genes were explored in both normal and tumor tissues of CRC, showing higher expression of BIRC5, TOP1, and UHRF1 in tumor tissues, while NCOA2 and NR3C2 were less expressed, consistent with the correlative efficiency of these model genes (Figure 2M). Patient age, gender, T stage, N stage, MSI status, and SUMOscore were included in univariate and multivariate Cox regression models to assess prognostic factors for survival outcomes (Figure 3A,3B). The analysis revealed that age [HR =1.933; 95% confidence interval (CI): 1.505–2.484; P<0.001], N stage (HR =1.942; 95% CI: 1.541–2.447; P<0.001), and SUMOscore (HR =8.724; 95% CI: 3.412–22.307; P<0.001) are reliable and independent risk factors for evaluating patient survival outcomes. Based on these independent prognostic markers, we developed a nomogram to assess 1-, 3-, and 5-year survival rates for CRC patients (Figure 3C). Time-dependent calibration curves demonstrated that the nomogram predicts CRC prognosis with high accuracy (Figure 3D-3F). Decision curve analysis indicated that the nomogram possesses substantial predictive capability (Figure 3G,3H). In summary, we developed a SUMOscore model that effectively predicts the prognosis of CRC patients, linking higher SUMOscore with poorer outcomes, advanced disease stages, and MSS status.

Figure 3 A nomogram for survival prediction in CRC. (A) Univariate Cox regression analysis to identify risk factors in CRC survival. (B) Multivariate Cox regression analysis to identify independent prognostic indicators in CRC. (C) A nomogram composed of independent prognostic indicators (SUMOscore, N stage, age). (D-F) The calibration curve of nomogram at 1- (D), 3- (E), and 5-year (F), respectively. (G,H) Decision curve analysis of nomogram. The asterisks represented the statistical P value (***, P<0.001). CI, confidence interval; CRC, colorectal cancer; MSI-H, microsatellite instability-high; N, node; OS, overall survival; SUMO, small ubiquitin-like modifiers; SUMOscore, SUMOylation-based scoring system; T, tumor.

Differences in biological features between high and low SUMOscore groups

GSVA enrichment analysis on groups with high and low SUMOscore identified significant enrichment of WNT, MAPK, and ERBB signaling pathways in the high SUMOscore group (Figure 4A). These pathways are significantly associated with epithelial-mesenchymal transition (EMT), crucial for stromal infiltration in the TME. To confirm these findings, we used the ESTIMATE algorithm to analyze groups with different SUMOscore. Consistently, we observed that stromal scores, indicating the abundance of stromal cells within tumor samples, were significantly higher in the high SUMOscore group (Figure 4B,4C). Elevated stromal scores facilitate tumor growth and metastasis and also affect the efficacy of immunotherapy. Additionally, we applied xCell, MCPcounter, ssGSEA, and EPIC algorithms to elucidate the immune landscape of the TME in high and low SUMOscore groups. As shown in Figure 4D, a significant stromal component infiltration (such as endothelial cells and tumor-associated fibroblasts) was evident in the high SUMOscore group, with stromal component infiltration positively correlating with SUMOscore. Moreover, SUMOscore was significantly associated with the infiltration of immunosuppressive cells such as tumor-associated macrophages, MDSCs, and Tregs. To further characterize the features associated with SUMOylation, we examined the correlations between characteristics and SUMOscore. Heatmaps of the correlation matrix revealed positive correlations between SUMOscore and EMT and pan-fibroblast TGF-β response signatures, while negative correlations were noted with DNA repair and mismatch repair signatures (Figure 4E). Interestingly, ESTIMATE immune scores did not differ significantly between high and low SUMOscore groups, and there were no significant differences in classical cytotoxic cells like NK and CD8⁺ T cells (Figure 4B,4D). Differential analysis of immune activation, antigen presentation, and immune checkpoint-related markers between the high and low SUMOscore groups did not reveal significant differences (Figure 4F). Overall, our results suggest that the SUMOscore primarily reflects the infiltration of stromal components and immunosuppressive cells in the CRC TME, thus potentially associating high SUMOscore with poor patient prognosis.

Figure 4 Biological and TME infiltration characteristics of high and low SUMOscore groups. (A) Barplot depicting the GSVA score of representative KEGG pathways curated from MSigDB in high and low SUMOscore groups. (B) The comparison of immune score between high and low SUMOscore groups. (C) The comparison of stromal score between high and low SUMOscore groups. (D) The immune landscape in meta-cohort of high and low SUMOscore groups as depicted by ssGSEA, EPIC, MCPcounter, and Xcell. (E) Correlations between SUMOscore and the known biological gene signatures in meta-cohort using Spearman analysis. Negative correlation was marked with blue and positive correlation with orange. (F) The comparison of expression level of antigen presentation, immune-activation and immune-checkpoint gene sets between high and low SUMOscore groups. The asterisks represented the statistical P value (*, P<0.05). FC, fold change; GSVA, gene set variation analysis; ssGSEA, single-sample gene set enrichment analysis; SUMO, small ubiquitin-like modifiers; SUMOscore, SUMOylation-based scoring system; TME, tumor microenvironment.

Predictive role of SUMOscore in immune response and drug sensitivity

Current studies confirm that EMT is associated with resistance to chemotherapy and immunotherapy. Therefore, we predicted the drug responses of patients in high and low SUMOscore groups to guide clinical practice. Since chemotherapy, particularly 5-FU, is a common treatment for CRC, we evaluated its efficacy. Using a predictive model for this drug, we estimated the IC₅₀ for each sample in the meta-cohort. There was a statistically significant difference in 5-FU sensitivity between the low and high SUMOscore groups (P<0.001, Figure 5A). Correlation analysis showed that the IC₅₀ values of 5-FU were significantly positively correlated with the SUMOscore (Figure 5B). We built a drug-response prediction framework using the CTRP resource, paired with expression profiles across hundreds of cancer cell lines. Compounds with >20% missing AUC values and hematopoietic/lymphoid lineages were excluded; remaining not available values (NAs) were imputed by kNN, yielding 680 cell lines and 354 compounds for analysis. For each agent-sample pair, AUC was estimated from expression data via ridge regression (propredict), and associations with SUMOscore were evaluated using Spearman correlation. Consistent with prior findings, lower SUMOscore corresponded to greater 5-FU sensitivity (lower AUC) (Figure 5C,5D).

Figure 5 Prediction of drug sensitivity according to SUMOscore. (A) The differences of response to 5-FU between the high and low SUMOscore groups. (B) The correlation between SUMOscore of patients and the estimated IC₅₀ value of 5-FU. (C) Results of Spearman correlation analysis between SUMOscore and CTRP-derived compounds. (D) CTRP derived compounds in the high and low SUMOscore groups of AUC values. (E,F) The landscape of genetic mutation in the low (E) and high (F) SUMOscore group. (G) The difference of TMB between the high and low SUMOscore groups. (H) The difference of IPS between the high and low SUMOscore groups. (I) The relative distribution of TIDE was compared between high and low SUMOscore groups. (J) Prediction of responses of CRC samples between high and low SUMOscore groups to immunotherapies. The asterisks represented the statistical P value (*, P<0.05; ***, P<0.001). 5-FU, 5-fluorouracil; AUC, area under the curve; CRC, colorectal cancer; CTRP, Cancer Therapeutics Response Portal; IC₅₀, half-maximal inhibitory concentration; IPS, immune phenotype score; SUMO, small ubiquitin-like modifiers; SUMOscore, SUMOylation-based scoring system; TIDE, Tumor Immune Dysfunction and Exclusion; TMB, tumor mutation burden.

Undoubtedly, anti-CTLA-4/PD-1 therapies have achieved significant breakthroughs in oncological treatments. However, in CRC, only a minority of patients with MSI-H/dMMR benefit from immunotherapy. Our previous findings indicated that patients with MSI-H/dMMR exhibit lower SUMOscore, suggesting that the SUMOscore may predict the efficacy of immunotherapy. Tumor mutation burden (TMB) is another predictor of immunotherapy efficacy, where patients with high TMB appear responsive and experience durable clinical outcomes. Using maftools, we compared somatic mutation landscapes between high- and low-SUMOscore groups in TCGA-CRC. Tumors with high SUMOscore showed a reduced mutation burden relative to low-SUMOscore tumors (Figure 5E-5G). We noted a significant increase of IPS in the low SUMOscore group (Figure 5H). Furthermore, TIDE, a newly recognized predictor of immune response, is widely applied and highly recommended alongside well-known markers such as TMB, PD-L1, and MSI. According to our analysis, TIDE values were significantly lower in the low SUMOscore group within the meta-cohort, indicating sensitivity to anti-CTLA-4 immunotherapy (Figure 5I,5J). Given the SUMOscore’s potential to predict sensitivity to immunotherapy, we subsequently investigated whether it could forecast responses to ICB therapy across four immunotherapy cohorts. In these ICB cohorts, patients were similarly stratified into high and low groups based on the same five model genes. Consistent with findings in the CRC cohorts, high SUMOscore in the ICB cohorts was associated with increased stromal components and infiltration of immunosuppressive cells (Figure 6A,6B). Moreover, the SUMOscore showed no significant correlation with classic CD8⁺ T cell infiltration, immune activation, immune checkpoints, or antigen presentation (Figure 6C). We also analyzed TMB and neoantigen load in high and low SUMOscore groups, discovering that the low SUMOscore group exhibited higher TMB and neoantigen content (Figure 6D,6E). Consistent with the CRC cohort, the ICB cohort also demonstrated lower TIDE values in the low SUMOscore group, suggesting potential sensitivity to anti-CTLA-4 treatment (Figure 6F,6G). Exploring immunotherapy sensitivity in the ICB cohort, we found that patients with low SUMOscore significantly benefited from ICB therapy, with a higher response rate (29% vs. 19%, as shown in Figure 6H). Figure 6I further indicated that patients achieving complete response (CR)/partial response (PR) had lower SUMOscore. Taken together, these findings indirectly demonstrate the crucial role of tumor SUMOylation patterns in responses to immunotherapy and chemotherapy.

Figure 6 Potential immunotherapy in high and low SUMOscore groups. (A) The immune landscape in ICB-cohort of high and low SUMOscore groups as depicted by ssGSEA, EPIC, MCPcounter, and Xcell. (B) The fraction of different signatures (immune-relevant signature, mismatch-relevant signature, and stromal-relevant signature) and TME score in ICB-cohort. The line in the box represented the median value. The bottom and top of the boxes are the 25^th and 75^th percentiles (interquartile range). The whiskers encompassed 1.5 times the interquartile range. (C) The comparison of expression level of antigen presentation, immune-activation and immune-checkpoint gene sets between high and low SUMOscore groups in ICB-cohort. (D) The difference of TMB between the high and low SUMOscore groups in ICB-cohort. (E) The difference of neoantigen level between the high and low SUMOscore groups in ICB-cohort. (F) The relative distribution of TIDE in ICB-cohort was compared between high and low SUMOscore groups. (G) Prediction of responses of the ICB-cohort between high and low SUMOscore groups to immunotherapies. (H,I) The fraction of patients with immunotherapy response (ICB-cohort) in low and high SUMOscore group (H). The SUMOscore of CR/PR and PD/SD patients in ICB-cohort (I). The asterisks represented the statistical P value (*, P<0.05; **, P<0.01; ***, P<0.001; ns, P>0.05). CR, complete response; CRC, colorectal cancer; EMT, epithelial-mesenchymal transition; FC, fold change; ICB, immune checkpoint blockade; NER, nucleotide excision repair; PD, progressive disease; PR, partial response; SD, stable disease; ssGSEA, single-sample gene set enrichment analysis; SUMO, small ubiquitin-like modifiers; SUMOscore, SUMOylation-based scoring system; TBR, transforming growth factor-β response; TIDE, Tumor Immune Dysfunction and Exclusion; TMB, tumor mutation burden; TME, tumor microenvironment.

BIRC5 serves as an immune-associated prognostic biomarker

Among several key candidate genes, we ultimately selected the anti-apoptotic protein BIRC5 for subsequent investigations. Survival analyses across four independent cohort consistently revealed that patients with high BIRC5 expression exhibited significantly poorer OS compared to those with low expression (Figure 7A-7D), suggesting its prognostic relevance. To investigate the immunological context of BIRC5 expression, ssGSEA analysis was performed. The results indicated a markedly higher infiltration of CD8⁺ T cells and CD4⁺ T cells in the BIRC5-low subgroup, highlighting a more active antitumor immune microenvironment (Figure 7E). Furthermore, GSVA pathway enrichment analysis demonstrated that BIRC5-low tumors were positively associated with cell death-related pathways, including apoptosis, ferroptosis, and p53 signaling. In contrast, BIRC5-high tumors showed significant enrichment of pathways involved in tumor progression and stromal remodeling, such as Wnt signaling, PI3K-AKT signaling, TGF-β signaling, and the RAS pathway (Figure 7F). The ESTIMATE algorithm further supported these findings, showing that the BIRC5-high group exhibited a significantly higher stromal score (Figure 7G). Drug sensitivity analysis revealed that tumors with high BIRC5 expression displayed higher resistance to 5-FU and paclitaxel (Figure 7H,7I). Furthermore, compound screening from CTRP and PRISM databases also suggested that BIRC5-low tumors may be more sensitive to 5-FU, supporting its potential utility as a biomarker for treatment stratification (Figure 7J-7M). To validate the transcriptional profiles, qPCR was conducted in normal intestinal epithelial cells and various CRC cell lines. The results confirmed upregulated expression of BIRC5, UHRF2, and TOP1, and downregulated expression of NCOA2 and NR3C2 in tumor cells (Figure 7N). Additionally, immunohistochemistry revealed that BIRC5 protein levels were significantly elevated in CRC tissues compared to normal controls, corroborating the transcriptional findings (Figure 7O). In summary, BIRC5 is a robust prognostic biomarker in CRC that is associated with distinct immune infiltration patterns, oncogenic signaling pathways, and differential drug sensitivities. These findings highlight its potential as both a predictive and therapeutic stratification marker in clinical practice.

Figure 7 Comprehensive analysis of BIRC5 expression in CRC and its association with prognosis, immune infiltration, signaling pathways, drug sensitivity, and gene expression. (A-D) Kaplan-Meier survival curves comparing OS between high and low BIRC5 expression groups in meta-cohort (A), TCGA-cohort (B), GSE39582-cohort (C), and GSE103479-cohort (D); statistical significance was evaluated using the log-rank test. (E) The fraction of TME cell infiltration of high and low BIRC5 expression groups using ssGSEA algorithm. The top end, median line, and bottom end of the box represent the 25%, 50%, and 75% value, respectively. The black dots showed outliers. (F) Barplot depicting the GSVA score of representative KEGG pathways curated from MSigDB in the BIRC5 high vs. low expression groups. (G) The immune score and stromal score of BIRC5 high and low expression groups using ESTIMATE algorithm. (H,I) The differences of response to 5-FU (H) and paclitaxel (I) between the BIRC5 high and low expression groups. (J,K) The results of Spearman’s correlation analysis of 12 CTRP-derived compounds (J) and the AUC value between the BIRC5 high and low expression groups (K). (L,M) The results of Spearman’s correlation analysis of 12 PRISM-derived compounds (L) and the AUC value between the BIRC5 high and low expression groups (M). (N) Relative mRNA expression levels of BIRC5 and selected genes (UHRF2, NCOA2, NR3C2, and TOP1) across CRC cell lines and normal colon epithelial cells (NCM460). (O) Expression of BIRC5 in tumor and adjacent normal tissues as examined by immunohistochemistry. The asterisks represented the statistical P value [*, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001; ns, not significant (P>0.05)]. 5-FU, 5-fluorouracil; AUC, area under the curve; CRC, colorectal cancer; CTRP, Cancer Therapeutics Response Portal; GSVA, gene set variation analysis; IC₅₀, half-maximal inhibitory concentration; KEGG, Kyoto Encyclopedia of Genes and Genomes; mRNA, messenger RNA; OS, overall survival; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; TME, tumor microenvironment.

scRNA-seq identifies BIRC5 as a key modulator of the TME and therapeutic sensitivity

To comprehensively characterize the cellular heterogeneity and identify key regulatory molecules within the TME, we performed scRNA-seq analysis. Unsupervised clustering and dimensionality reduction using t-distributed stochastic neighbor embedding (t-SNE) visualized the cellular landscape, identifying 11 distinct cell clusters, including B cells, CD4⁺ T cells, CD8⁺ T cells, endothelial cells, epithelial cells, fibroblasts, and myeloid lineage cells (macrophages, neutrophils, etc.) (Figure 8A). We subsequently focused on BIRC5. Assessment of BIRC5 activity using AUCell scoring revealed a heterogeneous expression pattern across the scRNA-seq dataset (Figure 8B). Feature plots further localized high BIRC5 expression to specific proliferative clusters, primarily within the epithelial compartment (Figure 8C). To investigate the biological consequences of BIRC5 heterogeneity, we stratified the samples into BIRC5-high and BIRC5-low groups. Notably, the BIRC5-high group showed an increased epithelial cells fraction (27% vs. 14%) and a relative reduction in CD8⁺ T cells, suggesting a correlation between BIRC5 levels and tumor cell burden or TME infiltration (Figure 8D).

Figure 8 Single-cell and functional validation of BIRC5 linking an immune-excluded TME to enhanced 5-FU resistance in CRC. (A) t-SNE plot of the scRNA-seq dataset showing major annotated cell populations. (B) AUCell-based activity score distribution for BIRC5 program cells above the indicated threshold (AUC >0.23) were classified as BIRC5-high (n=1,609 cells). (C) t-SNE plot colored by AUCell AUC scores. (D) Stacked bar plot comparing the relative cellular composition between high- and low-BIRC5 groups. (E) GSVA comparing pathway activities between BIRC5-high and BIRC5-low groups. (F) Ligand-target inference analysis prioritizing candidate ligands from Treg and CD8⁺ T-cell compartments, with ligand activity and expression in sender cells and predicted downstream target genes in receiver cells; dot size and color indicate regulatory potential and expression metrics, respectively. (G) Network visualization of inferred signaling relationships among key immune regulators and downstream genes. (H) qRT-PCR validation of BIRC5 knockdown efficiency in HCT116 and HT29 cells using two independent shRNAs (Sh1 and Sh2) compared with NC. (I,J) 5-FU dose-response curves and IC₅₀ estimation in HT29 (I) and HCT116 (J) cells following BIRC5 knockdown vs. NC. 5-FU, 5-fluorouracil; AUC, area under the curve; AUPR, area under the precision-recall curve; CRC, colorectal cancer; GSVA, gene set variation analysis; IC₅₀, half-maximal inhibitory concentration; NC, normal control; NK, natural killer; qRT-PCR, quantitative real-time polymerase chain reaction; scRNA-seq, single-cell RNA sequencing; shRNA, short hairpin RNA; t-SNE, t-distributed stochastic neighbor embedding; TME, tumor microenvironment; Treg, regulatory T cell.

To uncover the molecular mechanisms driving the BIRC5-high phenotype, we performed GSVA. This analysis revealed a significant enrichment of oncogenic signaling pathways in the BIRC5-high group, with a particular prominence of the MAPK signaling pathway (Figure 8E). This suggests that BIRC5 upregulation may be intrinsically linked to hyperactive MAPK signaling, contributing to tumor aggressiveness. We next sought to identify the upstream extracellular signals modulating this transcriptional state. NicheNet ligand-receptor analysis prioritized several ligands with high regulatory potential. Notably, IFN-gamma (IFNG) emerged as a top-ranked upstream ligand predicted to drive the observed gene expression changes (Figure 8F). To further elucidate the intracellular regulatory logic, we constructed a target gene-ligand regulatory network. This network highlighted a direct regulatory potential of IFNG on downstream targets, including the chemokine CCL5 and the integrin ITGAL. Crucially, the network structure underscored the involvement of MAPK8 (JNK1) as a central hub connecting IFNG signaling to these downstream effectors (Figure 8G), thereby establishing a putative IFNG-MAPK-BIRC5 regulatory axis within the TME.

Finally, to functionally validate the therapeutic implications of these findings, we targeted BIRC5 in HCT116 and HT29 CRC cell lines. We established stable knockdown models using two independent shRNAs, which effectively depleted BIRC5 mRNA levels (Figure 8H). We then assessed whether BIRC5 suppression could sensitize cancer cells to chemotherapy. Viability assays demonstrated that BIRC5 knockdown significantly potentiated the cytotoxicity of 5-FU. Compared to the negative control (NC), BIRC5-silenced cells exhibited a marked reduction in cell viability and lower IC₅₀ values upon 5-FU treatment (Figure 8I,8J). Together, the single-cell and functional data implicate BIRC5 as a key determinant of a proliferative tumor state associated with dampened immune signaling and increased therapeutic tolerance.

Discussion

Increasing research focuses on the role of SUMOylation in the oncogenesis, progression, and metastasis of CRC (35,36). Many SUMOylation-related proteins are overexpressed in tumor tissues, and drugs targeting these aberrantly expressed components of the SUMOylation process represent potential future cancer therapies. As mentioned earlier, SUMOylation is associated with IFNs and the cGAS/STING pathway, suggesting its involvement in regulating the TME. Furthermore, numerous studies have confirmed that SUMOylation plays a role in modulating the TME pattern (37-39). Thus, therapies targeting SUMOylation might not only directly impair tumor growth and survival but also help overcome tumor heterogeneity and enhance the benefits of immunotherapy. However, research on the relationship between SUMOylation and the CRC TME and its predictive value for treatment strategies remains largely unexplored. Elucidating the role of SUMOylation patterns in the immune infiltration of the CRC TME will help clarify the mechanisms of SUMOylation in anti-CRC immune responses and aid in devising effective therapeutic strategies.

In this study, we initially categorized CRC patients into two distinct clusters based on the expression levels of SRGs. Although there were significant differences in prognosis between the two clusters, with cluster C2 exhibiting higher immune cell infiltration, C2 also showed higher stromal cell infiltration, representing an immune-exclusion phenotype that contradicts the better prognosis observed in C2 patients. Therefore, we concluded that stratifying CRC patients solely based on SRG expression levels was insufficient, and a more in-depth classification was necessary.

To provide more accurate guidance for individual treatment strategies, we developed a quantification system called SUMOscore based on five key SRGs to identify distinct SUMOylation patterns. High SUMOscore is significantly associated with poor prognosis and advanced stages in CRC and are independent prognostic biomarker. Additionally, high SUMOscore correlates with elevated infiltration of TGF-β and EMT pathways, as well as stromal components in the TME. Substantial evidence links the TME, particularly immune and stromal components, to tumor progression and immunotherapy outcomes. Robust infiltration by CD4⁺/CD8⁺ T cells correlates with response, whereas dense stroma can confine immune cells to peri-tumoral nests rather than the parenchyma, thereby dampening antitumor immunity (40). Recent studies also indicate that SUMOylation is associated with the activation of EMT and TGF-β pathways, which may hinder lymphocyte infiltration into the tumor parenchyma (41-45). Our findings corroborate these insights, validating our classification of two SUMOylation immune phenotypes. In the high SUMOscore group, activation of EMT and TGF-β pathways, a higher stromal score, a higher TIDE score, and a lower TMB suggest potential resistance to immunotherapy and chemotherapy (34,46,47). Indeed, in the ICB cohort, the SUMOscore proved valuable in predicting responses to immunotherapy, showing statistically significant differences between responsive and non-responsive groups.

5-FU is an antimetabolite drug where the hydrogen at the C-5 position of uracil is replaced by fluorine. Studies indicate that 5-FU-based combination therapies can significantly prolong the survival in CRC patients (48,49). However, response rates to 5-FU-based chemotherapy remain modest, and chemoresistance frequently undermines efficacy (50,51). Accordingly, identifying and validating biomarkers predictive of 5-FU-based chemotherapy could improve CRC outcomes. Emerging evidence also implicates SUMOylation in multidrug resistance, including resistance to 5-FU (52). In our study, we observed a positive association between SUMOscore and 5-FU IC₅₀, indicating that higher SUMOscore corresponds to reduced 5-FU sensitivity. The specific mechanisms may relate to the higher infiltration of endothelial cells and CAFs in the TME of high SUMOscore group, where the stromal components have been proven to be associated with chemotherapy resistance (53-55). SUMOylation is a central PTM coordinating replication stress and DNA damage response (DDR) signaling, thereby shaping cellular tolerance to genotoxic therapies (56). Notably, SUMOylation of ATRIP enhances assembly and output of the ATR pathway (57), suggesting that an activated SUMO program may promote checkpoint-mediated survival under fluoropyrimidine-induced replication stress; consistent with this concept, ATR inhibition can potentiate 5-FU-associated cytotoxicity in cancer models (58). Among the five genes composing our SUMO-acidification signature, BIRC5 represents a plausible driver. BIRC5 supports mitosis and apoptosis evasion and has been linked to chemoresistance, while targeting BIRC5 can reverse 5-FU resistance in CRC (59,60). Concordantly, we observed that BIRC5 knockdown reduced the IC50 of 5-FU. In parallel, our single-cell analysis indicated that the BIRC5-high subgroup exhibited an epithelial-enriched and CD8⁺ T-cell-deplete TME, and BIRC5 expression can be promoted by RAS-MAPK/ERK signaling, a pathway known to remodel the TME toward immune escape (61,62). Finally, because our score captures acidification-related programs, it is notable that tumor acidosis can impair antitumor immunity, and buffering tumor acidity improves responses to checkpoint blockade in vivo (63). Together, these findings support that the SUMOscore may modulate therapy response through both DDR/replication-stress tolerance and apoptosis regulation and immune activation status of the TME.

To facilitate clinical translation, the SUMOscore should ultimately be converted into a practical and analytically robust assay compatible with routine specimens. One feasible route is to lock the current five-gene algorithm into a small targeted expression panel measurable from formalin-fixed paraffin-embedded (FFPE) biopsy/resection tissue using clinically deployable platforms such as qRT-PCR or hybridization-based digital counting assays (64,65). This development path is supported by precedent in CRC, where multigene assays have been translated from discovery signatures into standardized clinical tests with defined workflows and validated risk stratification (66,67). Key challenges are as follows. First, pre-analytical variability that can affect gene-expression measurements, particularly in FFPE samples (68,69). Second, intratumoral heterogeneity and tumor purity, which may confound bulk transcriptomic readouts and signature stability across sampling regions (26,70). Third, cross-platform calibration and cut-off transferability, requiring rigorous normalization, external controls, and independent cohort validation under REMARK-style reporting frameworks (71). Future work should therefore proceed through analytical validation, followed by retrospective and prospective clinical validation to establish clinical utility and to define how the SUMOscore complements existing clinicopathologic and molecular markers.

This study had several limitations. First, our findings lacked validation from additional clinical datasets beyond public data sources, which would help further substantiate our conclusions. Second, due to the retrospective nature of the ICB cohort and the absence of an appropriate CRC cohort receiving ICB therapy, the level of evidence provided by our study was relatively limited. Third, the current research was confined to mRNA levels, whereas the SUMOylation process depends on protein interactions, introducing potential inaccuracies. Lastly, the specific mechanisms underlying the interaction between SUMOylation patterns and immune cell infiltration in the TME remained unclear and warranted further investigation.

Conclusions

In summary, our work enhanced the understanding of the regulatory mechanisms of SUMOylation signatures on cell infiltration in the CRC TME. Diverse SUMOscore provided a solid foundation for explaining the heterogeneity and complexity of individualized TMEs, guiding more effective strategies in immunotherapy, adjuvant chemotherapy, and targeted treatments.

Acknowledgments

We would like to thank Shanshan Zhang (Ruijin Hospital, Shanghai Jiao Tong University School of Medicine) for assistance with thoughtful discussion.

Footnote

Reporting Checklist: The authors have completed the TRIPOD and MDAR reporting checklists. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2194/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2194/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2194/prf

Funding: This work was supported by the National Natural Science Foundation of China (No. 82373237, to B.F.; No. 82103207, to S.Z.) and the Noncommunicable Chronic Diseases-National Science and Technology Major Project (No. 2023ZD0501600, to B.F.).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2194/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Ruijin Hospital, Shanghai Jiao Tong University School of Medicine [No. 2022 (138)]. Informed consent was taken from all the patients.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Siegel RL, Giaquinto AN, Jemal A. Cancer statistics, 2024. CA Cancer J Clin 2024;74:12-49. [Crossref] [PubMed]
2. Han B, Zheng R, Zeng H, et al. Cancer incidence and mortality in China, 2022. J Natl Cancer Cent 2024;4:47-53. [Crossref] [PubMed]
3. Chen L, Hu H, Yuan Y, et al. CSCO guidelines for colorectal cancer version 2024: Updates and discussions. Chin J Cancer Res 2024;36:233-9. [Crossref] [PubMed]
4. André T, Shiu KK, Kim TW, et al. Pembrolizumab in Microsatellite-Instability-High Advanced Colorectal Cancer. N Engl J Med 2020;383:2207-18. [Crossref] [PubMed]
5. Cercek A, Lumish M, Sinopoli J, et al. PD-1 Blockade in Mismatch Repair-Deficient, Locally Advanced Rectal Cancer. N Engl J Med 2022;386:2363-76. [Crossref] [PubMed]
6. Chao J, Fuchs CS, Shitara K, et al. Assessment of Pembrolizumab Therapy for the Treatment of Microsatellite Instability-High Gastric or Gastroesophageal Junction Cancer Among Patients in the KEYNOTE-059, KEYNOTE-061, and KEYNOTE-062 Clinical Trials. JAMA Oncol 2021;7:895-902. [Crossref] [PubMed]
7. Luke JJ, Rutkowski P, Queirolo P, et al. Pembrolizumab versus placebo as adjuvant therapy in completely resected stage IIB or IIC melanoma (KEYNOTE-716): a randomised, double-blind, phase 3 trial. Lancet 2022;399:1718-29. [Crossref] [PubMed]
8. O'Brien M, Paz-Ares L, Marreaud S, et al. Pembrolizumab versus placebo as adjuvant therapy for completely resected stage IB-IIIA non-small-cell lung cancer (PEARLS/KEYNOTE-091): an interim analysis of a randomised, triple-blind, phase 3 trial. Lancet Oncol 2022;23:1274-86. [Crossref] [PubMed]
9. Diaz LA Jr, Shiu KK, Kim TW, et al. Pembrolizumab versus chemotherapy for microsatellite instability-high or mismatch repair-deficient metastatic colorectal cancer (KEYNOTE-177): final analysis of a randomised, open-label, phase 3 study. Lancet Oncol 2022;23:659-70. [Crossref] [PubMed]
10. Ganesh K, Stadler ZK, Cercek A, et al. Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol 2019;16:361-75. [Crossref] [PubMed]
11. Lizardo DY, Kuang C, Hao S, et al. Immunotherapy efficacy on mismatch repair-deficient colorectal cancer: From bench to bedside. Biochim Biophys Acta Rev Cancer 2020;1874:188447. [Crossref] [PubMed]
12. Weng J, Li S, Zhu Z, et al. Exploring immunotherapy in colorectal cancer. J Hematol Oncol 2022;15:95. [Crossref] [PubMed]
13. Eng C, Kim TW, Bendell J, et al. Atezolizumab with or without cobimetinib versus regorafenib in previously treated metastatic colorectal cancer (IMblaze370): a multicentre, open-label, phase 3, randomised, controlled trial. Lancet Oncol 2019;20:849-61. [Crossref] [PubMed]
14. Westcott PMK, Sacks NJ, Schenkel JM, et al. Low neoantigen expression and poor T-cell priming underlie early immune escape in colorectal cancer. Nat Cancer 2021;2:1071-85. [Crossref] [PubMed]
15. Venne AS, Kollipara L, Zahedi RP. The next level of complexity: crosstalk of posttranslational modifications. Proteomics 2014;14:513-24. [Crossref] [PubMed]
16. Han ZJ, Feng YH, Gu BH, et al. The post-translational modification, SUMOylation, and cancer Int J Oncol 2018;52:1081-94. (Review). [Crossref] [PubMed]
17. Yang Y, He Y, Wang X, et al. Protein SUMOylation modification and its associations with disease. Open Biol 2017;7:170167. [Crossref] [PubMed]
18. Hickey CM, Wilson NR, Hochstrasser M. Function and regulation of SUMO proteases. Nat Rev Mol Cell Biol 2012;13:755-66. [Crossref] [PubMed]
19. Antonczyk A, Krist B, Sajek M, et al. Direct Inhibition of IRF-Dependent Transcriptional Regulatory Mechanisms Associated With Disease. Front Immunol 2019;10:1176. [Crossref] [PubMed]
20. Liu S, Yang B, Hou Y, et al. The mechanism of STING autoinhibition and activation. Mol Cell 2023;83:1502-1518.e10. [Crossref] [PubMed]
21. Samson N, Ablasser A. The cGAS-STING pathway and cancer. Nat Cancer 2022;3:1452-63. [Crossref] [PubMed]
22. Charoentong P, Finotello F, Angelova M, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18:248-62. [Crossref] [PubMed]
23. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 2017;18:220. [Crossref] [PubMed]
24. Yeo JG, Wasser M, Kumar P, et al. The Extended Polydimensional Immunome Characterization (EPIC) web-based reference and discovery tool for cytometry data. Nat Biotechnol 2020;38:679-84. [Crossref] [PubMed]
25. Becht E, Giraldo NA, Lacroix L, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 2016;17:218. [Crossref] [PubMed]
26. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
27. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
28. Ghandi M, Huang FW, Jané-Valbuena J, et al. Next-generation characterization of the Cancer Cell Line Encyclopedia. Nature 2019;569:503-8. [Crossref] [PubMed]
29. Rees MG, Seashore-Ludlow B, Cheah JH, et al. Correlating chemical sensitivity and basal gene expression reveals mechanism of action. Nat Chem Biol 2016;12:109-16. [Crossref] [PubMed]
30. Yang C, Huang X, Li Y, et al. Prognosis and personalized treatment prediction in TP53-mutant hepatocellular carcinoma: an in silico strategy towards precision oncology. Brief Bioinform 2021;22:bbaa164. [Crossref] [PubMed]
31. Hugo W, Zaretsky JM, Sun L, et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell 2016;165:35-44. [Crossref] [PubMed]
32. Cho JW, Hong MH, Ha SJ, et al. Genome-wide identification of differentially methylated promoters and enhancers associated with response to anti-PD-1 therapy in non-small cell lung cancer. Exp Mol Med 2020;52:1550-63. [Crossref] [PubMed]
33. Rose TL, Weir WH, Mayhew GM, et al. Fibroblast growth factor receptor 3 alterations and response to immune checkpoint inhibition in metastatic urothelial cancer: a real world experience. Br J Cancer 2021;125:1251-60. [Crossref] [PubMed]
34. Mariathasan S, Turley SJ, Nickles D, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 2018;554:544-8. [Crossref] [PubMed]
35. Peng C, Tan Y, Yang P, et al. Circ-GALNT16 restrains colorectal cancer progression by enhancing the SUMOylation of hnRNPK. J Exp Clin Cancer Res 2021;40:272. [Crossref] [PubMed]
36. Liu T, Wang H, Chen Y, et al. SENP5 promotes homologous recombination-mediated DNA damage repair in colorectal cancer cells through H2AZ deSUMOylation. J Exp Clin Cancer Res 2023;42:234. [Crossref] [PubMed]
37. Wu Z, Huang H, Han Q, et al. SENP7 senses oxidative stress to sustain metabolic fitness and antitumor functions of CD8+ T cells. J Clin Invest 2022;132:e155224. [Crossref] [PubMed]
38. Xia QD, Sun JX, Xun Y, et al. SUMOylation Pattern Predicts Prognosis and Indicates Tumor Microenvironment Infiltration Characterization in Bladder Cancer. Front Immunol 2022;13:864156. [Crossref] [PubMed]
39. Gu Y, Fang Y, Wu X, et al. The emerging roles of SUMOylation in the tumor microenvironment and therapeutic implications. Exp Hematol Oncol 2023;12:58. [Crossref] [PubMed]
40. Topalian SL, Taube JM, Anders RA, et al. Mechanism-driven biomarkers to guide immune checkpoint blockade in cancer therapy. Nat Rev Cancer 2016;16:275-87. [Crossref] [PubMed]
41. Li SN, Wu JF. TGF-β/SMAD signaling regulation of mesenchymal stem cells in adipocyte commitment. Stem Cell Res Ther 2020;11:41. [Crossref] [PubMed]
42. Brabletz S, Schuhwerk H, Brabletz T, et al. Dynamic EMT: a multi-tool for tumor progression. EMBO J 2021;40:e108647. [Crossref] [PubMed]
43. Chanda A, Ikeuchi Y, Karve K, et al. PIAS1 and TIF1γ collaborate to promote SnoN SUMOylation and suppression of epithelial-mesenchymal transition. Cell Death Differ 2021;28:267-82. [Crossref] [PubMed]
44. Tauriello DVF, Palomo-Ponce S, Stork D, et al. TGFβ drives immune evasion in genetically reconstituted colon cancer metastasis. Nature 2018;554:538-43. [Crossref] [PubMed]
45. Bai X, Yi M, Jiao Y, et al. Blocking TGF-β Signaling To Enhance The Efficacy Of Immune Checkpoint Inhibitor. Onco Targets Ther 2019;12:9527-38. [Crossref] [PubMed]
46. Voron T, Marcheteau E, Pernot S, et al. Control of the immune response by pro-angiogenic factors. Front Oncol 2014;4:70. [Crossref] [PubMed]
47. Vértessy BG, Tóth J. Keeping uracil out of DNA: physiological role, structure and catalytic mechanism of dUTPases. Acc Chem Res 2009;42:97-106. [Crossref] [PubMed]
48. Grothey A, Sargent D, Goldberg RM, et al. Survival of patients with advanced colorectal cancer improves with the availability of fluorouracil-leucovorin, irinotecan, and oxaliplatin in the course of treatment. J Clin Oncol 2004;22:1209-14. [Crossref] [PubMed]
49. Van Cutsem E, Rivera F, Berry S, et al. Safety and efficacy of first-line bevacizumab with FOLFOX, XELOX, FOLFIRI and fluoropyrimidines in metastatic colorectal cancer: the BEAT study. Ann Oncol 2009;20:1842-7. [Crossref] [PubMed]
50. Vodenkova S, Buchler T, Cervena K, et al. 5-fluorouracil and other fluoropyrimidines in colorectal cancer: Past, present and future. Pharmacol Ther 2020;206:107447. [Crossref] [PubMed]
51. Sethy C, Kundu CN. 5-Fluorouracil (5-FU) resistance and the new strategy to enhance the sensitivity against cancer: Implication of DNA repair inhibition. Biomed Pharmacother 2021;137:111285. [Crossref] [PubMed]
52. Qin Y, Bao H, Pan Y, et al. SUMOylation alterations are associated with multidrug resistance in hepatocellular carcinoma. Mol Med Rep 2014;9:877-81. [Crossref] [PubMed]
53. Schuth S, Le Blanc S, Krieger TG, et al. Patient-specific modeling of stroma-mediated chemoresistance of pancreatic cancer using a three-dimensional organoid-fibroblast co-culture system. J Exp Clin Cancer Res 2022;41:312. [Crossref] [PubMed]
54. Vienot A, Pallandre JR, Renaude E, et al. Chemokine switch regulated by TGF-β1 in cancer-associated fibroblast subsets determines the efficacy of chemo-immunotherapy. Oncoimmunology 2022;112144669.
55. Zeng Q, Mousa M, Nadukkandy AS, et al. Understanding tumour endothelial cell heterogeneity and function from single-cell omics. Nat Rev Cancer 2023;23:544-64. [Crossref] [PubMed]
56. Cremona CA, Sarangi P, Zhao X. Sumoylation and the DNA damage response. Biomolecules 2012;2:376-88. [Crossref] [PubMed]
57. Wu CS, Ouyang J, Mori E, et al. SUMOylation of ATRIP potentiates DNA damage signaling by boosting multiple protein interactions in the ATR pathway. Genes Dev 2014;28:1472-84. [Crossref] [PubMed]
58. Ito SS, Nakagawa Y, Matsubayashi M, et al. Inhibition of the ATR kinase enhances 5-FU sensitivity independently of nonhomologous end-joining and homologous recombination repair pathways. J Biol Chem 2020;295:12946-61. [Crossref] [PubMed]
59. Wheatley SP, Altieri DC. Survivin at a glance. J Cell Sci 2019;132:jcs223826. [Crossref] [PubMed]
60. Cao Y, Tang H, Wang G, et al. Targeting survivin with Tanshinone IIA inhibits tumor growth and overcomes chemoresistance in colorectal cancer. Cell Death Discov 2023;9:351. [Crossref] [PubMed]
61. Chang WH, Liu Y, Hammes EA, et al. Oncogenic RAS promotes MYC protein stability by upregulating the expression of the inhibitor of apoptosis protein family member Survivin. J Biol Chem 2023;299:102842. [Crossref] [PubMed]
62. Hamarsheh S, Groß O, Brummer T, et al. Immune modulatory effects of oncogenic KRAS in cancer. Nat Commun 2020;11:5439. [Crossref] [PubMed]
63. Pilon-Thomas S, Kodumudi KN, El-Kenawi AE, et al. Neutralization of Tumor Acidity Improves Antitumor Responses to Immunotherapy. Cancer Res 2016;76:1381-90. [Crossref] [PubMed]
64. Bustin SA, Ruijter JM, van den Hoff MJB, et al. MIQE 2.0: Revision of the Minimum Information for Publication of Quantitative Real-Time PCR Experiments Guidelines. Clin Chem 2025;71:634-51. [Crossref] [PubMed]
65. Li J, Fu C, Speed TP, et al. Accurate RNA Sequencing From Formalin-Fixed Cancer Tissue To Represent High-Quality Transcriptome From Frozen Tissue. JCO Precis Oncol 2018;2018:PO.17.00091.
66. Yothers G, O'Connell MJ, Lee M, et al. Validation of the 12-gene colon cancer recurrence score in NSABP C-07 as a predictor of recurrence in patients with stage II and III colon cancer treated with fluorouracil and leucovorin (FU/LV) and FU/LV plus oxaliplatin. J Clin Oncol 2013;31:4512-9. [Crossref] [PubMed]
67. Maak M, Simon I, Nitsche U, et al. Independent validation of a prognostic genomic signature (ColoPrint) for patients with stage II colon cancer. Ann Surg 2013;257:1053-8. [Crossref] [PubMed]
68. Agrawal L, Engel KB, Greytak SR, et al. Understanding preanalytical variables and their effects on clinical biomarkers of oncology and immunotherapy. Semin Cancer Biol 2018;52:26-38. [Crossref] [PubMed]
69. Lin Y, Dong ZH, Ye TY, et al. Optimization of FFPE preparation and identification of gene attributes associated with RNA degradation. NAR Genom Bioinform 2024;6:lqae008. [Crossref] [PubMed]
70. Dunne PD, Alderdice M, O'Reilly PG, et al. Cancer-cell intrinsic gene expression signatures overcome intratumoural heterogeneity bias in colorectal cancer patient classification. Nat Commun 2017;8:15657. [Crossref] [PubMed]
71. Altman DG, McShane LM, Sauerbrei W, et al. Reporting Recommendations for Tumor Marker Prognostic Studies (REMARK): explanation and elaboration. PLoS Med 2012;9:e1001216. [Crossref] [PubMed]