Introduction

Colorectal cancer (CRC) ranks as the third most common and second most fatal cancer worldwide1. Despite surgery being the primary treatment for metastatic CRC, the cure rate remains limited to 20–30%, with a post-surgery recurrence rate of 50–75%2. Although advances in surgical techniques, chemotherapy, and comprehensive treatment strategies have improved patient survival and quality of life3, the prognosis remains poor due to late diagnosis, rapid progression, and early metastasis4,5. Multidisciplinary approaches have reduced mortality, with annual decreases of 3% in the 2000s and 1.8% from 2012 to 20172. Nevertheless, understanding the tumor characteristics and identifying prognostic markers are crucial for the development of personalized treatments.

Copper (Cu) is an indispensable trace element in eukaryotes, playing key roles in iron transport, detoxification of reactive oxygen species, and mitochondrial respiration6. Dysregulation of its transmembrane transport can lead to intracellular copper accumulation, resulting in cytotoxicity and cell death7. However, the specific molecular mechanisms by which Cu induces cell death remained unclear until 2022, when Tsvetkov et al. introduced a novel pathway termed cuproptosis8. They found that when cells were exposed to 40 nM of the cuproptosis-inducing agent elesclomol, there was a significant increase in intracellular Cu levels, accompanied by marked cell death. This cell death mechanism is distinct from apoptosis9, necrosis10, autophagy11, necroptosis12, pyroptosis13, and ferroptosis14, marking a significant advancement in understanding copper-induced cytotoxicity.

Cuproptosis is a copper-dependent form of cell death. Copper binds to lipoylated components of the tricarboxylic acid (TCA) cycle, leading to an accumulation of these copper-bound mitochondrial proteins and a reduction in iron-sulfur (Fe-S) cluster proteins. The mitochondria are the primary targets, experiencing oxidative damage to their membranes and disrupted enzyme function in the TCA cycle15. Mitochondrial dysfunction may contribute to Cu-induced neurotoxicity, as evidenced by elevated mitochondrial Reactive Oxygen Species (ROS) levels, reduced pyruvate dehydrogenase activity, and impaired respiratory complex I in neuroblastoma cells exposed to Cu16. Additionally, Cu exposure collapses mitochondrial membrane potential (ΔΨm) in neuronal/glial cultures15. Elesclomol (ES), a copper ion carrier, has elucidated the role of Cu in cells. A 2012 study found that ES facilitated Cu transport in melanoma cells, decreasing mitochondrial protein levels, increasing ROS, and suppressing tumor cell proliferation17. Buccarelli et al. reported that elesclomol-induced intracellular Cu increases mitochondrial ROS production and decreases glutathione (GSH) levels18, with GSH depletion increasing cancer cell susceptibility to cuproptosis19. Studies suggest that elesclomol-copper (ES-Cu) complexes may block cells in the G1 phase, induce DNA damage, and affect mitochondrial membrane potential20. Tsvetkov et al. discovered that ferrotoxin 1 (FDX1), crucial for ES sensitivity, binds directly to ES-Cu, inhibiting Fe-S cluster formation21. These findings offer new insights into cuproptosis induced by ES-Cu complexes.

The association between copper (Cu) and colorectal cancer (CRC) has gained increasing recognition. Aubert et al. reported significant differences in copper levels between the serum and cancerous tissues of CRC patients compared to normal subjects22. Similarly, Hendrych et al. demonstrated that Cu ion carriers could enhance the antitumor efficacy of chemotherapeutic agents or overcome chemoresistance in various tumor models by modulating Cu levels23. Furthermore, genes related to cuproptosis, such as FDX1 and DLAT, are predominantly downregulated in CRC tissues24. Based on these findings, we selected colorectal cancer as the focus of our study to explore the relationship between cuproptosis and CRC.

The objective of this research was to comprehensively evaluate the regulation of cuproptosis in CRC. Through WGCNA and differential expression analysis, we identified a gene module consisting of 1,951 genes associated with cuproptosis progression in CRC. Additionally, six novel genes were screened using LASSO regression and random forest methods to establish a prognostic signature for cuproptosis-related mRNAs. Further analysis identified P4HA1 as a key gene modulating the impact of cuproptosis in CRC. This study enhances the accuracy of prognostic assessment and informs the selection of targeted treatments, potentially improving therapeutic outcomes and extending patient survival.

Materials and methods

Collection and preprocessing of data

We obtained 44 normal and 542 CRC samples from the TCGA database (accessed 13 April 2023) and converted the data to FPKM format. Additionally, CRC datasets GSE17538 and GSE29623 were downloaded from GEO (accessed 13 April 2023). The Robust Multichip Average technique was used to standardize and adjust baseline data25. Clinical profiles of all CRC patients are provided in Table 1. We also retrieved the single-cell sequencing (scRNA-seq) CRC dataset GSE231559, consisting of 11 normal and 15 tumor samples (GEO, accessed on 28 March 2024). Seurat objects were created using the “Seurat” R package, excluding cells with fewer than 1000 genes, fewer than 200 detected genes, or mitochondrial gene content exceeding 20%. After merging the datasets, log normalization was applied, and the 2000 most variable genes were identified using the “FindVariableFeatures” function. Dimensionality reduction was achieved using Uniform Manifold Approximation and Projection (UMAP). Differentially expressed genes within cell clusters were identified using the “FindAllMarkers” algorithm with criteria of |log2FC|≥ 0.25 and adjusted p < 0.05. Cell subpopulations were annotated using the “SingleR” software package26.

Table 1 Baseline characteristics of colorectal cancer patients from the TCGA database.
Full size table

Identification of hub cuproptosis-related genes

Using the “ConsensusClusterPlus” R package, we performed unsupervised clustering to classify 542 CRC patients into two groups27. Differentially expressed genes (DEGs) were screened using the “limma” R package with |logFC|> 0.585 and p < 0.0528. A cuproptosis gene set consisting of 12 CRGs was assembled, and cuproptosis scores for each patient were calculated using the “GSVA” R package29. We applied WGCNA to select genes most relevant to the cuproptosis score. Samples and genes with excessive missing values were filtered using the goodSampleGenes function in the “WGCNA” package. Hierarchical clustering was employed to eliminate outliers, and a soft threshold of 4 was selected using the pickSoftThreshold function. The Topological Overlap Matrix (TOM) was computed, and gene modules were identified by dynamic tree cutting with a minimum module size of 100 genes. Modules with similar expression patterns were merged at a threshold of 0.1. The turquoise module showed the highest correlation with clinical traits (correlation coefficient = 0.58, p = 1 × 10-50). We intersected genes from the DEGs and the turquoise module to identify potential cuproptosis-related genes. Prognosis-related genes were first identified using univariate Cox regression analysis. These genes were then further refined through LASSO regression30 and Random Forest analysis31, ultimately selecting six key genes. These six genes were subsequently used to construct the prognostic model through multivariate Cox regression analysis.

Creation of cuproptosis-related prognostic signature

Hazard ratios (HR) for each gene were calculated using the “survival” R package. The TCGA tumor data served as the training set, while GSE17538 and GSE29623 were used as independent validation sets. After identifying six key genes using LASSO regression and Random Forest analysis, a multivariate Cox regression was applied to construct the prognostic model. For each gene in the model, the coefficient (coef) was derived from the multivariate Cox regression analysis, reflecting the strength and direction of the association between gene expression and overall survival. This equation allows us to calculate an individualized risk score for each patient based on the expression levels of the six selected genes and the contribution of each gene to the overall survival, as determined by the Cox regression coefficients:

$$RiskScore={\sum }_{i=1}^{N}({exp}_{i}*{coef}_{i})$$

Patients were categorized into high- and low-risk groups based on the median risk score. Kaplan–Meier survival curves were used to compare overall survival (OS) between the two groups. ROC curves, generated by the “timeROC” R package, assessed model accuracy. Univariate and multivariate Cox regression analyses identified factors predicting survival. A nomogram predicting CRC patient survival was constructed using the “rms” R package.

Gene function enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the “clusterProfiler” and “enrichplot” R packages. A significance threshold of p < 0.05 was set for enrichment results32,33.

Analysis of immune infiltration, somatic mutation and drug sensitivity

Immune cell infiltration and functionality were assessed using single-sample genomic enrichment analysis (ssGSEA)34. The CIBERSORT algorithm estimated the infiltration of 22 immune cell types in CRC samples35. The ESTIMATE algorithm was employed to estimate immune, stromal, and ESTIMATE scores in malignant tissues based on gene expression data.

Somatic mutation data were obtained from TCGA, and mutation burden was compared between low- and high-risk groups using the “maftools” R package36.

The “pRRophetic” R package was used to predict the IC50 values of chemotherapy drugs, where IC50 represents the concentration at which 50% inhibition of biological activity is achieved37.

Analysis of cellular communication and receptor-ligand interaction

The R package “CellChat” was used to analyze CRC scRNA-seq data and identify common CellChat objects. The “mergeCellChat” function consolidated these objects using “CellChatDB.Human” as the reference. The “netVisual_diffInteraction” function was then applied to distinguish ligand-receptor interactions and visualize gene distribution.

Pseudo-temporal analysis of cell trajectories

Pseudotime analysis was conducted using the “Monocle 2” R package. Preprocessed scRNA-seq data were used to construct cell trajectories using the “DDRTree” algorithm, and pseudotime values were assigned to visualize cell state transitions and gene expression patterns across the trajectory.

Tissue specimens, RT‑qPCR, and western blot (WB)

Fresh colorectal cancer (CRC) tissues and adjacent normal tissues were collected from 20 patients at Jiangxi Cancer Hospital. The study was approved by the Jiangxi Cancer Hospital Ethics Committee (Approval No. 2022KY293). All procedures involving human participants were conducted in accordance with the relevant guidelines and regulations, including the Declaration of Helsinki. Informed consent was obtained from all participants or their legal guardians prior to sample collection.

RNA was extracted from these tissues and from CRC cell lines, followed by reverse transcription using the PrimeScript RT kit. Quantitative PCR (qPCR) was conducted to assess mRNA levels, with GAPDH as the normalization control. The sequences of the primers used are listed in Supplementary Table S1. Relative expression levels of mRNA were calculated using the 2-ΔΔCT method.

After SDS-PAGE separation, proteins (20 µg) were transferred to PVDF membranes. Antibodies were incubated with the membranes. Primary antibodies are as follows: rabbit P4HA1 (1:2000), rabbit HSP70 (1:2000), rabbit FDX1 (1:1000), and rabbit β-actin (1:2000). A secondary antibody (1:5000) was then incubated for 2 h on each sample.

Cell culture and siRNA treatment

Human normal colorectal epithelial cell line FHC and colorectal cancer cell lines SW620, SW480, HT8, HT29, HCT116, and DLD-1 were obtained from the Chinese Academy of Sciences (CAS) Stem Cell Bank (Shanghai, China). HCT116 cells were cultured in 5A medium (FuHeng, China) supplemented with 10% fetal bovine serum (FBS) at 37 °C in a humidified atmosphere containing 5% CO2. SW480 cells were maintained in L15 medium (Gibco, USA) with 10% FBS at 37°℃ in an atmosphere without CO2. The other cell lines were cultured in DMEM (Gibco, USA) with 10% FBS under the same conditions as HCT116. Transfections were performed using Lipofectamine RNAiMAX (Invitrogen, USA) and siRNA (RIBOBIO, China), which were diluted in reduced serum medium, incubated for 20 min, and added to the cell cultures for 48 h.

Detection of the malignant proliferation and invasion capabilities of CRC cells

Cell proliferation was assessed using the MTS assay (Promega, USA) following the manufacturer’s instructions. After siRNA transfection, HCT116 and SW480 colorectal cancer cells were seeded at a density of 5 × 103 cells/well in 96-well plates. Cell proliferation was measured at 0, 24, 48, 72, and 96 h after transfection, with OD490 recorded using a BIO-RAD microplate reader.

For the colony formation assay, cells were seeded at a density of 1000 cells/well into 6-well plates. After 14 days of incubation, cells were fixed with 4% paraformaldehyde for 20 min, stained with 0.1% Crystal Violet for 15 min, and then washed with PBS. Colonies were imaged under a microscope, and colony formation was quantified. In the wound healing assay, HCT116 and SW480 cells transfected with siP4HA1 were seeded in 6-well plates at a density that allowed them to reach 80% confluence. Standardized wounds were created using a 200 µL pipette tip. After washing with PBS to remove detached cells, the medium was replaced with FBS-free medium, and images were captured at 0 and 48 h using a 10 × microscope.

For the invasion assay, transfected HCT116 or SW620 cells were seeded in the upper chamber of a 24-well Transwell insert coated with Matrigel at a density of 1 × 105 cells/well in FBS-free medium. The lower chamber contained medium with 10% FBS as a chemoattractant. After 48 h of incubation, non-invading cells were removed, and invading cells on the lower surface of the membrane were fixed with 4% paraformaldehyde, stained with 0.1% Crystal Violet, and photographed in at least three randomly selected fields of view using an inverted microscope.

Assessment of copper levels

CRC cells were seeded in 6 cm plates at a density of 1 × 10^6 cells/plate and treated with 40 nM Elesclomol (ES) and 2 µM CuCl2 pulses for 2 h. The medium was then changed to remove excess Cu, and cells were incubated for an additional 24 h. After incubation, cells were harvested, resuspended in PBS, and lysed according to the manufacturer’s instructions (E-BC-K775-M, Elabscience, Wuhan, China). Intracellular copper (Cu) levels were measured using a colorimetric assay, and results were normalized to total protein concentration.

GSH assay

To measure intracellular glutathione (GSH) and oxidized glutathione (GSSG) levels, CRC cells were lysed and analyzed using the GSH/GSSG Assay Kit (Beyotime, S0053), following the manufacturer’s protocol. Cell lysates were incubated with a GSH detection reagent, and NADPH solution was added to initiate the reaction. OD values at 412 nm were recorded using a microplate reader, and the relative concentrations of GSH and GSSG were calculated according to the kit’s instructions. Cells were seeded at a density of 2 × 10^6 cells/well in 6-well plates prior to assay preparation.

Measurement of intracellular ROS and mitochondrial membrane potentials

HCT116 and SW480 cells were seeded on glass coverslips in a 12-well plate at 1 × 105 cells/well. Cells were incubated with 10 µM Dihydroethidium (DHE) for 30 min at 37 °C, followed by PBS washes. To stain the nuclei, cells were treated with DAPI (1 µg/mL) for 5 min at room temperature. After a final PBS wash, the coverslips were mounted and the fluorescence signals for ROS (DHE, red) and nuclei (DAPI, blue) were captured using an inverted fluorescence microscope.

Cells seeded on glass coverslips in 12-well plates were incubated with 10 µM JC-1 dye for 20 min at 37°℃. After PBS washes, DAPI (1 µg/mL) was applied for nuclear staining. Fluorescence from JC-1 (red for healthy mitochondria, green for depolarized mitochondria) and DAPI (blue for nuclei) was captured using an inverted fluorescence microscope.

Statistical analysis

Bioinformatic analysis was conducted using R software. Experimental data were analyzed with GraphPad Prism 9.0, with results presented as mean ± standard deviation (SD). All experiments were performed in triplicate. Data comparisons among groups were initially made using one-way analysis of variance (ANOVA) for single-factor experiments. For experiments involving two factors, two-way analysis of variance (ANOVA) was employed to assess the individual and interactive effects of both factors. Post hoc comparisons among groups were performed using Tukey’s Honestly Significant Difference (HSD) test. Comparisons between two groups were performed using the t-test. A p-value of less than 0.05 was considered statistically significant (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001).

Results

Expression and mutational landscape of 12 cuproptosis-related genes

The research flowchart is shown in Fig. 1. We analyzed 12 cuproptosis-related genes (CRGs) based on previous studies8,38. A total of 586 CRC samples (44 normal and 542 tumor) were obtained from the TCGA database to compare mRNA levels between normal and tumor tissues. The results showed that FDX1, UBE2D3, DLD, and PDE3B were downregulated in tumors, while UBE2D1, GLS, MAP2K1, and UBE2D2 were upregulated (Fig. S1A). To further investigate, we retrieved somatic mutation data for all CRC samples from the TCGA database. Chromosomal locations of copy number variation (CNV) mutations in all CRGs are presented in a circle diagram (Fig. S1B). A summary of somatic mutations in the 12 genes revealed that 70 samples (12.92%) exhibited mutations. ATP7A had the highest mutation frequency, followed by PDE3B and MAP2K1, while UBE2D1, UBE2D2, UBE2D3, and FDX1 showed no mutations (Fig. S1C). Most genes exhibited copy number deletions, with only MAP2K1 showing amplification (Fig. S1D). Kaplan–Meier survival analysis was used to evaluate the prognostic value of these genes (Fig. S1E). These findings suggest differential expression and prognostic relevance of certain CRGs in CRC, underscoring significant genetic heterogeneity.

Fig. 1
figure 1

Flow chart of this study.

Full size image

Identification of CRC Molecular Clusters and Gene Collection Based on CRGs

To explore the role of CRGs in CRC, we first conducted unsupervised clustering to identify molecular subtypes of CRC based on the expression profiles of 12 CRGs. The optimal number of clusters was determined at k = 2, categorizing the patients into Cluster A (284 cases) and Cluster B (258 cases) (Fig. 2A). Principal Component Analysis (PCA) confirmed the clear separation of these two clusters based on CRG expression (Fig. 2B). Importantly, patients in Cluster A exhibited better overall survival (OS) compared to those in Cluster B (Fig. 2C), suggesting that CRG expression is associated with patient prognosis. To further investigate the molecular differences between these two clusters, we performed a differential gene expression analysis using the limma R package, identifying 2,292 DEGs between Cluster A and Cluster B (logFC > 0.585, p < 0.05), as visualized in the volcano plot (Fig. 2D). These DEGs provided a preliminary list of genes potentially associated with the biological processes underlying CRC progression and prognosis.

Fig. 2
figure 2

Identification of cuproptosis modification patterns in colorectal cancer using consensus clustering and WGCNA. (A) Two clusters were identified using consensus clustering analysis (k = 2). (B) Principal component analysis demonstrates distinct distributions between the two subtypes. (C) Kaplan–Meier survival curves of patients exhibiting different cuproptosis modification patterns. (D) Volcano plot showing DEGs between Cluster A and Cluster B groups. (E) A dynamic tree-pruning was used to combine genes with comparable cuproptosis-score patterns into a single module, resulting in a hierarchically clustered tree. (F) Correlation coefficients and corresponding p-values between the cuproptosis score and each gene module are displayed in the boxes. (G-H) GO analysis and KEGG analysis of 1,951 genes strongly associated with cuproptosis.

Full size image

In parallel, we calculated cuproptosis scores for each patient sample using GSVA, which quantified the activity of cuproptosis-related pathways based on the expression of the 12 CRGs. Next, we applied WGCNA to identify gene modules associated with these cuproptosis scores. WGCNA revealed 15 mRNA co-expression modules, achieving a scale-free topology model fit exceeding 0.85 and a soft threshold power of 4 (β = 4). Among these, the turquoise module showed the highest correlation with cuproptosis scores (R = 0.84) and contained 7,977 mRNAs, making it the most relevant module for further analysis (Fig. 2E,F).

To refine our gene set and establish a more focused list for prognostic modeling, we took the intersection of the DEGs identified from the unsupervised clustering analysis and the genes within the turquoise module from WGCNA. This intersection yielded 1,951 genes that were strongly associated with cuproptosis progression in CRC (Supplementary Table S2). Next, these 1,951 genes were subjected to functional enrichment analyses. KEGG pathway enrichment analysis revealed their involvement in various cancer-related pathways, including Adherens junctions, p53 signaling, ErbB signaling, Hippo signaling, Sphingolipid signaling, MAPK signaling, and TNF signaling (Fig. 2G). Gene Ontology (GO) enrichment analysis further indicated that these genes are involved in critical biological processes (BP), cellular components (CC), and molecular functions (MF) relevant to CRC progression (Fig. 2H).

Construction and validation of the prognostic cuproptosis-related gene signature

To construct a robust prognostic model based on cuproptosis-related genes, we initially performed univariate Cox regression analysis on the 1951 cuproptosis-related genes identified in CRC patients from the TCGA database. This analysis identified 46 prognostic mRNAs significantly associated with overall survival. To further refine these prognostic markers, we employed LASSO regression and random forest machine learning techniques to narrow down the most relevant genes. This process resulted in the identification of six key hub genes—OSBPL1A, PSMD12, ADAM9, P4HA1, SPINK4, and USP53—which were determined to have the strongest association with cuproptosis in CRC (Fig. 3A–C). The correlation between these six hub genes and cuproptosis scores was further examined through scatter plots (Fig. S2), revealing significant correlations for five of the genes, with the exception of SPINK4. The univariate and multivariate Cox regression results for these genes are displayed in Fig. 3D,E, underscoring their prognostic value. Next, we developed a prognostic cuproptosis-related gene signature based on the multivariate Cox regression analysis of these six genes. This gene signature was expressed as a risk score formula:

Fig. 3
figure 3

Construction and validation of the prognostic cuproptosis-related gene signature. (A–C) Screening for cuproptosis-related genes using LASSO and RF. (D) Univariate and (E) multivariate Cox regression analysis identifies the hub cuproptosis-related genes. (F–H) Distribution of risk scores and OS rates in the training and testing cohorts. (I–K) Kaplan–Meier curves for survival status and survival time in the training and testing cohorts. (L–N) The ROC curve shows the potential of the prognostic cuproptosis-related gene signature in predicting 1-, 3-, and 5-year OS in the training and testing cohorts.

Full size image

Cuproptosis-related risk score = 0.3167*ExprOSBPL1A-0.9077*ExprPSMD12-0.3519*ExprIADAM9 + 0.7178*ExprP4HA1-0.1074*ExprSPINK4-0.2691*ExprUSP53.

Using this risk score formula, we stratified patients into high-risk and low-risk groups. The predictive performance of this signature was validated using the GSE17538 and GSE29623 datasets as independent training cohorts. In both cohorts, patients in the low-risk group exhibited significantly improved overall survival (OS) compared to those in the high-risk group (Fig. 3F–N), confirming the model’s predictive capability. The Kaplan–Meier survival curves further demonstrated that low-risk patients had consistently better OS rates than those in the high-risk group. The prognostic accuracy of the model was also validated using ROC curve analysis, where the area under the curve (AUC) values for 1-year, 3-year, and 5-year OS were above 0.7 in the training cohorts, and above 0.65 in the validation cohorts (Fig. 3L–N), indicating that the model had strong predictive power and robustness.

Association between prognostic cuproptosis-related gene signature and clinical features in colorectal cancer patients

Our findings indicate differences in the clinicopathological characteristics of CRC patients between high- and low-risk groups. Sankey diagrams illustrated correlations among clinical characteristics, risk groups, and patient survival, showing a higher mortality risk in the high-risk group (Fig. 4A). The box plot demonstrated that patients with advanced CRC had significantly higher risk scores compared to those in earlier stages (Fig. 4B). A heatmap displayed the distribution of the six genes in relation to clinical characteristics (Fig. 4C). Furthermore, univariate analysis identified this signature as a significant independent prognostic factor (HR = 1.541; p < 0.001), which was confirmed by multivariate Cox regression analysis (HR = 1.387; p < 0.001) (Fig. 4D,E). To enhance clinical prognosis predictions, we developed a comprehensive nomogram incorporating risk scores and clinicopathological characteristics (T-stage, risk score, age, and stage) (Fig. 4F). This nomogram proved to be a robust tool for predicting 1-, 3-, and 5-year survival, showing the highest AUC for survival prediction (Fig. 4G). In conclusion, our nomogram serves as a valuable tool for predicting the survival prognosis of CRC patients.

Fig. 4
figure 4

Clinical characteristics analysis of the prognostic cuproptosis-related gene Signature. (A) Sankey diagram representing the detailed connection between age, gender, CRGclutser, risk and the live status of patients with CRC. (B) Differences in risk scores in patients with different staging. (C) A heatmap of the distribution of 7 different clinicopathological characteristics of the risk groups of each patient based on the signature. (D-E) Results of the univariate and multivariate Cox regression analyses with regard to the OS of the prognostic cuproptosis-related gene Signature. (F) A nomogram that predicts overall survival at 1-, 3-, and 5 years in patients with colorectal cancer. (G) ROC curves for different clinical characteristics.

Full size image

Immune infiltration analysis

To explore immune cell differences between high-risk and low-risk CRC patient groups, we performed ssGSEA using TCGA data. The analysis revealed a notable enrichment of 15 out of 23 immune-related cell types in the high-risk group (Fig. 5A). Immune function scores were significantly higher in the high-risk group for all functions except MHC class I and Type II IFN Response (Fig. 5B). Further analysis using the CIBERSORT algorithm evaluated 22 tumor-infiltrating immune subpopulations in CRC and their correlation with the cuproptosis risk score. Pearson correlation analysis identified four immune cell types significantly associated with the cuproptosis risk score: Naïve B cells (R = -0.14, p = 0.025) and CD4 memory resting T cells (R = -0.28, p < 0.001) were negatively correlated, while CD8 T cells (R = 0.25, p < 0.001) and regulatory T cells (Tregs) (R = 0.19, p = 0.0028) were positively correlated (Fig. 5C). The six pivotal genes were differentially correlated with various immune cells (Fig. 5D). Additionally, the high-risk group exhibited significantly elevated immune, stromal, and ESTIMATE scores (p < 0.001) (Fig. 5E).

Fig. 5
figure 5

Variation in the immune environment between two different risk groups. (A) ssGSEA analysis of diverse risk groups for immune infiltrating cells and (B) Immune-related function. ssGSEA, single-sample gene set enrichment analysis. (C) Correlation analysis between the cuproptosis risk score and immune infiltrating cells. (D) Correlation between six hub cuproptosis-related genes and immune cells. (E) ESTIMATE analysis of diverse risk groups.

Full size image

TMB and drug susceptibility analysis

This analysis evaluated specific gene mutations in CRC and identified the twenty most significant driver genes. The three genes with the highest mutation frequencies in both risk groups were TP53, TTN, and MUC16. The high cuproptosis risk score group (93.57%) exhibited a broader spectrum of somatic mutations compared to the low-risk group (90.59%) (Fig. 6A,B). Nine drugs—Elesclomol39, Imatinib40, Nilotinib41, Roscovitine42, Shikonin43, Temsirolimus44, Nutlin-3a45, EHT.186446, and LFM-A1347—were selected based on their reported efficacy against colorectal cancer. Drug efficacy was compared between high- and low-risk groups, revealing that the high-risk group showed notable sensitivity to all nine chemotherapeutic agents (Fig. 6C–K). Notably, patients in the high-risk group demonstrated heightened sensitivity to the cuproptosis inducer Elesclomol.

Fig. 6
figure 6

Differences in TMB and drug susceptibility between the distinct risk group. (A-B) Top 20 genes with the highest mutation frequencies in the high-risk group and low-risk group. (C) Difference in sensitivity to Elesclomol between the two risk groups. (D) Imatinib. (E) Nilotinib (F) Roscovitine (G) Shikonin (H) Temsirolimus (I) Nutlin-3a (J) EHT.1864 (L) LFM-A13.

Full size image

Expression analysis of six model genes in scRNA-seq and CRC samples

The GSE231559 single-cell RNA sequencing dataset (11 normal and 15 tumor samples) was utilized to investigate potential cuproptosis-related genes. After gene filtering, normalization, and PCA, 13 distinct cell clusters were identified (Fig. 7A). Cells were categorized into seven types based on characteristic markers and genetic profiles: T cells, NK cells, monocytes, epithelial cells, B cells, endothelial cells, and smooth muscle cells (Fig. 7B). Violin plots showed that OSBPL1A, ADAM9, P4HA1, and USP53 were highly expressed in smooth muscle cells, while PSMD12 was expressed in smooth muscle cells, endothelial cells, and monocytes. SPINK4 was primarily expressed in epithelial cells (Fig. 7C). The epithelial cell fraction was isolated from the scRNA-seq data to compare gene expression levels between non-malignant and malignant epithelial cells. OSBPL1A was markedly downregulated in tumor epithelial cells, while P4HA1 was upregulated. No significant differences were observed for PSMD12, ADAM9, SPINK4, and USP53 between normal and tumor cells (Fig. 7D). Analysis of CRC patient tissue samples revealed elevated expression of ADAM9 (p < 0.05) and P4HA1 (p < 0.01) in tumor tissues. In contrast, OSBPL1A (p < 0.001) and USP53 (p < 0.01) were reduced in adjacent normal tissues, consistent with scRNA-seq data, while PSMD12 and SPINK4 showed no significant changes (Fig. 7E).

Fig. 7
figure 7

Distribution of expression in single-cell profiles and sample tissues. (A) UMAP plot of cell subtypes. (B) UMAP plot of the annotation results of CRC cell subgroups. (C) The violin plots show the expression levels of model-built genes in all cell types. (D) Expression of OSBPL1A, PSMD12, ADAM9, P4HA1, SPINK4, and USP53 in Normal and Tumor Epithelial Cells from Single-Cell Transcriptome Data. (E) Expression of OSBPL1A, PSMD12, ADAM9, P4HA1, SPINK4 and USP53 in CRC and adjacent normal tissues.

Full size image

Given its highest hazard ratio (HR = 2.050, p < 0.001) among the identified genes in multivariate analysis (Fig. 3E), as well as its strong correlation with the cuproptosis score (R = 0.41, p < 0.001, Fig. S2B), P4HA1 emerged as a critical gene of interest. Furthermore, single-cell RNA sequencing and PCR analyses confirmed significant differential expression of P4HA1 between tumor and normal tissues, with the smallest p-value among the six hub genes (Fig. 7D,E). These findings highlight P4HA1 as a potential high-risk prognostic factor and justify its selection for further in-depth investigation into its role in colorectal cancer progression and its involvement in the cuproptosis pathway.

Single-cell transcriptomic analysis of P4HA1

Epithelial cells from the GSE231559 tumor samples were categorized as malignant, while those from normal samples were considered non-malignant. Pseudotemporal analysis revealed that as non-malignant cells transitioned to malignant cells, the percentage of cells with high P4HA1 expression increased (Fig. 8A). Using the “CellChat” R package, we analyzed variations in cellular communication between high and low P4HA1 expression groups in tumor samples. The high P4HA1 expression group exhibited a greater number of interactions compared to the low expression group (Fig. 8B). Notably, there was an increase in outgoing and incoming interaction strength between smooth muscle myocytes and endothelial cells in the high expression group (Fig. 8C). Signaling pathways such as CDH1, DESMOSOME, HGF, and NRG were significantly elevated in the epithelial cells of the high expression group (Fig. 8D).

Fig. 8
figure 8

Overall analysis of P4HA1 in single-cell datasets. (A) Pseudotime plots exploring changes in gene expression with progression of non-malignant to malignant cells (B) The network diagram shows the number of interactions between cell types in low and high expression groups. (C) Bubble plots illustrate the incoming and outgoing cell interaction strength in the low expressing patients and high expressing patients. (D) Heatmap of the overall signaling pathway in patients with high and low expression. (E) KEGG enrichment analysis results. (F) GO enrichment analysis results. (G) Scatter plot of correlation between P4HA1 expression level and cuproptosis score.

Full size image

Differentially expressed genes (logFC > 0.25, p < 0.05) (Supplementary Table S3) between the high and low P4HA1 expression groups were identified and subjected to KEGG and GO enrichment analyses. KEGG analysis suggested that these genes may be involved in pathways associated with neurological diseases (Parkinson’s, Huntington’s, Alzheimer’s), ribosome function, oxidative phosphorylation, chemical carcinogenesis (reactive oxygen species), and cell adhesion molecules, among others (Fig. 8E). GO enrichment analysis indicated associations with cytoplasmic translation, oxidative phosphorylation, aerobic respiration, electron transport chain, and cellular respiration (BP); ribosomal subunit, cytosolic ribosome, mitochondrial protein-containing complex, and cell-substrate junction (CC); as well as ribosome structure, cadherin binding, electron transfer activity, enzyme inhibitor activity, and cell adhesion mediator activity (MF) (Fig. 8F).

Furthermore, a significant correlation was observed between P4HA1 expression in tumor epithelial cells and the cuproptosis score derived from 12 cuproptosis-related genes (Fig. 8G). These findings suggest that P4HA1 may play a role in tumorigenesis and progression, potentially influencing tumor cells through mechanisms involving cuproptosis.

The relationship between P4HA1 expression and tumor cell phenotypes in colorectal cancer

Initially, we extracted eight pairs of colorectal tissues from the sample bank to further validate the expression of P4HA1 in CRC using Western blot analysis (Fig. 9A). P4HA1 expression was significantly higher in various CRC cell lines compared to fetal human colorectal cells, particularly in HCT116 and SW620 cells (Fig. 9B). To elucidate the role of P4HA1 in CRC cells, we used two different siRNAs to knock down P4HA1 expression in HCT116 and SW620 cells. The efficacy of the knockdown was confirmed through Western blot and quantitative PCR assays (Fig. 9C, Fig. S3). Subsequent experiments were conducted to assess whether P4HA1 knockdown affected malignant behavior in CRC cells. Growth curve analysis showed that cells with reduced P4HA1 expression exhibited significantly slower growth rates compared to control cells (Fig. 9D). Migration and invasion abilities were measured by changes in gap distance and the number of migrated cells, revealing that P4HA1 knockdown effectively suppressed these processes (Fig. 9E,F). Consistent with these findings, the colony formation assay demonstrated that P4HA1 knockdown significantly inhibited the colony-forming ability of the cells (Fig. 9G). These results suggest that P4HA1 plays a crucial role in CRC cell proliferation, migration, and invasion.

Fig. 9
figure 9

Validation of P4HA1 affecting colorectal tumor phenotypes. (A) P4HA1 protein expression level in CRC tissues and adjacent normal tissues. (B) The relative expression of P4HA1 mRNA in CRC cells and normal colon epithelial cells. (C) The expression of P4HA1 protein of HCT116 and SW620 cells after treating with two siRNAs. (D) Growth curves of cells without knockout and cells with knockout of P4HA1. (E–G) Wound healing, colony formation, and transepithelial migration assays in CRC cells with and without P4HA1 knockout.

Full size image

Elesclomol induces cuproptosis and correlates with oxidative stress

Previous research has shown that Elesclomol (ES), as a Cu ionophore, can induce cuproptosis in cancer cells8. To investigate this in CRC cells, we treated the cells with ES and CuCl2 (Cu) for 2 h, followed by a medium change and continued culture for 24 h, as described in earlier studies.

Our results indicate that ES + Cu significantly inhibits CRC cell viability, an effect that is reversed by the Cu ion chelator Tetrathiomolybdate (TTM) (Fig. 10A). To confirm that this effect is linked to cuproptosis, we measured intracellular Cu ion content and found that ES + Cu markedly increases Cu ion concentration within tumor cells, while TTM significantly reduces free Cu ions (Fig. 10B). This suggests that the reduced cell viability induced by ES + Cu is associated with elevated Cu ion levels, leading to cuproptosis. Given the established connection between cuproptosis and oxidative stress, we assessed the GSH/GSSG ratio, a key indicator of intracellular redox balance, in CRC cells treated with ES + Cu. The treatment significantly decreased the GSH/GSSG ratio, an effect reversed by TTM (Fig. 10C). Additionally, we evaluated intracellular ROS levels and mitochondrial membrane potential. ES + Cu treatment significantly increased ROS production and decreased mitochondrial membrane potential, effects that were mitigated by TTM (Fig. 10D–G).

Fig. 10
figure 10

Elesclomol induces cuproptosis and is associated with oxidative stress in CRC. (A) Cell viability of HCT116 and SW6620 cells following treatment with elesclomol (40 nM), CuCl2(2 µM), or TTM (20 µM). (B) Measurement of intracellular copper ion concentration in HCT116 and SW620 cells after treatment with elesclomol (40 nM), CuCl2(2 µM), or TTM (20 µM). (C) GSH/GSSG ratio of CRC cells after treatment with elesclomol (40 nM), CuCl2(2 µM), or TTM (20 µM). (D, E) Assessment of mitochondrial membrane potential in CRC cells using JC-1 staining after treatment with elesclomol (40 nM), CuCl2(2 µM), or TTM (20 µM). (F, G) Detection of ROS in CRC cells utilizing the fluorescent probe DHE following treatment with elesclomol (40 nM), CuCl2(2 µM), or TTM (20 µM).

Full size image

In conclusion, our findings suggest that ES + Cu induces cuproptosis in CRC cells, potentially promoting intracellular oxidative stress through elevated Cu ion levels.

Cuproptosis regulation by P4HA1 in CRC

Given the high molecular weight of P4HA1 and its association with cuproptosis, we explored its potential role in regulating this process. The MTS assay revealed a significant decrease in cell viability in both colorectal cancer cell lines with P4HA1 knockdown after treatment with ES and Cu (Fig. 11A). Furthermore, knockdown of P4HA1 further enhanced the increase in intracellular Cu levels induced by ES and Cu (Fig. 11B). Since cuproptosis is characterized by the loss of Fe-S cluster-containing proteins and proteotoxic stress, we examined the levels of FDX1 and HSP70. Western blot results demonstrated that cuproptosis activity was significantly increased in the si-P4HA1 group following a 2-h pulse treatment with ES and Cu, suggesting that P4HA1 may promote cuproptosis in CRC cells (Fig. 11C,D).

Fig. 11
figure 11

Knockdown of the P4HA1 gene potentially enhances cuproptosis triggered by Elesclomol + CuCl2. (A) Cell viability of HCT116 and SW620 cells after si-P4HA1 transfection. (B) Detection of intracellular copper ion concentration in HCT116 and SW620 cells after si-P4HA1 transfection. (C) The expression levels of P4HA1, HSP70, and FDX1 proteins in HCT116 cells following si-P4HA1 transfection (D) The expression levels of P4HA1, HSP70, and FDX1 proteins in SW620 cells following si-P4HA1 transfection.

Full size image

To further investigate the relationship between P4HA1 and oxidative stress induced by cuproptosis, we found that concurrent induction of cuproptosis and P4HA1 knockdown significantly decreased the cellular GSH/GSSG ratio (Fig. 12A), increased ROS generation (Fig. 12C,D), and led to mitochondrial membrane potential collapse (Fig. 12B,E). Notably, these adverse effects were reversed by TTM treatment. Thus, our results suggest that P4HA1 is a critical regulator of the cuproptosis pathway, and its downregulation enhances cellular sensitivity to cuproptosis, likely through the amplification of cuproptosis-induced oxidative stress. These findings identify P4HA1 as a potential therapeutic target in CRC treatment.

Fig. 12
figure 12

Knockdown of the P4HA1 gene exacerbates oxidative stress induced by Elesclomol + CuCl2, which can be reversed by the copper ion chelator TTM. (A) GSH/GSSG ratio in CRC cells after P4HA1 knockdown. (B, E) Mitochondrial membrane potential in HCT116 cells after P4HA1 knockdown. (C, D) ROS levels in HCT116 cells following P4HA1 knockdown.

Full size image

Discussion

In this study, we first identified 12 key marker genes—FDX1, UBE2D1, UBE2D3, DLD, DLAT, GLS, MAP2K1, ATP7A, PDK1, PDE3B, CCS, and UBE2D2—that are strongly associated with cuproptosis, based on findings from previous literature. Using the “ConsensusClusterPlus” R package, we performed unsupervised clustering of colorectal cancer (CRC) samples from the TCGA database, stratifying them into two distinct subgroups based on the expression profiles of these genes. These two subgroups showed significant differences in survival outcomes, suggesting that cuproptosis-related gene expression may have prognostic value. We then performed differential expression analysis, identifying 2292 differentially expressed genes (DEGs) between the two subgroups. To further explore the role of cuproptosis in CRC, we calculated cuproptosis scores for each sample using Gene Set Variation Analysis (GSVA), quantifying cuproptosis pathway activity across the samples. Concurrently, we applied Weighted Gene Co-Expression Network Analysis (WGCNA) to identify gene modules associated with these cuproptosis scores, resulting in the identification of 15 gene co-expression modules, of which the turquoise module showed the strongest correlation with cuproptosis. By taking the intersection of the 2292 DEGs and the genes from the turquoise module, we identified 1951 genes that were both differentially expressed and closely related to cuproptosis progression in CRC. Functional enrichment analyses using GO and KEGG revealed that these 1951 genes are involved in several critical pathways related to copper metabolism, including the EGFR48, PI3K-AKT49 and p53 signaling pathways50, which are known to influence tumor progression. Our GO enrichment analysis further revealed significant correlations with ATP and GTP depletion, consistent with the fact that both copper ion uptake and efflux require energy51 and that cuproptosis involves mitochondrial respiration8. These results suggest that cuproptosis is regulated by copper-dependent mechanisms that influence key cancer-related processes, highlighting the pivotal role of copper in regulating CRC progression through these molecular pathways.

We then applied univariate Cox regression to identify genes significantly associated with patient prognosis, followed by LASSO regression and random forest analysis to refine the selection. This process led to the identification of six key cuproptosis-related hub genes—OSBPL1A, PSMD12, ADAM9, P4HA1, SPINK4, and USP53. These six genes were then incorporated into a multivariate Cox regression model, enabling the stratification of patients into high-risk and low-risk groups with distinct survival outcomes. This approach highlights the potential of these cuproptosis-related genes as reliable prognostic biomarkers for CRC. All of these genes have been previously reported to be closely associated with tumors. OSBPL1A is recognized as an immune- and cancer-associated fibroblast (CAF)-related gene in colorectal cancer52, SPINK4 inhibits ferroptosis and promotes tumor growth53, P4HA1 protects nasopharyngeal carcinoma cells from erastin-induced ferroptosis by activating HMGCS154, PSMD12 may activate the MEK-ERK pathway via KIF15 upregulation, thereby promoting tumor progression55, ADAM9 is highly expressed in various cancers, including non-small cell lung, pancreatic, gastric, breast, ovarian, and colorectal cancers56, and USP53 induces apoptosis in hepatocellular carcinoma cells through the deubiquitination of cytochrome c57. These findings indicate that the identified model genes play significant roles in the progression of various cancers. However, there have been no reports linking these genes specifically to cuproptosis. Thus, we used these six genes to build a prognostic model associated with cuproptosis, stratifying patients into high- and low-risk groups based on their risk scores. It was observed that the prognosis for the high-risk group was considerably poorer compared to the low-risk group. We then integrated the model with clinical characteristics to construct a nomogram for predicting patient survival stages. Preliminary results suggest that the nomogram has significant predictive value, as evidenced by ROC analysis.

The composition and degree of immune cell infiltration in the tumor microenvironment are key factors influencing tumor development and prognosis. Regulatory T cells (Tregs) tend to promote tumor evasion from immune surveillance58. In contrast, CD4 memory T cells enhance macrophage and natural killer cell activity by secreting pro-inflammatory cytokines such as IFN-γ and TNF-α, thereby increasing tumor cell clearance59. In this study, patients with higher risk scores had elevated levels of Tregs, lower levels of CD4 memory T cells, and higher StromalScore and ImmuneScore, indicating a more complex immune environment in the high-risk group, which may facilitate tumor invasion and metastasis.

Predictive drug sensitivity analysis revealed that high-risk patients were more responsive to nine chemotherapeutic agents: Elesclomol, Imatinib, Nilotinib, Roscovitine, Shikonin, Temsirolimus, Nutlin-3a, EHT.1864, and LFM-A13. Consequently, risk scores can guide physicians in offering more personalized treatment plans. Notably, the significant efficacy of Elesclomol, a cuproptosis inducer, in the high-risk group suggests a promising direction for the future development of chemotherapeutic drugs.

Subsequently, we analyzed the expression of the six identified cuproptosis-related hub genes across different cell types, isolating epithelial cells to compare gene expression between normal and tumor epithelial cells. Among these, both OSBPL1A and P4HA1 demonstrated differential expression consistent with quantitative PCR findings, with P4HA1 showing the most pronounced and statistically significant difference. This distinction, coupled with P4HA1’s highest hazard ratio (HR = 2.050, p < 0.001) in multivariate regression analysis, highlighted its potential as a high-risk prognostic marker. Additionally, P4HA1 had a high coefficient in the cuproptosis-related model and a strong correlation with the cuproptosis score (R = 0.41, p < 0.001), suggesting its critical involvement in CRC progression. Furthermore, despite its potential significance, studies on P4HA1 in colorectal cancer remain limited, and its role in cuproptosis has not been explored. These findings collectively justify our decision to prioritize P4HA1 for further in-depth investigation into its role in colorectal cancer and its connection to the cuproptosis pathway. Initially, pseudo-temporal trajectory analysis of epithelial cells revealed an increase in the number of cells with high P4HA1 expression as non-malignant cells transitioned to malignant ones, indicating the critical role of P4HA1 in tumorigenesis. Differences in cellular communication were then investigated between groups with high and low P4HA1 expression, revealing a significant increase in the number and strength of intercellular networks in the high-expression group. Notably, the HGF and NRG signaling pathways, both crucial in tumor progression, were markedly more active in the high-expression group compared to the low-expression group60. KEGG and GO enrichment analyses of differentially expressed genes between the high and low P4HA1 expression groups showed that these genes were predominantly enriched in neurological disease pathways, including Parkinson disease, Huntington disease, Alzheimer disease, and amyotrophic lateral sclerosis, which have strong associations with copper61. Additionally, these genes were enriched in pathways related to oxidative phosphorylation and chemical carcinogenesis involving reactive oxygen species, consistent with the mechanism of cuproptosis affecting mitochondrial respiration and inducing ROS production. GO analysis further showed enrichment in oxidative phosphorylation, aerobic respiration, the electron transport chain, and mitochondrial protein-containing complexes, which aligns with our previous findings.

Therefore, this study aims to further elucidate the role of P4HA1 in colorectal cancer and its potential effect on cuproptosis through experimental investigations. We found that P4HA1 knockdown in colorectal cancer cells led to decreased migration, invasion, and proliferation in vitro, which aligns with our initial hypotheses. Furthermore, gene knockdown significantly enhanced cellular cuproptosis following the exogenous addition of Elesclomol and copper, as evidenced by a marked reduction in cell viability, loss of the iron-sulfur cluster protein FDX1, elevated HSP70 levels, and a substantial increase in intracellular copper content. These changes were accompanied by reduced glutathione levels, increased reactive oxygen species production, and decreased mitochondrial membrane potential. Consequently, it is plausible that P4HA1 knockdown increases the susceptibility of colorectal cancer cells to cuproptosis and associated oxidative stress. This heightened sensitivity to cuproptosis could be therapeutically exploited, offering a novel approach to colorectal cancer treatment.

However, our study has certain limitations. Firstly, we focused on the detailed analysis of only one gene, P4HA1, and did not explore the roles of the other five genes. Secondly, our research was confined to the use of public databases and in vitro experiments, without conducting any in vivo studies. Thirdly, although we identified the potential regulatory role of P4HA1 in cuproptosis within colorectal cancer, the specific mechanisms underlying this process were not fully elucidated, which will be the focus of our future research. Thus, this study explored the possible presence of cuproptosis-related molecules in colorectal cancer and their impact on tumor progression.

Conclusion

In conclusion, we developed a robust prognostic prediction model comprising six genes using public databases, demonstrating high accuracy in predicting clinical outcomes. Among these genes, P4HA1 was selected for further analysis using single-cell RNA sequencing data. In vitro experiments confirmed that P4HA1 knockdown inhibits colorectal cancer cell proliferation and enhances sensitivity to cuproptosis, potentially mediated by oxidative stress induced by this process. The model presented in this study shows promise as a potential biomarker and therapeutic target for colorectal cancer, with significant clinical implications. Additionally, targeting P4HA1 to increase colorectal cancer cell sensitivity to cuproptosis may offer a novel therapeutic strategy for managing colorectal cancer in the future.