Prognostic models of immune-related cell death and stress unveil mechanisms driving macrophage phenotypic evolution in colorectal cancer

Abstract

Background

Tumor microenvironment (TME), particularly immune cell infiltration, programmed cell death (PCD) and stress, has increasingly become a focal point in colorectal cancer (CRC) treatment. Uncovering the intricate crosstalk between these factors can enhance our understanding of CRC, guide therapeutic strategies, and improve patient prognosis.

Methods

We constructed an immune-related cell death and stress (ICDS) prognostic model utilizing machine learning methodologies. Furthermore, we performed enrichment analyses and deconvolution algorithms to elucidate the complex interactions between immune cell infiltration and the processes of PCD and stress within a substantial array of transcriptomic data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus data base (GEO) related to CRC. Single-cell sequencing and biochemical experiments were used to validate the interaction between the model genes and programmed cell death in tumor cells.

Results

The ICDS prognostic model exhibited robust predictive performance in seven independent cohorts, revealing an inverse correlation between model scores and patient prognosis. Meanwhile, the ICDS index was positively correlated with clinical stage. Model analysis indicated that patient subgroups with low ICDS index exhibited heightened immune activation features and elevated activity in PCD and stress pathways. Single-cell analysis further revealed that macrophages were the central drivers of immune characteristics underlying prognostic differences within the ICDS prognostic model. Pseudotime analysis and cellular experiments indicated that the model gene GAL3ST4 promotes the transition of macrophages toward an M2 pro-tumor phenotype. Furthermore, cell communication analysis and experimental validation revealed that the cuproptosis in tumor cells suppress GAL3ST4 expression, thereby inhibiting M2-like macrophage polarization.

Conclusion

In summary, we constructed the ICDS prognostic model and uncovered the mechanism by which tumor cells downregulate GAL3ST4 expression via cuproptosis to inhibit M2-like macrophage polarization, providing new targets and biomarkers for CRC treatment and prognosis evaluation.

Introduction

Colorectal cancer (CRC) is the third most common cancer worldwide and the second leading cause of cancer-related deaths globally [1]. The treatment for CRC currently includes a variety of methods, such as chemotherapy, targeted therapy and immunotherapy [2, 3]. However, many patients fail to benefit from these therapies. Consequently, scientists worldwide are committed to exploring innovative diagnostic and therapeutic strategies and developing novel drugs to expand the range of treatment options [4,5,6]. Notably, cancer treatments can induce various types of programmed cell death (PCD), which significantly influence therapeutic outcomes through distinct mechanisms [7, 8]. However, these mechanisms remain poorly studied and insufficiently understood.

Tumor microenvironment (TME) is composed of various immune cells, including macrophages, fibroblasts, T cells, and neutrophils, along with cytokines and biochemical components such as the extracellular matrix, providing critical support for tumor initiation and progression [9,10,11]. The remodeling of TME is often accompanied by significant alterations in the states of immune cells, profoundly impacting cancer treatment outcomes, particularly immunotherapy [9, 12].The remodeling of TME is regulated by multiple factors, with recent studies revealing extensive crosstalk between the tumor immune microenvironment, PCD and cellular stress [13,14,15]. Therefore, elucidating the mechanisms by which PCD and cellular stress contribute to immune microenvironment remodeling will provide fresh perspectives into how PCD influences cancer therapy.

Certain forms of PCD, such as pyroptosis and necroptosis, are referred to as immunogenic cell death due to their ability to significantly impact immune cell states and reshape TME [16]. In CRC, TME remodeling induced by immunogenic cell death not only enhances the efficacy of chemotherapy and radiotherapy but also significantly increases the response rate of microsatellite-stable (MSS) colorectal cancer patients to immunotherapy [17,18,19]. Despite recent advances in understanding the role of immunogenic cell death in regulating TME remodeling in CRC, such as immunogenic cell death triggering immune responses and reshaping the immune microenvironment through the release of damage-associated molecular patterns (DAMPs) interacting with pattern recognition receptors (PRRs) on immune cells [20,21,22,23], the specific mechanisms of different forms of immunogenic cell death in this process and their precise interactions with immune cell subsets in CRC remain largely unknown.

In this study, we developed a prognostic model that integrates immunogenic cell death features and cellular stress characteristics using immune cell infiltration data. The model demonstrated excellent predictive performance across multiple independent cohorts, with model scores significantly negatively correlated with patient prognosis while revealing differences in immune responses and immunogenic cell death-related biological processes among patients with varying prognoses. Single-cell transcriptomics analysis revealed a strong association between model grouping and differences in GAL3ST4 expression levels as well as macrophage polarization states. Further experimental validation demonstrated that GAL3ST4 promoted the polarization of macrophages toward an M2 pro-tumor phenotype, a process regulated by tumor cell-mediated cuproptosis. By integrating multi-omics data and experimental analyses, we discovered that tumor cells suppress GAL3ST4 expression through the process of cuproptosis, thereby inhibiting the polarization of M2-like macrophages. These findings provided critical theoretical foundations and potential applications for understanding the mechanisms of CRC immune microenvironment remodeling and related therapeutic strategies.

Methods

Processing of Bulk RNA-seq data

CRC RNA-seq data and clinical information from the TCGA database were retrieved via the UCSC Xena browser (http://xena.ucsc.edu/). Gene expression data obtained from the GEO database (GSE29621, GSE39582, GSE17536, GSE17537, GSE33382, and GSE72970) were processed with log2(x + 1) transformation to acquire samples with normalized RNA-seq or microarray data.

Single-cell RNA-seq data collection and analysis

A single-cell RNA sequencing (scRNA-seq) dataset comprising 23 patients was obtained from the GEO database (GSE132465). Cells with > 1,000 unique molecular identifier (UMI) counts > 200 genes < 6,000 genes, and < 10% of mitochondrial gene expression in UMI counts were filtered. Following the standard dimensionality reduction and clustering procedures, clusters were annotated with classical cell markers.

For the sub-clustering analysis of myeloid cells and macrophages, we extracted cell populations initially annotated as myeloid cells and macrophages and repeated the aforementioned workflow. Myeloid cells were annotated based on classical myeloid markers, while macrophages were further annotated and classified according to their function.

To differentiate between malignant and non-malignant epithelial cells, the ‘InferCNV’ R package was used to deduce large-scale chromosomal copy-number variations (CNVs) from single-cell transcriptome profiles. The R package ‘CellChat’ with default parameters was employed to compute the communication intensity between cuproptosis malignant epithelial cells and other cell types, while the ‘Monocle2’ package, also with default settings, was used to delineate the phenotypic transition trajectory of macrophages.

Construction of the ICDS model gene set

Selection of Immune-Related Genes: single-sample gene set enrichment analysis (ssGSEA) via the R package GSVA was employed to quantify the relative infiltration of 28 immune cell types in the TCGA-CRC cohort (97 genes for 28 immune cells downloaded from the Molecular Signatures Database (MSigDB)). Subsequently, clustering was performed using the ‘Consensus Cluster Plus’R package based on the infiltration profiles of these immune cells. Thereafter, a comprehensive assessment utilizing the consensus score matrix and CDF curve was conducted to determine the optimal number of clusters. The ‘WGCNA’R package was then employed to perform unsigned weighted gene co-expression network analysis (WGCNA) to identify immune-related gene modules (n = 4295).

Selection of programmed cell death and stress genes: The R package ‘GSVA’ was used for ssGSEA to evaluate the enrichment of four pathways associated with programmed cell death and stress responses in the TCGA-CRC cohort. Clustering was then carried out with the ‘Consensus Cluster Plus’ R package based on the enrichment data of these pathways. The method for determining the optimal number of clusters is the same as described above. Finally, differential analysis was performed using the R package ‘limma’ (|logFC|> 0.2), resulting in the identification of 2,769 differentially expressed genes associated with programmed cell death and stress.

Construction of the ICDS Model Gene Set: The final ICDS model gene set (n = 579) was established as the intersection of programmed cell death and stress genes with immune-related genes, deemed a reliable representation.

Construction of the ICDS machine learning model

We employed a three-step strategy to construct the ICDS machine learning model, with each patient's final score considered the ICDS index. Step 1: Univariate Cox Regression: Prognostic genes were identified through univariate Cox regression analysis, defining genes with hazard ratios (HR) and 95% confidence intervals (CI) with p < 0.05. Step 2: LASSO Cox Regression: To prevent overfitting, LASSO Cox regression was applied using the ‘glmnet’ R package to narrow down the candidate genes. Genes with the minimal lambda value were retained in the model. Step 3: Stepwise Cox Regression: Stepwise Cox regression was ultimately used to select the genes most suitable for model construction. The model with the minimal Akaike Information Criterion (AIC) for survival genes was designated as the final:

$$ICDS={\sum }_{i=1}^{n}Coefi*Expi$$

where n denotes the number of model genes, Coefi denotes each model gene risk coefficient, and Expi denotes each model gene mRNA level.

Prognostic evaluation of ICDS index

Kaplan–Meier survival curves were generated using the ‘survival’ package in R, followed by univariate and multivariate Cox regression analyses. Forest plots were created using the ‘forestplot’ package. Clinical nomograms were constructed utilizing the ‘survival’ and ‘rms’ packages.

Analysis of immune cell and cytokine infiltration

Utilizing the R package ‘GSVA’, ssGSEA was employed to quantify the relative abundance of various TME-infiltrating cells and cytokines within ICDS subgroups [24].

Predictive analysis of immunotherapy

Utilizing the R package ‘easier’, we computed immunotherapy response scores for ICDS subgroups based on cellular components, pathway activity, transcription factor activity, ligand-receptor weights, and intercellular interactions.

Pathway enrichment analysis

Through the R package ‘clusterProfiler’, pathways from the Gene Ontology (GO) Biological Process Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) were subjected to over-representation analysis (ORA) and gene set enrichment analysis (GSEA). Differential genes in the single-cell dataset were identified using the FindAllMarkers function (p < 0.01, avg_Log2FC ≥ 0.25), while bulk transcriptomic differential genes were analyzed with the R package ‘limma’ (p < 0.05, |logFC|< 0.5).

Single-cell ICDS differential gene scoring

Quantitative analysis of differential gene expression between ICDS subgroups in the TCGA-CRC cohort was conducted using the R package ‘AUCell’. Visualization of cell type-specific scores was achieved through the use of violin plots and ridge plots, produced by the R packages ‘ggplot2’ and ‘ggridge’, respectively.

Cell culture

The human colon cancer cell lines DLD1 and HCT116 were obtained from the American Type Culture Collection (ATCC). The HCT116 cell line was cultured in DMEM medium (Gibco, Carlsbad, USA) supplemented with 10% fetal bovine serum, while the DLD1 cell line was maintained in RPMI-1640 medium (Gibco, Carlsbad, USA) with 10% fetal bovine serum. THP-1 cells were purchased from Procell (Wuhan, China). THP-1 cells were maintained in RPMI 1640 medium (Gibco, Carlsbad, USA) augmented with 0.05 mM 2-mercaptoethanol (Sigma-Aldrich) and 10% fetal bovine serum. All cell lines were kept in a cell incubator at 37 °C with a 5% CO2 atmosphere.

Macrophage induction and polarization

Initially, THP-1 cells were exposed to 100 nM PMA (Sigma-Aldrich) fo 24 h to facilitate their differentiation into M0 phase macrophages. Subsequently, the M0 macrophage culture medium was supplemented with 20 ng/ml IL-4 (PeproTech) and 20 ng/ml IL-13 (PeproTech), followed by a 48-h incubation to induce polarization into M2 macrophages.

Collection of cuproptosis cell fragments

The HCT116 and DLD1 colon cancer cell lines were uniformly plated in six-well plates and cultured until they reached 80% confluence. Then, 0.5 µM Elesclomol (Selleck, S1052) was added, and the cells were incubated for 48 h. After 48 h, extensive cell death was observed. The cell culture medium was then subjected to centrifugation at 1000 g for 10 min, and the supernatant was collected for co-culturing with M2 macrophages.

Real-time quantitative PCR (qRT-PCR) assay

The reverse transcription of total RNA was performed utilizing Hifair® III Reverse Transcriptase (Yeasen Biotechnology Co., Ltd., Shanghai, China) according to the manufacturer's specifications. Subsequently, the resultant product was immediately stored at −80 °C for future use. Quantitative real-time polymerase chain reaction (qRT-PCR) was conducted using a QuantStudio 3 Real-Time PCR System (Thermo Fisher) in conjunction with Hieff® qPCR SYBR Green Master Mix (Yeasen Biotechnology Co., Ltd., Shanghai, China). The qRT-PCR protocol entailed an initial denaturation step at 95 °C for 5 min, followed by 40 amplification cycles, each comprising denaturation at 95 °C for 10 s, primer-specific annealing at 60 °C for 20 s, and extension at 72 °C for 20 s. Primer sequences are provided in Table S1. The relative expression levels of the target genes were determined using the 2 − ΔΔCt method, with GAPDH employed as the internal reference gene.

Western blot

Proteins were extracted from cell lines using RIPA buffer, and equivalent protein samples (30 μg, quantified by the BCA Protein Assay Kit) were separated on SDS-PAGE gels. The proteins were subsequently transferred onto PVDF membranes. Post-transfer, the membranes were blocked with TBST containing 5% BSA, followed by overnight incubation at 4 °C with primary antibodies targeting GAL3ST4 (1:3000). The membranes were then hybridized with secondary antibodies at 37 °C for 1 h, after which band intensities were detected using the ECL-plus kit.

Immunofluorescence of cells

Inoculate the prepared cells into a 24-well plate with coverslips, cultivating until a suitable cell density is attained; aspirate the culture medium and wash with pre-cooled PBS for 5 min, repeated three times; fix the cells with 4% paraformaldehyde, followed by washing with PBS for 5 min, repeated three times; treat with a permeabilization solution containing 0.1% Triton X-100 for 15 min, followed by washing with PBS for 5 min, repeated three times; treat with a permeabilization solution containing 0.1% Triton X-100 for 15 min, followed by washing with PBS for 5 min, repeated three times; aspirate the blocking solution and add the primary antibody (1:200), incubating overnight at 4 °C; the next day, aspirate the primary antibody, wash with PBS for 5 min, repeated three times, then add the fluorescent secondary antibody, incubating in the dark at room temperature for 1 h, followed by PBS washing for 5 min, repeated three times; use a pipette tip to apply an antifade mounting medium containing DAPI onto a glass slide, and gently place the coverslip on top; observe the cells under a fluorescence microscope.

Immunofluorescence Assay for Tissue Sections

Subject the paraffin sections to a deparaffinization agent (Wuhan, Seville) until they achieve complete hydration; apply an antigen retrieval buffer (Wuhan, Seville) to the prepared tissue sections to restore antigenicity, subsequently rinsing with PBS for three repetitions of 5 min; administer a 3% BSA solution onto the sections for a blocking period of 30 min; prepare a mixed solution containing two primary antibodies targeting distinct protein markers, dispense the solution onto the sections, and incubate overnight at 4 °C in a humidified chamber, followed by washing with PBS for three cycles of 5 min each; introduce the appropriate secondary antibodies, incubating for 1 h at room temperature while shielded from light, then perform PBS washes for three repetitions of 5 min each; apply DAPI staining solution to the sections and incubate at room temperature for 10 min, subsequently washing with PBS for three repetitions of 5 min; apply DAPI staining solution to the sections and incubate at room temperature for 10 min, subsequently washing with PBS for three repetitions of 5 min; utilize an antifade mounting agent for the sections, and subsequently acquire images with a fluorescence scanner (Pannoramic MIDI, 3DHISTECH).

Immunohistochemistry

Immerse paraffin sections in a deparaffinization solution (Wuhan, Seville) until fully rehydrated; apply an antigen retrieval solution (Wuhan, Seville) to the processed tissue sections for antigen restoration, followed by three washes with PBS for 5 min each; incubate the sections in a 3% hydrogen peroxide solution at room temperature, shielded from light, for 25 min to inhibit endogenous peroxidase activity, followed by three washes with PBS for 5 min each; incubate the sections in a 3% hydrogen peroxide solution at room temperature, shielded from light, for 25 min to inhibit endogenous peroxidase activity, followed by three washes with PBS for 5 min each; incubate the sections in a 3% hydrogen peroxide solution at room temperature, shielded from light, for 25 min to inhibit endogenous peroxidase activity, followed by three washes with PBS for 5 min each; introduce the corresponding secondary antibody and incubate at room temperature for 1 h, followed by three washes with PBS for 5 min each; dispense freshly prepared DAB chromogen, controlling the color development time under a microscope, with positive staining appearing brown, and subsequently rinse with tap water to halt the reaction; subject the sections to hematoxylin counterstaining for approximately 3 min, followed by rinsing under running water; then apply a differentiation solution for a few seconds, rinse again, and use a bluing solution followed by another rinse; sequentially immerse the sections in 75% ethanol for 5 min, 85% ethanol for 5 min, absolute ethanol I for 5 min, absolute ethanol II for 5 min, n-butanol for 5 min, and xylene I for 5 min for dehydration and clearing; then remove the sections from xylene, allow them to air dry briefly, and mount with adhesive; observe the sections under a bright-field microscope.

Cell transfection

Transfection is conducted when the cell density reaches 70%. siRNA is directly mixed with the transfection reagent (Advanced DNA RNA Transfection Reagent, AD600150) in a 1:1 ratio. The mixture is gently pipetted 10–15 times to ensure thorough homogenization and left to stand at room temperature for 10–15 min. The nucleic acid complexes are then added to the cell culture medium for transfection. After 48 h, the mRNA expression levels are assessed using PCR. The si-RNA (si-GAL3ST4) sequences are detailed in Table S2.

Cell coculture and cell metastasis assays

Thereafter, DLD1 and HCT116 cells were introduced into Transwell chambers (8-μm pore membrane, Corning, NY, USA) at a density of 2 × 10^5 cells per chamber, with varying treatments applied according to experimental requirements. Following a 48-h incubation, the chambers were fixed using 4% paraformaldehyde, and cells adhering to the upper side of the membrane were carefully removed by gentle wiping. The membranes were subsequently stained with crystal violet and examined under an optical microscope.

The CCK8 assay

The cells were initially seeded into 96-well culture plates at a density of 2000 cells per well. After a 24-h incubation, the cells underwent various treatments. Cell viability was evaluated at 0, 24, 48, 72, and 96 h post-treatment using the CCK8 assay, following the manufacturer's protocol.

Human samples

For the internal CRC sample cohort, three patients with a CRC diagnosis were recruited, and their clinical data were presented in Table S3. Tissue samples were obtained from the First Affiliated Hospital of Zhengzhou University, and all procedures were performed in accordance with the approval granted by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University (Approval No: 2023-KY-0353-003). All participating patients provided informed consent.

Statistical Analysis

All data analyses were conducted using R software version 4.2.2. Pearson's correlation test was used for correlation analysis. Comparisons between the two groups were calculated using unpaired Student’s t-test or unpaired nonparametric Mann–Whitney test. Each assay was independently replicated more than three times, with results displayed as means ± standard error of the mean (SEM) in the figures. Multiple time point comparisons between two groups, as well as comparisons among more than two groups, were computed using analysis of variance (ANOVA).

Results

Construction of the immune-related cell death and stress (ICDS) prognosis model

To identify the immune-related cell death and stress (ICDS) genes, we initially conducted enrichment scoring of immune cells, programmed cell death (PCD) and stress pathways using the Gene Set Variation Analysis (GSVA) algorithm. Subsequently, the Consensus Clustering method was used to distinguish between high-expression and low-expression groups (Fig. S1A-F). Following this, the Weighted Correlation Network Analysis (WGCNA) (Fig. S1G-J) was utilized to screen for cell death and stress related genes, while Linear Models for Microarray Data (Limma) was applied to identify immune-related genes (Fig. S1K). This analysis revealed a total of 4,295 immune-related genes and 2,769 cell death and stress related genes. By filtering for genes common to both categories, we identified 579 the ICDS genes (Fig. 1A). To further screen clinically relevant ICDS genes, we performed univariate Cox analysis and Lasso regression analysis (Fig. 1B-C), resulting in the identification of 29 ICDS genes. Then, based on Stepwise Cox Proportional Hazards Regression, the ICDS model composed of 9 key genes was finally determined, and the ICDS index was calculated (Fig. 1D). Univariate Cox analysis of overall survival time in the TCGA-CRC cohort indicated that CXXC5, ATP6V1B2, FAS, and MMP1 were associated with improved survival for colorectal cancer patients. Conversely, BST2, TRIB2, ENO2, GAL3ST4, and FLT1 were identified as adverse factors affecting overall survival (Fig. 1E).

Fig. 1

Construction of the immune-related cell death and stress (ICDS) prognostic model and its clinical features. A Gene selection process for the ICDS model. B LASSO coefficient curve for 29 predictor variables. C Identification of the optimal lambda penalty in the LASSO regression model. D Nine model genes and their coefficients selected through LASSO Cox regression analysis. E Forest plot showing the univariate Cox values of the model genes. Red dots indicate negative effects, blue dots indicate positive effects. The points represent estimates, and the error bars indicate 95% confidence intervals. F–H Boxplots showing the association between the ICDS index and variations in clinical indicators within the TCGA-CRC cohort. I  Display the differential expression of signature genes across ICDS groups in the TCGA-CRC cohort, with the ICDS index and clinical factors presented as figure legends. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001

The clinical staging results from the TCGA-CRC cohort demonstrated a significant correlation between the ICDS index and the T, N and the clinical stage. Specifically, in the T stage, patients in the T3 and T4 categories exhibited higher ICDS index (Fig. 1F). In the N stage, individuals with N1 and N2 classifications also showed elevated ICDS index (Fig. 1G). Furthermore, in the clinical stage, patients in stages III and IV had higher ICDS index compared to those in stages I and II (Fig. 1H), which were linked to poorer clinical outcomes. The relationship between these clinical traits and the ICDS index had also been corroborated by other datasets (Fig. S2B-E). Among these cohorts, no significant correlation between the ICDS index and the M stage was detected (Fig. S2A). A further analysis was conducted to investigate the differential expression of modeling genes between the high and low ICDS index groups within the TCGA-CRC cohort. The results revealed that, except for CXXC5, all other genes exhibited significant differences in expression (Fig. 1I). Moreover, our findings reveal correlations between the ICDS index and model genes with T stage, N stage, and clinical stage, suggesting their potential utility in clinical prognosis. (Fig. 1I and Fig. S2F-G).

Prognostic value of the ICDS index in multiple colorectal cancer cohorts

To assess the applicability of the ICDS index in predicting the survival probability of CRC patients, we categorized patients into high ICDS index and low ICDS index subgroups based on quartiles across seven datasets. Kaplan–Meier survival analysis revealed that CRC patients in the high ICDS index subgroup exhibited a significantly elevated risk of death across all cohorts (p < 0.05) (Fig, 2A-F right and Fig. S2H). The comparison of receiver operating characteristic (ROC) curves revealed that the ICDS index showed strong predictive accuracy for overall survival in CRC cohorts at 1, 3, and 5 years. (Fig, 2A-F left and Fig. S2I). We performed both univariate and multivariate Cox analyses to thoroughly evaluate the prognostic significance of the ICDS index within the TCGA-CRC cohort. The results of both indicate that the ICDS index is an independent risk factor in the TCGA-CRC cohort (p < 0.05). (Fig. 2G and Fig. S2J). In addition, we constructed a predictive nomogram model that integrates multiple clinical factors, including the ICDS index, to predict overall survival (OS) at 1, 3, and 5 years (Fig. 2H). The C-index of the nomogram models was 0.820 (95% CI 0.70–0.872). The bootstrap-corrected calibration curve showed that the nomogram model for overall survival exhibits good predictive efficiency (Fig. S2K). In summary, our findings indicated that the ICDS index could robustly predict the prognosis of patients with CRC.

Figure2

Prognostic validation of the ICDS model. A-F ROC curves (1-, 3-, and 5-year) and overall survival Kaplan–Meier curves for high and low ICDS index groups across datasets (TCGA-CRC, GSE29621, GSE38832, GSE39582, GSE17536, GSE17537). G Forest plot showing the univariate Cox regression hazard ratio (HR) and 95% confidence interval (CI) for overall survival, based on clinical characteristics and ICDS index. H Nomogram developed to predict overall survival probability based on clinical characteristics and ICDS index. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001

Immune cell infiltration, programmed cell death and cell stress characteristics of the ICDS index

We evaluated the impact of the ICDS index on immune cell infiltration based on our model features. In comparison to the high-ICDS index group, the low-ICDS index group exhibited increased levels of activated CD4 + T cells, activated dendritic cells (DCs), central memory CD8 T cells, effector memory CD4 T cells, cytotoxic T cells, macrophages, mast cells, myeloid-derived suppressor cells (MDSCs), as well as helper T cells of types 1, 2, and 17 (Fig. 3A). For these cells exhibiting significant differences, aside from MDSCs, which might play a role in immunosuppression, the remaining cells primarily contributed to immune antigen presentation, phagocytosis, or direct tumor cell killing, while also promoting immunity through the secretion of various cytokines. This phenomenon might suggest that the low-ICDS index subgroup exhibited a highly activated immune state, which could be associated with its favorable prognosis. We further analyzed the cytokine level differences between the two subgroups, and the results showed significantly elevated expression levels of pro-inflammatory cytokines, including IL1B, IL6, IL18, TNF, OSM, CSF2, CCL2, CCL3, CCL4, and IL10, in the subgroup with a lower ICDS index (Fig. 3B). These findings further indicated that the subgroup with a lower ICDS index displayed a heightened immune-promoting state, consistent with the observed differences in immune cell profiles. In the comparison of immune checkpoints, ICOSLG (CD275) and SIGLEC15 (CD33) were found to be up-regulated in the high ICDS index group. Conversely, the low ICDS index group demonstrated high expression of immune-suppressive checkpoints such as CD274, CD70, CTLA4, PDCD1LG2 (PD-L2), and IDO1, as well as immune-promoting checkpoints including ICOS and TNFRSF9 (CD137) (Fig. 3C). These findings indicated that the two patient groups had disparate responses to immunotherapy, with the low-ICDS index group likely exhibiting superior therapeutic effectiveness. This conclusion was corroborated by immunotherapy predictions derived from the Easier package (Fig. 3D). According to the Estimate algorithm, patients with low ICDS index exhibited higher immune scores and lower tumor purity scores (Fig. 3E-F), suggesting that low-ICDS index patients have greater immune infiltration, which aligns with the scoring of immune cells and cytokines. Interestingly, the stromal scores displayed no significant differences among the various ICDS index subgroups (Fig. 3G), suggesting that variations between ICDS subgroups were largely linked to immune cell differences. To investigate the functional differences between the high and low-ICDS index groups, we conducted GO functional pathway enrichment analysis on differential genes. The low ICDS index group was mainly enriched in immune-related pathways, while the high ICDS index group was predominantly enriched in pathways associated with stromal growth (Fig. 3H). The gene mutation waterfall plot revealed that the most common mutated genes in both subgroups were APC, TP53, TTN, and KRAS (Fig. 3I), with only TP53 showing a significant difference (Fig. S3A). Meanwhile, genomic changes (FGA), genomic amplifications (FGG), and genomic deletions (FGL) (Fig. S3B-D), along with tumor mutation burden (TMB), single nucleotide polymorphisms (SNPs), and insertion and deletion mutations (INDELs) between the two groups showed no significant differences (Fig. S3E-G). Genomic findings indicated that the immune disparities between the two patient groups did not correlate with tumor mutational burden, suggesting that changes in adjuvant-like factors were the primary drivers of the observed differences in immune infiltration.

Fig. 3

Characteristics and differences between the ICDS subgroups. A-C Differences between the ICDS subgroups in 28 immune and stromal cells (A), immune factors (B), and immune checkpoints (C). D Boxplot showing the differences in immunotherapy response scores (Easier Score) between the ICDS subgroups. E–G Boxplots displaying differences in Immune Score (E), Tumor Purity (F), and Stromal Score (G) between ICDS subgroups. H Gene Ontology (GO) pathway enrichment analysis showing significant differences in differential genes between the ICDS subgroups. I Tumor mutation burden (TMB) distribution and mutation waterfall plot for the ICDS subgroups. J-K Gene set enrichment analysis (GSEA) of cell death (J) and cell stress (K) showing relative enrichment differences between the ICDS subgroups. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001

Further investigation into the differences in cell death and stress between the two groups using GSEA analysis showed that the low-ICDS index group exhibited significant enrichment in multiple cell death pathways, especially cuproptosis (Fig. 3J), along with endoplasmic reticulum stress and oxidative stress (Fig. 3K). In conclusion, these findings validated that ICDS subgroups exhibited marked differences in immune cell infiltration, programmed cell death and cell stress.

Single-cell analysis revealed the key differential cell subpopulations and regulatory genes of the ICDS model

To investigate the critical cellular subpopulations contributing to ICDS subgroup differences, we performed single-cell analysis. After stringent quality control and filtering of the GSE132465 single-cell dataset, we retained cells originating from tumor tissues, resulting in 27,292 single cells. Using unsupervised graph-based clustering on all cells and identifying them based on classical cell markers, we identified six main cell types: B cell, Endothelial cell, Epithelial cell, Fibroblast, T/NK cell, and Myeloid cell (Fig. 4A-B). Then, AUCell was used to calculate the enrichment scores of ICDS subgroup differential genes across each cell type, revealing that myeloid cells were highly enriched with ICDS subgroup differential genes (Fig. 4C-D). To further identify the key subpopulations contributing to the differences, we conducted additional unsupervised graph-based clustering on myeloid cells, classifying them into three cell types (Dendritic cell, Macrophage, Monocyte) using classical markers (Fig. 4E and Fig. S4A). Subsequent AUCell calculations revealed that the differential gene scores were highest in macrophages (Fig. 4F). Macrophages are broadly classified into pro-inflammatory M1-type and anti-inflammatory M2-type based on their functions. We classified the three macrophage clusters (Fig. 4G-H) as M1-like and M2-like types based on their functional roles. Cluster 1 was primarily enriched in functions related to synapse modification, lipoprotein catabolism, collagen degradation, and regulation of lipid transport, which are associated with wound healing and tissue repair (Fig. S4B). Thus, we defined this cluster as M2-like macrophages. Conversely, Cluster 2 (Fig. S4C) and Cluster 3 (Fig. S4D) were enriched in functions related to leukocyte aggregation, neutrophil chemotaxis, inflammatory response and MHC-II complex assembly, associated with pro-inflammatory and antigen presentation roles. Therefore, we defined these cells as M1-like macrophages. Interestingly, we found significant differences between M1-like and M2-like macrophages when we recalculated enrichment scores using differential genes (Fig. 4I). These findings suggest that variations in macrophage types could play a major role in driving the immune differences observed between the ICDS subgroups.

Fig. 4

Single-cell analysis revealed the key differential cell subpopulations and regulatory genes of ICDS. A The t-SNE plot of unsupervised clustering of cells from the single-cell dataset. B Heatmap of characteristic genes for the cell types. The dot size reflects the percentage of cells expressing these genes, and dot color intensity indicates the average expression levels of the genes. C Ridge plot displaying differences in Area Under the Curve (AUCell) scores for differential genes derived from the ICDS subgroup across cell types. D The t-SNE plot reveals differences in the density distribution of enrichment scores for differential genes from the ICDS subgroup across cells. E Unsupervised clustering of myeloid cells, represented as a t-SNE plot. F The box plot illustrates differences in AUCell scores for the ICDS differential genes among three myeloid cell types. G-H Unsupervised clustering of macrophages resulted in three clusters (G), which were further categorized into two functional subtypes (H), as shown in the t-SNE plot. I The box plot illustrates differences in AUCell scores for ICDS differential genes between two macrophage populations. J The heatmap shows the relative expression values of model genes across different cell subpopulations. K-L The t-SNE plots display the differences in density distribution of GAL3ST4 in myeloid cells (K) and macrophages (L). M Gene preferences of different macrophages estimated by the STARTRAC-dist index (+ +  + , Ro/e > 1; +  + , 0.8 < Ro/e ≤ 1; + , 0.2 ≤ Ro/e ≤ 0.8; + / − , 0 < Ro/e < 0.2; − , Ro/e = 0). N qRT-PCR analysis of GAL3ST4 expression levels in distinct polarization states of macrophages. O Western blotting reveals GAL3ST4 activity in distinct polarization states of macrophages. P Display of human colon cancer tissue sections stained with DAPI/CD163/GAL3ST4. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001

To further investigate whether model genes play a driving role in these differences, we examined their distribution across all cell populations. The model genes ATP6V1B2 and GAL3ST4 were predominantly expressed in myeloid cells (Fig. 4J and Fig. S4E-M), with GAL3ST4 specifically expressed in macrophages (Fig. 4K), while ATP6V1B2 showed widespread expression across the three main myeloid cell types (Fig. S4N). Notably, GAL3ST4 demonstrated high expression in M2-like macrophages (Fig. 4L). This observation was further validated by the STARTRAC-dist index, which highlighted the distributional difference (Fig. 4M). Moreover, violin plots revealed that the expression differences of GAL3ST4 between the groups were statistically significant (Fig. S4O). The differential expression of GAL3ST4 in macrophages of varying states suggested that it may play a key regulatory role in the transition between different macrophage polarization states. In the TCGA-CRC cohort, we observed higher levels of GAL3ST4 expression in the group with a poor prognosis and high ICDS index (Fig. 1I and Fig. S4P). To validate the specific expression of GAL3ST4 in M2-like macrophages, we induced THP-1 cells to differentiate into various macrophage states (M0, M1, M2). Our results showed elevated GAL3ST4 expression in M2-like macrophages, demonstrated through both mRNA and protein analyses using qRT-PCR and Western blot assays (Fig. 4N-O). Additionally, immunofluorescence staining on human colon cancer sections showed co-localization of GAL3ST4 with CD163-high expressing macrophages (M2-like) (Fig. 4P and Fig. S5A-B).

In conclusion, our results demonstrated the specific overexpression of GAL3ST4 in M2-like macrophages, with the differing functional states of macrophages being the primary factors contributing to the prognostic differences in the ICDS subgroup. Therefore, the prognostic disparities in ICDS patients might stem from the tumor-promoting functions of GAL3ST4-driven M2-like macrophages.

GAL3ST4 impacts the M2-like phenotype and functional characteristics of macrophages

Due to the notable differences in GAL3ST4 expression observed between M1-like and M2-like macrophages, we hypothesized that this gene could play a crucial role in regulating the switching of macrophage phenotypes. To explore this, we conducted the monocle analysis and identified a pseudotime trajectory from M2-like to M1-like macrophage transformation (Fig. 5A-B). The expression of GAL3ST4 progressively decreased as M2-like macrophages transitioned to M1-like macrophages (Fig. 5C), implying that GAL3ST4 could be pivotal in M2-like macrophage transformation. To validate the bioinformatics results, we established GAL3ST4 knockdown THP-1 cell lines (Fig. 5D). Immunofluorescence staining revealed that CD163 expression was reduced in GAL3ST4 knockdown M2-like macrophages, suggesting that GAL3ST4 indeed affects the M2-like macrophage phenotype (Fig. 5E). Through monocle BEM Heatmap analysis, we observed that pro-tumor immune factors or molecules (CCL8, CCL2, C1QC, A2M) progressively increased with M2-like transition. In contrast, anti-tumor immune factors or molecules (IL1B, TNFRSF11B, CCR7, IL1R2) increased with M1-like transition (Fig. 5F and Fig. S6A-G). qRT-PCR results validated that CCL8, CCL2, C1QC, A2M, GPNMB, LAPTM5, and SERPING1 were significantly downregulated in GAL3ST4 knockdown M2-like macrophages (Fig. 5G and Fig. S6H-I). CCK8 assays demonstrated a marked decrease in the tumor-promoting capacity of M2-like macrophages following GAL3ST4 knockdown (Fig. 5H-I). Co-culturing M2-like macrophages with colon cancer cells in transwell chambers revealed that the ability of M2-like macrophages to promote colon cancer cell invasion was significantly reduced after GAL3ST4 knockdown (Fig. 5J and Fig S8A).

Fig. 5

The impact of GAL3ST4 on the M2-like phenotype and functional properties of macrophage. A-C Time-course analysis of macrophage differentiation, with assessments of macrophage differentiation states (A), pseudotime (B), and GAL3ST4 expression (C). D Constructing GAL3ST4 knockdown M2 macrophages. E Immunofluorescence staining of M2 macrophages, showing CD163 (red) and DAPI (blue), with a merged image of CD163 (red) and DAPI (blue). F Heatmap of the relative expression of representative genes across different differentiation states inferred from trajectory analysis. G qRT-PCR results of A2M, CCL2, CCL8, and C1QC in GAL3ST4 knockdown macrophages and control macrophages. H-I CCK8 viability assay performed at 12, 24, 48, and 72 h after treating DLD1 (H) and HCT116 (I) cells with conditioned medium from different differentiation states of macrophages and GAL3ST4 siRNA-transfected macrophages. J Transwell migration assay of DLD1 and HCT116 cells co-cultured with different differentiation states of macrophages and GAL3ST4 siRNA-transfected macrophages. si#1:si-GAL3ST4#1; si#2:si-GAL3ST4#2; si#3:si-GAL3ST4#3. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001

These findings indicated that GAL3ST4 is essential for sustaining the M2-like phenotype and pro-tumor functions of macrophages.

Tumor cells undergoing cuproptosis influence GAL3ST4 expression in macrophages and inhibit their tumor-promoting functions

According to the previous GSEA enrichment analysis, cuproptosis exhibited the most significant enrichment differences among the ICDS index subgroups (Fig. 6A, Fig. 3I). Prior studies had also suggested a potential link between cuproptosis and macrophage phenotype expression [25]. Correlation analysis revealed a significant negative correlation between high ICDS index subgroup differential gene enrichment scores and cuproptosis enrichment scores (Fig. 6B), while a significant positive correlation was observed in the low ICDS index subgroup (Fig. 6C). Since the differences in ICDS index subgroups were mainly driven by macrophage differences (Fig. 4E-I), we hypothesized that cuproptosis might play a role in regulating macrophage functions. To validate this inference, we used inferCNV to identify malignant epithelial cells (Fig. S7A-B). We then categorized these cells into two groups: high cuproptosis malignant epithelial cells (HC EpiT) and low cuproptosis malignant epithelial cells (LC EpiT), based on their cuproptosis scores. Using CellChat to analyze intercellular communication between HC EpiT and other cells, we found that it had a strong interaction strength with macrophages (Fig. 6D). Further CellChat analysis of the interactions between HC epithelial cells and M1-like and M2-like macrophages revealed that M2-like macrophages exhibited a stronger interaction weight with HC epithelial cells (Fig. 6E), suggesting that cuproptosis may regulate the pro-tumor phenotype of macrophages. Next, we co-cultured Cu-treated CRC cell fragments with M2-like macrophages to investigate the effect of cuproptosis CRC cells on M2-like macrophages. Transwell and CCK8 assays indicated that cuproptosis CRC cells inhibited the tumor-promoting invasion and activity of M2-like macrophages (Fig. 6F-H and Fig. S8B). Immunofluorescence staining showed that CD163 expression in M2-like macrophages was reduced after treatment with cuproptosis CRC cell fragments (Fig. 6I). Additionally, pro-tumor factors, including CCL8, CCL2, C1QC, A2M, GPNMB, LAPTM5, and SERPING1, were significantly downregulated (Fig. 6J and Fig. S7C-E). The results showed that cuproptosis in tumor epithelial cells suppressed the phenotype and pro-tumor functions of M2-like macrophages. Interestingly, we found that GAL3ST4 expression was significantly reduced in M2-like macrophages following treatment with cuproptosis CRC cell fragments (Fig. 6K), aligning with the previous findings on GAL3ST4's role in regulating macrophage function. Moreover, GAL3ST4 showed a significant negative correlation with several key cuproptosis regulatory genes (Fig. S7F).

Fig. 6

Cuproptosis in relation to macrophages and GAL3ST4. A Box plot showing differences in GSVA enrichment scores for cuproptosis in the ICDS subgroups. B-C Correlation analysis between cuproptosis GSVA enrichment scores and the ICDS subgroups characteristic gene GSVA enrichment scores. D The circular plot illustrates the interactions in the cell–cell communication network between high cuproptosis malignant epithelial cells (HC EpiT) and other cells. E Bar chart representing the interaction strength between macrophages and cuproptosis malignant epithelial cells. F-G CCK8 viability assay was performed at 12, 24, 48, and 72 h on DLD1 (E) and HCT116 (G) cells treated with conditioned media from macrophages co-cultured with cuproptosis CRC cell fragments. H Transwell migration assay of DLD1 and HCT116 cells co-cultured with macrophages and macrophages co-cultured with cuproptosis CRC cell fragments. I Immunofluorescence staining of macrophages and macrophages co-cultured with cuproptosis CRC cell fragments, showing CD163 (red) and DAPI (blue), with a merged image of CD163 (red) and DAPI (blue). J qRT-PCR results of CCL8, CCL2, C1QC and A2M in macrophages and macrophages co-cultured with cuproptosis CRC cell fragments. K qRT-PCR results for GAL3ST4 mRNA in macrophages and macrophages co-cultured with cuproptosis CRC cell fragments. CU: cuproptosis CRC cell fragments. *p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001

These results suggested that CRC cells undergoing cuproptosis could influence the expression of the GAL3ST4 gene in macrophages. This interaction subsequently affected the phenotype and pro-tumor functions of M2-like macrophages, ultimately impacting the prognosis of CRC patients.

Discussion

Programmed cell death (PCD) and stress are biological processes that are highly correlated with immune cell infiltration within the tumor microenvironment [20]. By inducing PCD and stress in the immune microenvironment, it is possible to remodel the immune microenvironment and improve immunotherapy [26]. Although this biological phenomenon has been applied in various novel therapeutic strategies with encouraging efficacy [27], the complexity of the tumor immune microenvironment still presents significant challenges in exploring the crosstalk mechanisms between tumor cell death and the tumor immune microenvironment. Unveiling the mechanisms of crosstalk between PCD, stress, and the immune microenvironment is crucial for identifying new therapeutic targets, optimizing treatment strategies, and improving patient prognosis. In this study, we identified genes closely associated with PCD, stress, and immune cells in colorectal cancer (CRC) and constructed the ICDS prognostic model using machine learning as a foundation for further investigation into the crosstalk between PCD, stress, and the immune microenvironment.

Through the evaluation of multiple CRC cohorts, the ICDS prognostic model demonstrated robust predictive performance. Univariate and multivariate Cox regression analyses further indicated that the ICDS index has potential as an independent prognostic indicator. Subsequently, we used the Estimate and ssGSEA algorithms to analyze differences in immune cell infiltration within the tumor microenvironment of ICDS subgroups. Specifically, as the ICDS index increased, the infiltration of various immune cells (such as macrophages and T cells) in the tumor microenvironment (TME) of CRC patients decreased significantly, and patient prognosis gradually worsened. Combined with previous studies, these changes in immune and prognostic characteristics suggest that TME of patients transitions from an inflammatory infiltration state to an immune exhaustion state, indicating the potential value of the ICDS model in predicting immune infiltration status [24, 28, 29].Our results showed that the activation levels of different PCD pathways are closely associated with the ICDS index, with subgroups having lower ICDS indices exhibiting significantly higher activation levels, which are associated with better prognosis. These findings further highlight the potential biological value of the ICDS model. Remarkably, we found that cuproptosis exhibited the most significant differential enrichment scores between high and low ICDS index subgroups, as previously reported, cuproptosis might play a crucial role in the crosstalk between PCD and immune cells in CRC [30, 31]. Although previous studies have created prognostic models that connect cuproptosis and immune-related features to outcomes in CRC patients, these models lack depth. They focus solely on clinical characteristics and fail to explore the underlying cellular interactions that contribute to the observed differences in prognosis [32, 33].

The rapid development of single-cell technologies has provided an essential tool for deciphering interactions between different cell populations within TME. To explore the underlying mechanisms further, we used the AUCell algorithm to map bulk transcriptomic data to specific cell types and study the interactions among PCD, stress responses, and immune cells. Based on single-cell analysis results, we found that differentially expressed genes between high-ICDS index and low-ICDS index subgroups were predominantly enriched in M2-like and M1-like macrophages. Previous studies have shown that macrophages exhibit multiple functional states, where M1-like macrophages promote the infiltration of inflammatory cells in the tumor microenvironment, whereas M2-like macrophages suppress inflammatory responses and facilitate immune exhaustion [34, 35]. Combining our single-cell analysis results, we speculated that the transition in macrophage polarization states is a crucial biological mechanism driving the immune predictive power of the ICDS model. Interestingly, single-cell pseudotime analysis revealed that the evolutionary trajectory of the model gene GAL3ST4 closely aligns with the trajectory of M2-like macrophages. However, few studies have explored the specific regulatory role of GAL3ST4 in the evolution of M2-like macrophages. Through biological experiments, we further validated the regulatory effect of GAL3ST4 on the pro-tumor capabilities of M2-like macrophages. This regulatory mechanism may represent the key biological basis underlying prognostic differences between ICDS index subgroups.

Additionally, in transcriptomic data cohorts, we found that GAL3ST4 is significantly negatively correlated with several key genes related to cuprotosis [36,37,38,39]. This intriguing finding prompted us to further investigate the regulatory relationship between tumor cell cuprotosis and macrophage phenotype evolution. A recent study showed that copper-induced cell death in tumor cells could enhance the efficacy of PD-1 monoclonal antibody immunotherapy and reduce the proportion of M2-like macrophages in TME. However, the specific mechanism underlying the reduction in M2-like macrophages was not clarified in that study [25]. Through further single-cell analysis and biological experiments, we validated that cuprotosis in tumor cells inhibits the expression of GAL3ST4, significantly suppressing the pro-tumor functions of M2-like macrophages. This discovery provides a new perspective on the mechanisms linking cuprotosis and immune microenvironment regulation.

Our study combines machine learning and single-cell sequencing technologies to reveal that cuproptosis in tumor cells suppresses the pro-tumor function of M2 macrophages by downregulating GAL3ST4 expression in macrophages. However, recent studies have shown that the relationship between PCD and TME is widespread across various cancer types. Therefore, it is necessary to further investigate the applicability of this model in other cancers [40, 41]. Meanwhile, the conclusions of this study are primarily based on omics analyses and in vitro experiments, highlighting the need for further in vivo validation of the mechanism. For instance, constructing GAL3ST4-specific knockout mouse models would help elucidate its precise role in macrophage polarization and the overall regulation of the TME. Furthermore, GAL3ST4 could serve as a therapeutic target for drug development, and its efficacy as a standalone therapy or in combination with cuprotosis inducers could be assessed for synergistic effects, followed by validation of its therapeutic potential in animal models. Additionally, the ICDS prognostic model needs optimization to enhance its clinical applicability and practical utility.

Conclusion

We developed a prognostic model focused on immune-related cell death and stress (ICDS) that accurately predicts survival outcomes for colorectal cancer (CRC) patients and highlights immune response variations among distinct prognostic groups. Additionally, our research reveals that tumor cells undergoing cuproptosis can suppress M2-like macrophage polarization by downregulating GAL3ST4. In summary, this study provides a reliable CRC prognostic model, uncovering the crosstalk between the cuproptosis in tumor cells and macrophages, and providing potential targets and biomarkers for CRC treatment and prognosis evaluation.