Identification of Novel Cuproptosis-Related Genes Mediating the Prognosis and Immune Microenvironment in Cholangiocarcinoma

Abstract

Background

Cuproptosis is a novel type of mediated cell death strongly associated with the progression of several cancers and has been implicated as a potential therapeutic target. However, the role of cuproptosis in cholangiocarcinoma for prognostic prediction, subgroup classification, and therapeutic strategies remains largely unknown.

Methods

A systematic analysis was conducted among 146 cuproptosis-related genes and clinical information based on independent mRNA and protein datasets to elucidate the potential mechanisms and prognostic prediction value of cuproptosis-related genes. A 10-cuproptosis-related gene prediction model was constructed, and its effects on cholangiocarcinoma prognosis were significantly connected to poor patient survival. Additionally, the expression patterns of our model included genes that were validated with several cholangiocarcinoma cancer cell lines and a normal biliary epithelial cell line.

Results

First, a 10-cuproptosis-related gene signature (ADAM9, ADAM17, ALB, AQP1, CDK1, MT2A, PAM, SOD3, STEAP3, and TMPRSS6) displayed excellent predictive performance for the overall survival of cholangiocarcinoma. The low-cuproptosis group had a significantly better prognosis than the high-cuproptosis group with transcriptome and protein cohorts. Second, compared with the high-risk and low-risk groups, the 2 groups displayed distinct tumor microenvironments, reduced proportions of endothelial cells, and increased levels of cancer-associated fibroblasts based on CIBERSORTx and EPIC analyses. Third, patients’ sensitivities to chemotherapeutic drugs and immune checkpoints revealed distinctive differences between the 2 groups. Finally, in replicating the expression patterns of the 10 genes, these results were validated with quantitative real-time polymerase chain reaction results validating the abnormal expression pattern of the target genes in cholangiocarcinoma.

Conclusions

Collectively, we established and verified an effective prognostic model that could separate cholangiocarcinoma patients into 2 heterogeneous cuproptosis subtypes based on the molecular or protein characteristics of 10 cuproptosis-related genes. These findings may provide potential benefits for unveiling molecular characteristics and defining subgroups could improve the early diagnosis and individualized treatment of cholangiocarcinoma patients.

Introduction

Cholangiocarcinoma (CCA) is a rare and lethal disease accounting for nearly 3% of digestive neoplasias and 15% of primary liver cancers.^1–3 However, the global incidence and mortality rates of this fatal disease have increased in recent decades.^4–6 Due to the absence of obvious clinical symptoms or characteristics at the early stages, most CCA patients are diagnosed at advanced stages with limited therapeutic options, resulting in a disappointing prognosis with a median overall survival (OS) duration of less than 18 months.^7–10 Despite advances in curative therapeutics for CCA, surgery is not suitable for patients at advanced stages, and the operative resection rate is no more than 30%, with a high recurrence rate after resection.^3,9 Although individualized medicine and precise therapies have improved for various tumors in recent decades, the characteristics of the biological mechanisms in CCA remain largely unknown.

Previous studies have indicated that different forms of cell death are related to distinctive clinical outcomes of cancer. Necroptosis and pyroptosis have been reported to play important roles in increasing antitumor functions, such as restraining tumorigenesis and metastasis.^11–13 In addition, therapy-resistant cancer cells are vulnerable to detecting ferroptosis signals, especially in the mesenchymal state, and are prone to metastasis. ¹⁴ Badgley et al showed that a genetically modified ferroptosis sensitivity mouse model dramatically decreased the formation and progression of pancreatic cancer. ¹⁵ Recently, cuproptosis was shown to be a novel regulated form of cell death that depends on excess intracellular Cu ions, which induces accumulation of lipoylated dihydrolipoamide S-acetyltransferase and disruption of the mitochondrial tricarboxylic acid (TCA) cycle, resulting in proteotoxic stress. ¹⁶ Interestingly, the process of cuproptosis in various cancer cell lines was not repressed when necroptosis, pyroptosi,s and ferroptosis were inhibited via genetic and pharmacological methods.^16,17 Previous studies have reported that increased Cu ion concentrations were pervasively identified in tumor tissue and serum in numerous kinds of cancers.^18,19 Compared with healthy cell lines, cancer cell lines have a higher requirement for Cu ionophores to maintain propagation and metastasis. ²⁰ Although cuproptosis and cuproptosis-related genes (CRGs) are involved in the activation of cancer cell proliferation- and metastasis-related pathways, few studies have implicated the therapeutic functions and molecular characterization of cuproptosis in the progression of CCA. Hence, illuminating the key mechanisms of cuproptosis in oncology and its potential association with clinical outcomes may facilitate the treatment and prognostic prediction of CCA patients.

In the current study, we aimed to comprehensively investigate the prognostic and therapeutic values of CRGs in CCA. The expression patterns of CRGs in tumor tissues and adjacent normal specimens were examined with The Cancer Genome Atlas (TCGA) and 2 independent datasets from the Gene Expression Omnibus (GEO). Additionally, based upon the expression data of CRGs and the clinical information of CCA patients from the National Omics Data Encyclopedia (NGDC) dataset and 2 other independent datasets, we constructed and validated a 10-CRG signature to stratify CCA patients into 2 heterogeneous cuproptosis subtypes with specific prognostic outcomes. Notably, these 2 groups presented discrepancies in functional and biological processes, the immune microenvironment and drug sensitivity. Furthermore, independent cohort analysis and quantitative real-time polymerase chain reaction (qRT-PCR) were performed to verify the expression patterns of the genes involved in the prognostic model. Our findings presented the prognostic value of a novel 10-CRG prediction model and provided novel insights into the potential mechanism as well as therapeutic opportunities in CCA.

Materials and Methods

Data and Resources

The count-based gene expression profiles and associated clinical information of the TCGA-CCA cohort (n = 45), which included 9 adjacent normal tissue samples and 36 tumor tissue samples, were downloaded from the UCSC Xena browser (https://gdc.xenahubs.net). ²¹ Two high-throughput sequencing datasets, GSE107943 and GSE76297, were collected from GEO as validation datasets (https://www.ncbi.nlm.nih.gov/geo/). GSE107943 contains 30 CCA tissue samples and 27 adjacent normal tissue samples, and GSE76297 includes paired CCA tumor tissues and adjacent controls from 91 CCA patients. However, OS information was not reported in the GSE76297 dataset. The clinical information and high-throughput sequencing datasets of patients included in this study were complete and without missing data. To develop and replicate the accuracy and sensitivity of the novel prognostic model, the transcriptomic and proteomic profiles and associated clinical information of 255 CCA patients were acquired from the NGDC dataset (https://www.biosino.org/node/project/detail/OEP001105). To ameliorate potential false-positive results, lowly expressed genes with average counts less than 5 were excluded from downstream analyses. Additionally, to ensure the accuracy of the results and limited the false positive drive from heterogeneity among samples, an integrative bioinformatics analyses were applied to assess the essential trait-relevant CRGs with discovery cohort and validation cohort in parallel. Therefore, we identified and verified an accurate and sensitive 10-CRG signature to predict the OS of CCA patients.

Identification of Survival-Related Genes and Construction of a Prognostic Signature for CCA

A total of 146 CRGs were selected from previous articles and the Molecular Signature Database (MsigDB) v7.0 dataset (http://www.gsea-msigdb.org/gsea/msigdb/) to assess the expression patterns between tumor and adjacent tissue.^16,18,22,23 We performed differential expression gene (DEG) analyses using the limma package (v. 3.48.3). ²⁴ Genes with P values <.05 and log2-transformed fold changes with concordant direction in the GSE76297, GSE107943, and TCGA-CCA datasets were defined as candidate CRGs to perform subsequent bioinformatics analyses. In addition, expression data from OEP001105 were used to identify the association between the abovementioned CRGs and OS. A univariate Cox proportional hazards regression model was conducted to screen prognosis-related genes through the R package glmnet (v 4.1.2). ²⁵ P values <.05 were considered statistically significant. The intersecting CRGs in both DEG analysis and univariate Cox regression were selected for further analyses. To decipher the potential biological functions of the overlapping genes, protein‒protein interaction network analysis with selected targets was performed using the STRING dataset (v.10.5) with default parameters (https://string-db.org).

A LASSO-penalized Cox regression model with 10-fold cross-validation was performed to develop the model with the best predictive performance. According to the linear combination of the regression coefficients from LASSO, a risk score was calculated for each patient. The patients were divided into high- and low-risk groups based on the median value of all risk scores. Kaplan-Meier analysis (K-M) and time-dependent receiver operating characteristic (ROC) curves were employed to assess the sensitivity and specificity of the constructed model. A P value <.05 in the K-M analysis was considered significant. To test the prognostic risk score and assess its predictive precision and sensitivity in different populations, 30 patients in GSE107943 and 36 patients in the TCGA-CCA dataset were classified into low- and high-risk groups based on the median risk score. The K-M and ROC analyses were conducted by using the R packages timeROC (v 0.4), survival (v 3.2-11), survminer (v 0.4.9), and survivalROC (v 1.0.3) in both the training and validation datasets. Additionally, the relationships between the expression levels of the genes in the prognostic signature, including mRNAs, and the OS rate were analyzed in the datasets. The K-M survival curves were drawn with a modified drawing function “ggsurvplot” in survminer (v 0.4.9).

Pathway Enrichment Analysis

To illuminate the potential pathologies of CCA, the R package ClusterProfiler (v. 3.14.3) was applied to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. ²⁶ Items with false discovery rate (FDR) P values <.05 were considered significantly enriched pathways.

Assessment of Immune Cell Infiltration and Stromal Cells

The CIBERSORTx algorithm was applied to evaluate the relative abundance of 22 types of tumor-infiltrating immune cells based on gene expression in the OEP001105 dataset.^27,28 The categories of immune cell types included naive B cells, memory B cells, plasma cells, resting memory CD4⁺ T cells, activated memory CD4⁺ T cells, naive memory CD4⁺ T cells, CD8⁺ T cells, follicular helper T cells, regulatory T cells, gamma delta T cells, resting natural killer cells, activated NK cells, monocytes, M0-M2 macrophages, resting mast cells, activated MCs, resting dendritic cells, activated dendritic cells, eosinophils, and neutrophils. Because of the limited number of samples in the separate datasets, we combined the 2 datasets to identify the immune-related cell proportions. In addition, the proportions of cancer-associated fibroblasts (CAFs) and ECs were calculated using the online EPIC algorithm (https://gfellerlab.shinyapps.io/EPIC_1-1/). ²⁹

Drug Sensitivity Differences in the Risk Score Classification

Drug sensitivity to chemotherapeutic drugs and novel therapeutic drugs for CCA patients was evaluated according to the Genomics of Drug Sensitivity in Cancer database (https://www.cancerrxgene.org/). The half-maximal inhibitory concentration (IC₅₀) of each drug for CCA was predicted by the R package OncoPredict (v. 0.2) with the RNA-seq expression matrix. ³⁰ The Wilcoxon test was applied to calculate drug sensitivity differences between the low-risk and high-risk groups in the training and validation cohorts.

Independent Prognostic Value Analysis

Univariate and multivariate Cox regression analyses of the risk score of the 10-CRG prediction model and conventional clinical features, such as tumor size diameter, carbohydrate antigen 19-9 (CA19-9), carcinoembryonic antigen (CEA), tumor stage, liver fluke infection, sex, and age, were applied to evaluate independent risk factors for CCA patients. Only P values <.05 were considered statistically significant.

Single-Cell RNA Sequencing Analysis of the CCA Tumor Microenvironment

We performed single-cell analysis with 7 CCA patients extracted from the GSE189903 dataset to evaluate the expression levels of 10 CRG-related markers in T cells, B cells, hepatocytes, cholangiocytes, tumor-associated macrophages (TAMs), CAFs, and tumor-associated endothelial cells (TECs). ³¹ The Seurat package (v. 4.3.0) was applied to remove the genes expressed in at least 3 cells and cells with more than 300 genes.

Cell Culture and Reverse Transcription-Quantitative PCR

Human intrahepatic biliary epithelial cells (BeNa Culture Collection Co. Ltd) and human CCA cell lines (HuCCT1, HCCC9810, and RBE; HaoKE Boitechnology Co. Ltd) were cultured in RPMI-1640 containing 10% fetal bovine serum, 1% penicillin G, and 100 μg/mL streptomycin solution under standard conditions (37 °C in 5% CO₂). These cell lines were negative for mycoplasma contamination.

The aberrant expression levels of the genes included in the prediction model were validated by using qRT-PCR. Total RNA was extracted from each cell line using TRIzol (Thermo Fisher Scientific Inc), and reverse transcription reactions were performed with the PrimeScript^TM II first Strand cDNA Synthesis Kit (TaKaRa) following the manufacturer's protocol. The cDNA quality was assessed with an Agilent 2100. Low-quality sequences were rejected for further analysis. The RT-PCR was performed in triplicate with PowerUp SYBR-Green Master Mix (Thermo Fisher Scientific Inc). The qRT-PCR was performed with UltraSYBR Mixture (CWBIO). The primers were purchased from Shanghai Sangon Biological Engineering Technology, and the sequences are listed in Supplemental Table 1. The transcriptome expression levels were assessed relative to the expression of β-actin by the 2^−ΔΔCt method.

Statistical Analysis

All analyses and creation of graphs were performed in the R package Tableone (v. 0.13.0) or SPSS (v. 22.0). Statistically significant differences were evaluated by Wilcoxon or Student t tests. The log-rank test was carried out to evaluate differences in survival between low-risk and high-risk groups with these independent cohorts. The correlation matrix was constructed with R according to the Pearson correlation coefficient.

Result

Identification of Candidate CRGs

In the present study, the analysis processes were implemented according to the workflow shown in Figure 1. In total, 146 CRGs were included for further analyses. Compared with tumor tissues and adjacent controls, 114, 53, and 98 pyroptosis-related genes were aberrantly expressed in the GSE107943, GSE76297, and TCGA-CCA cohorts, respectively (Supplemental Table 2). Among these targets, 47 significant CRGs overlapped in a consistent direction in all datasets (Figure 2A). As shown in Figure 2B, a univariate Cox regression analysis was conducted using 255 CCA patients from the OEP001105 cohort to inspect prognostic biomarkers. A total of 50 prognostic genes were detected to be significantly associated with OS (Supplemental Table 3). Among the cuproptosis-related DEGs and prognostic markers, 16 overlapping genes, 8 upregulated and 8 downregulated, were annotated as candidate CRGs for CCA (Figure 2B). To gain better insight into the molecular mechanisms involved in CCA, we constructed a protein‒protein network and investigated the orchestration of the candidate CRGs for CCA (Figure 2C). In addition, we evaluated the correlations among the candidate genes’ transcriptome expression profiles using a correlation matrix (Figure 2D). Interestingly, these candidate CRGs exhibited strong correlations with each other, which indicated a close interaction and exerted their effects by producing costimulatory or inhibitory signals in biological processes.

Figure 1.

Workflow presenting the steps of constructing the novel 10-gene prognostic model of CCA in the study. CCA, cholangiocarcinoma.

Figure 2.

Identification of specific DEGs in bile duct carcinoma. (A) Venn diagram analysis of DEGs in CCA with the training dataset and validation dataset. DEGs with P < .05 were considered significant. (B) Forest plot revealing the hazard ratio of candidate CRGs. (C) A protein‒protein interaction network of 17 CGRs in the STRING database. (D) Correlation analysis was implemented to identify the relationship based on transcriptomes. CCA, cholangiocarcinoma; CRG, cuproptosis-related gene; DEG, differential expression gene.

Construction of the Prognostic Prediction Model

To assess the joint effect of the prognostic significance of the 16 candidate CRGs on CCA patient survival, we constructed a multigene prognostic model by performing LASSO penalized regression analysis in the OEP001105 cohort. Consequently, LASSO modeling identified a 10-DEG-based signature with the minimum optional lambda value (Figure 3A and B), which included ADAM metallopeptidase domain 9 (ADAM9), ADAM metallopeptidase domain 17 (ADAM17), albumin (ALB), aquaporin 1 (AQP1), cyclin-dependent kinase 1 (CDK1), metallothionein 2A (MT2A), peptidylglycine alpha-amidating monooxygenase (PAM), superoxide dismutase 3 (SOD3), STEAP3 metalloreductase (STEAP3), and transmembrane serine protease 6 (TMPRSS6). The risk score for CCA was calculated as follows: (6.83 × 10⁶ × expression level of ADAM9) + (1.21 × 10⁵× expression level of ADAM17) + (−4.27 × 10⁷ × expression level of ALB) + (−1.43 × 10⁵ × expression level of AQP1) + (0.0898670 × expression level of CDK1) + (1.76 × 10⁵× expression level of MT2A) + (1.6696174 × expression level of PAM) + (0.7596718 × expression level of SOD3) + (0.1891796 × expression level of STEAP3) + (0.1891796 × expression level of TMPRSS6).

Figure 3.

Effectiveness of the prognostic signature for CCA patients. (A and B) LASSO Cox regression identified 10 CRGs related to OS. (C) A heatmap showing the expression level and distribution between the high-risk and low-risk groups with the OEP001105, GSE107943, and TCGA-CCA datasets. (F-H) OS curves comparing the high-risk group and low-risk group defined by the prognostic model in CCA in the training dataset and validation dataset. (I-K) ROC curve analysis of the 10-CRG prognostic model for predicting OS in the training and validation datasets. CCA, cholangiocarcinoma; CRG, cuproptosis-related gene; OS, overall survival; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic.

All CCA patients were divided into high-risk and low-risk groups for the discovery (high-risk: n = 127, low-risk: n = 128) and validation datasets (GSE107943: high-risk: n = 15, low-risk: n = 15; TCGA-CCA: high-risk: n = 18, low-risk: n = 18) according to corresponding median values of the risk score. The expression patterns of the 10-CRG signature in the high- and low-risk groups of the OEP001105 dataset were visualized (Figure 3C-E). The K-M curves were generated to compare survival conditions in the 2 risk groups. The K-M analysis indicated that the low-risk group had a better OS rate than the high-risk group in the training dataset (Figure 3F, P < .0001). In addition, in line with the K-M analysis in the training dataset, the high-risk group was associated with a worse OS rate in the validation datasets (Figure 3G, GSE107943: P = .0092; Figure 3H, TCAG-CCA: P = .04). Receiver operating characteristic curves were applied to assess the precision of the 10-CRG prognostic model in predicting patient survival. Remarkably, the ROC curve analysis in the discovery cohort showed that the areas under the ROC curve (AUC) values at 1, 2, and 3 years were up to 0.755, 0.76, and 0.743, respectively (Figure 3I). In the validation dataset, the AUC values at 1, 2, and 3 years were 0.912, 0.786, and 0.767 for the GSE107943 dataset and 0.838, 0.845, and 0.794 for the TCGA-CCA dataset, respectively (Figure 3J and K). Therefore, these results indicated that the CRG signature could predict CCA patients with high sensitivity and specificity.

Assessment of Prognostic Factors in CCA

To better understand the relevance and underlying mechanisms of the novel prognostic model in CCA, we applied univariate and multivariate Cox regression analyses to evaluate whether the 10-gene model was an independent prognostic factor of OS for CCA patients. Associations between the clinicopathological features and the risk score of the 10-CRG scoring system were assessed in the OEP001105 cohort. The distribution of the risk score, patient survival status, and survival time of CCA patients is visualized in Figure 4A and B. To analyze the clinical characteristics of CCA, tumor size diameter as well as the expression levels of CA19-9 and CEA were significantly higher in the high-risk group than in the low-risk group (P = .015, P = .009, P = .007, Figure 4C). Furthermore, compared with the early stage and advanced stage of low-risk CCA patients, the risk score was dramatically upregulated in the high-risk group (P < .0001, Figure 4D). Additionally, to test whether the developed 10-CRG signature risk score was independent of other clinicopathological features of CCA, we performed univariate and multivariate Cox regression analyses for risk score, including fluke infection, tumor size diameter, CA19-9, CEA, TNM stage, sex, and age (Figure 4E and F). Both univariate and multivariate Cox regressions revealed that tumor size diameter and risk score were independent prognostic predictors for CCA.

Figure 4.

Establishment of a cuproptosis-related risk score for predicting prognosis. (A and B) Distribution of the cuproptosis-related risk score and survival status in CCA patients. (C) Swarm plots showing the correlation of tumor size, expression level of CA19-9 and expression level of CEA and different groups in CCA. (D) Correlation of tumor stage and risk score in CCA. (E and F) Univariate and multivariate Cox regression analyses of 10 CRGs and clinical features. CCA, cholangiocarcinoma; CRG, cuproptosis-related gene; CA19-9, carbohydrate antigen 19-9; CEA, carcinoembryonic antigen.

The Risk Model Predicted the Infiltration of Immune Cells in CCA

To recognize the indicative role of prognostic signals in the tumor microenvironment (TME), we implemented CIBERSORT to evaluate the infiltration of 22 immune-related cells in CCA in the OEP001105 cohort (Figure 5A). Then, the differential distribution of the detected immune-related cells between the high-risk and low-risk groups was examined by the Wilcoxon rank-sum test. Compared with those in the low-risk group, the proportions of activated memory CD4T cells, follicular helper T cells, monocytes, M0 macrophages, activated mast cells, and neutrophils were significantly increased in the high-risk group (P = 3.48 × 10⁻⁰³, P = 6.80 × 10⁻⁰⁴, P = 4.02 × 10⁻⁰², P = 8.30 × 10⁻⁰⁴, P = 3.80 × 10⁻⁰⁷, P = 1.42 × 10⁻⁰³, respectively). Naive B cells, CD8 T cells, resting memory CD4 T cells, activated NK cells, M1 macrophages, and resting dendritic cells were significantly highly infiltrated in the low-risk group (P = 5.40 × 10⁻⁰⁴, P = 2.41 × 10⁻⁰², P = 1.90 × 10⁻⁰⁶, P = 2.10 × 10⁻⁰⁶, P = 2.96 × 10⁻⁰³, P = 8.27 × 10⁻⁰³, respectively). In addition, as shown in Figure 5B, based on the EPIC analysis, CAF scores were significantly higher and endothelial scores were decreased in the high-risk group compared with the low-risk group (P = .01, P = .003).

Figure 5.

Analysis of immune cell infiltration in CCA. (A) The proportion of 22 immune-related cells in CCA was assessed through the CIBERSORTx algorithm in the CCA cohort. Significant differences between the 2 groups were evaluated using the Wilcoxon test. Asterisks indicate significance, and the Y-axis indicates the percentage of immune cells in CCA. Red triangle: continual increase; blue triangle: continual decrease. (B) Expression of cancer-associated fibroblast score and endothelial cell score in different groups. (C-E) Association of cytokinesis-, immune checkpoint-, and chemotherapy-related genes with the 10-CRG prediction model (*P < .05; **P < .01; ***P < .001). (F) Correlation matrix of cytokinesis-, immune checkpoint-, and chemotherapy-related genes and immune cells in CCA. CCA, cholangiocarcinoma; CRG, cuproptosis-related gene.

Potential Therapeutic Targeting and Tertiary Lymphoid Structure Analyses in CCA

To identify the drug sensitivity differences of chemotherapeutic drugs and novel inhibitors for CCA patients between the high-risk and low-risk groups, we investigated the expression patterns of hub genes for chemokines, immune checkpoints, and chemotherapeutic regimens in the field of CCA treatment. A total of 17, 8, and 7 hub targets for chemokines, immune checkpoints, and chemotherapeutic regimens were acquired from previous studies to evaluate the therapeutic effects between the 2 groups.^32,33 Among the 17 chemokine-specific markers, the expression levels of 5 genes (CCL17, CCL19, CCL21, CCL7, and CCL8) were significantly decreased, and those of 2 markers (CCL13 and CCL24) were significantly increased in the high-risk group compared with the low-risk group (Figure 5C). Moreover, immune checkpoint and chemotherapeutic regimen analyses showed that HAVCR2, IDO1, LAG3, PD-L1, CDK1, CLSPN, PARP1, and RAD51 were significantly higher in the high-risk group than in the low-risk group (Figure 5D and E). Subsequently, to infer the potential functions of the hub genes, we performed correlation analysis between the immune-related infiltrating cells and the hub genes for chemokines, immune checkpoints, and chemotherapeutic regimens (Figure 5F).

In addition, the mRNA expression matrix was used for estimation with OncoPredict (v. 0.2). ³⁰ Bonferroni correction was used to define significant drugs (ie, P = .05/198 = 2.5 × 10⁻³), where 198 is the total number of drugs evaluated for the software (Supplemental Table 4). As a result, AGI-6780, AZD1208, AZD8055, cediranib, ibrutinib, ML323, OSI-027, PAK-5399, PCI-34051, PD325901, PFI3, picolinic-acid, rapamycin, SCH772984, selumetinib, and Wnt-C59 were lower in the low-risk group than in the high-risk group, with P values ranging from 1.80 × 10⁻⁴ to 1.90 × 10⁻⁸, illuminating that the low-risk group was more prone to sensitivity to these drugs than the high-risk group (Figure 6A-P).

Figure 6.

Estimated IC₅₀ value differences of chemotherapeutic drugs and novel inhibitors between the high-risk and low-risk score groups (A-P).

Validation of the Risk Score at the Protein Expression Level

To further replicate the sensitivity and accuracy of our prediction model, we performed a similar analysis with the target proteins. A total of 8320 proteins were detected in the OEP001105 cohort, and 8 of the 10-CRG signatures were found in this dataset (ADAM9, ADAM17, ALB, AQP1, CDK1, MT2A, PAM, SOD3, and STEAP3; Supplemental Table 5). The scoring system of CCA was conducted to assess the prognosis of CCA in this cohort dependent on these 8 gene protein expression levels. Of note, based on the K-M analysis results, the survival probability was significantly lower in the high-risk group than in the low-risk group (Figure 7A). The AUCs at 1, 2, and 3 years of survival were 0.713, 0.698, and 0.723, respectively (Figure 7B). In addition, Pearson correlation between the biomarker mRNA and protein expression levels was performed to test the relationship. Among the 8 genes, the mRNA and protein expression levels were strongly correlated with each other (ADAM9: r = 0.77, P < 2.2 × 10⁻¹⁶; ADAM17: r = 0.32, P = 4.0 × 10⁻⁰⁶; AQP1: r = 0.79, P < 2.2 × 10⁻¹⁶; CDK1: r = 0.80, P < 2.2 × 10⁻¹⁶; MT2A: r = 0.32, P = 4.9 × 10⁻⁰⁶; PAM: r = 0.62, P < 2.2 × 10⁻¹⁶; SOD3: r = 0.53, P = 1.1 × 10⁻¹⁵; STEAP3: r = 0.75, P < 2.2 × 10⁻¹⁶; Figure 7C). Furthermore, we performed differential expression protein analysis between the high-risk and low-risk groups via the limma package. A total of 131 significantly upregulated and 107 downregulated proteins were identified for subsequent analysis (Figure 7D). ClusterProfiler was used to reveal the underlying differences in functions between the low-risk and high-risk groups based on KEGG pathways. Fourteen CCA-related pathways, such as the PI3K-Akt signaling pathway, chemical carcinogenesis DNA adducts and bile secretion, were significantly associated with the risk score after FDR correction (Figure 7E). These results indicated that our 10-CRG indicator playing an effect on CCA development is closely related to the immune response and hepatobiliary functions.

Figure 7.

Validation of the 10-CRG signal based on the protein level in CCA patients. (A) K-M curve analysis comparing survival time between the high-risk group and low-risk group defined by the prognostic model in CCA. (B) ROC curve analysis of the 10-CRG prognostic model for predicting OS. (C) Scatter plot representing the correlation between mRNA and protein expression levels of ADAM9, ADAM17, ALB, AQP1, CDK1, MT2A, PAM, SOD3, and STEAP3. (D) Volcano plot of protein expression in the OEP001105 cohort. Plotted along the x-axis is the mean of the log2-fold-change, and along the y-axis is the negative logarithm of the log2 P values. Light blue indicates the proteins with significant P values (FDR < 0.05), and dark orange denotes these significant proteins. The horizontal line represents the threshold of the significant P value. (E) Pathway enrichment analysis for differentially expressed proteins using KEGG analysis. CCA, cholangiocarcinoma; CRG, cuproptosis-related gene; OS, overall survival; K-M, Kaplan-Meier; ROC, receiver operating characteristic; KEGG, Kyoto Encyclopedia of Genes and Genomes; FDR, false discovery rate.

Expression Levels of CRG-Related Genes in CCA

To explore the expression patterns of CRG-related genes in different types of infiltrating immune cells, a total of 47,261 cells from 7 CCA patients were classified into T cells, B cells, hepatocytes, cholangiocytes, TAMs, CAFs, and TECs according to 22 markers (Figure 8A and B). Our evidence revealed that MT2A, ADAM9, PAM, and ADAM17 pervasively existed in all kinds of cells, and SOD3, AQP1, and ALB were highly expressed in CAFs, TECs, and hepatocytes, respectively (Figure 8C-L).

Figure 8.

Single-cell transcriptomic pattern of CCA. (A) UMAP shows the T cell, B cell, tumor-associated macrophage, cancer-associated fibroblast, tumor-associated endothelial cell hepatocyte, and cholangiocyte subsets in all samples. (B) Dot plot revealing the expression levels of 22 markers in different subgroups of cells. (C-L) UMAP plot revealing ADAM9, ADAM17, ALB, AQP1, CDK1, MT2A, PAM, SOD3, STEAP3, and TMPRSS6 expression profiles. CCA, cholangiocarcinoma.

Validation of the Prediction Model Including Gene Expression Patterns

According to the qRT-PCR results, compared with the normal cell line, the majority of the target genes were aberrantly expressed in at least one CCA cell line, except PAM and STEAP (P < .05; Figure 9A-J). The expression levels of AQP1, MT2A, and SOD3 were significantly upregulated in the 3 CCA cell lines. ADAM9 and ADAM17 were evaluated in RBE and HCCC9810 cell lines, and CDK1 was upregulated in HuCCT1 and HCCC9810 cells. In addition, TMPRSS6 and ALB were significantly decreased in the HuCCT1 and RBE cell lines and the HCCC9810 cell line, respectively. However, inconsistent results for TMPRSS6 and ALB were presented in the remaining cell lines.

Figure 9.

Expression patterns of the genes included in the prognostic model at the mRNA level evaluated via qRT-PCR. Relative expression differences of (A) ADAM9, (B) ADAM17, (C) ALB, (D) AQP1, (E) CDK1, (F) MT2A, (G) PAM, (H) SOD3, (I) STEAP3, and (J) TMPRSS6 between the human normal intrahepatic biliary epithelial cell line (HIBEC) and 3 human cholangiocarcinoma cell lines (HuCCT1, HCCC9810, and RBE) via qRT-PCR (*P < .05; **P < .01; ***P < .001; ***P < .0001). qRT-PCR, quantitative real-time-polymerase chain reaction.

Discussion

Cholangiocarcinoma is a highly fatal digestive disease with poor diagnosis and prognosis. Molecular prognostic markers, as an effective and sensitive method, may be applied as a beneficial supplement to traditional clinicopathological parameters for bile duct cancer patients to increase prognostic prediction, early diagnostic precision, subgroup identification, and personalized treatment. Due to the tumor heterogeneity of CCA, the molecular characterization and physiological characteristics of CCA are still unclear, and the effect of the predictive model is unsatisfactory. Numerous studies have indicated that cell death is one of the most important mechanisms engaged in tumors, and dysregulation of cell death-related genes is directly associated with prognosis in several types of cancers. As a novel type of programmed cell death, cuproptosis-related pathway signatures and genes have gained increasing prominence in several kinds of cancer.^16,17,34 However, no study has investigated the prognostic and therapeutic values of CRGs for CCA patients. In this study, we systemically detected the classification potential, immune cell infiltration, and prognostic values of CRGs for CCA patients with an integrative bioinformatic analysis and validation through several independent cohorts. However, most of the prognostic prediction models did not perform validation for the included targets. Overall, we implemented complicated analyses for identifying cuproptosis-related subgroups and extending the knowledge about the potential mechanisms of CRGs, which suggests that CRGs may play essential roles and serve as novel targets for therapy.

In the current study, 146 CRGs were selected to evaluate the expression profiles between tumor and adjacent tissues. We identified 47 replicated DEGs overlapping in 3 independent cohorts of CCA, and 17 of them were significantly associated with OS. Furthermore, a novel 10-gene-based signature was constructed by LASSO Cox regression to predict the OS of CCA. Because heterogeneity among different patients has been widely reported to influence the biological mechanisms, management, and pathological features of CCA patients, ⁸ constructing a prognostic biomarker with high precision and sensitivity as well as validation in an independent cohort is urgently needed to improve prognosis. ³⁵ In our study, the prognostic signature was established and validated in several CCA transcriptome and proteome sequencing datasets with high sensitivity and accuracy. Of note, limited studies have validated the association between predictive models and OS at the proteome level. Compared with the protein expression profile, the constructed prognostic prediction model via the transcriptome profile as a dynamic tool fluctuating with tumor progression to reflect the prognosis is much easier to affect by other factors. Strikingly, consistent with the analysis results at the mRNA level, the low-cuproptosis group separated by the novel prognostic model, which was constructed based on the protein level, had a better prognosis than the high-cuproptosis group. Furthermore, the expression levels of 8 out of 10 targets (ADAM9, ADAM17, ALB, AQP1, CDK1, MT2A, SOD3, and TMPRSS6) were verified in at least one CCA cell line compared with normal bile duct cells. However, ALB and TMPRSS6 present conflicting results, possibly triggered by heterogeneity among different CCA cell lines.

Four of these 10 prognostic genes (ADAM17, ALB, AQP1, and CDK1) were previously reported to be associated with CCA prognosis. ALB, as one of the markers of hepatocytes and cholangiocytes, plays oncogenic roles in hepatocyte-specific processes, such as complement activation, detoxification, fatty acid catabolic process, and bile acid secretion. ³⁶ CDK1 regulates the hub genes of cytokinesis-related genes, which cooperate in the process of cellular proliferation, cytoskeleton organization, and cytokinesis.^37,38 AQP1 is a small hydrophobic integral transmembrane protein that is ubiquitously expressed in multiple human cancers. Increasing evidence indicates that AQP1 overexpression accelerates migration, extravasation, tumor angiogenesis, and metastasis in vitro and in vivo.^39,40 Julia and colleagues suggested that TNFR1-dependent tumor cell-induced endothelial cell death, tumor cell extravasation, TNF-induced necroptosis, and subsequent metastatic seeding rely on the activation of the membrane-bound metalloprotease ADAM17. ⁴¹ In addition, all of these targets serve as prognostic biomarkers for many other types of cancers, such as ovarian cancer, colorectal cancer, hepatocellular carcinoma, and kidney cancer.^42–45 However, the biological function of these genes in CCA is still unclear, and further investigation is needed to explore the potential mechanisms of these target genes in the future.

The tumor immune microenvironment is considered a key regulator of carcinogenesis, angiogenesis, and tumor growth and affects the efficacy of radiotherapy, chemotherapy, and immune checkpoint therapy.^46,47 Bile duct carcinoma is characterized by an intricate microenvironment in which the stroma is composed of TECs, fibroblasts, macrophages, neutrophils, NK cells, and T cells. ⁴⁸ According to the infiltration of immune cells in CCA, we detected more naive B cells, CD8 T cells, resting memory CD4 T cells, activated NK cells, M1 macrophages, and resting dendritic cells in the low-risk group. In addition, the risk score for CCA was positively correlated with activated memory CD4 T cells, follicular helper T cells, monocytes, M0 macrophages, activated mast cells, and neutrophils in CCA patients. M0 macrophages are nonactivated macrophages that differentiate into M1/2 macrophages in the presence of specific conditions. M1 macrophage infiltration promotes an inflammatory microenvironment by increasing the expression of IL-1, IL-6 and IL-12, and M2 macrophage infiltration promotes inflammatory escape by increasing the expression of PD-L1 and IL-10.^49,50 A previous study suggested that patients with a higher infiltration of CD8+ T cells and M1 macrophages in the TME have a better prognosis and high immune checkpoint gene expression than those infiltrated by M0/2 macrophages.^51,52 In our results, the microenvironment of the CCA low-risk group had significantly enriched M1 macrophage and CD8+ T-cell infiltration. Previous studies reported that the abundance of CAFs was a major TME with tumorigenic properties, which were positively correlated with tumor growth and poor prognosis.^53,54 Cancer-associated fibroblasts regulate innate immunity by suppressing NK cell activation to promote immunosuppression. ⁵⁵ Moreover, NK cell-induced antibody-dependent cellular cytotoxicity increased cancer cell death, and infusion of NK cells in CCA mouse xenograft models resulted in tumor regression.^56,57 Immune checkpoint blockade with monoclonal antibody therapies has been implemented as a novel plan for treating numerous malignancies with remarkable and durable response rates, ⁵⁸ while clinical trials of immunotherapies in CCA have displayed less notable success. ⁸ NK cells are currently widely targeted in immunotherapy for cancer patients without major histocompatibility complex (MHC) class I. ⁵⁹ When the MHC I expression level is downregulated or silenced in cancer cells, NK cells are activated to remove tumor cells. Laura et al demonstrated that activated mast cells could induce biliary hyperplasia, hepatic fibrosis, and vascular bed dysfunction, and knockout or inhibition of mast cells was likely to be a better therapy for patients with cholangiopathies. ⁶⁰ Interestingly, based on EPIC analysis, this study detected that CAF scores were significantly increased, and endothelial scores were reduced in the high-risk group. In summary, these results strongly highlight that it is necessary to explore the therapeutic value of the TME of M1 macrophages and NK cells in inhibiting CAF therapy for CCA.

For the drug sensitivity analyses, 16 kinds of therapeutic drugs, AGI-6780, AZD1208, AZD8055, cediranib, ibrutinib, ML323, OSI-027, PAK-5399, PCI-34051, PD325901, PFI3, picolinic-acid, rapamycin, SCH772984, selumetinib, and Wnt-C59, revealed more sensitivity in the low-risk group than in the high-risk group. Selumetinib, an inhibitor of mitogen-activated protein and extracellular signal-regulated kinase, can be targeted to suppress tumor growth and reduce IL-6 expression in several different types of cancer.^61,62 Moreover, selumetinib has highly effective effects on inhibiting tumor growth by reducing tumor cell proliferation and accumulating apoptosis in a K-Ras-driven CCA mouse model. ⁶³ In our results, the IC₅₀ value of selumetinib in the low-cuproptosis group was significantly lower than that in the high-cuproptosis group, suggesting that selumetinib may be more effective in the low-risk group. Thus, the novel prognostic model can serve as a personalized treatment indicator for CCA patients.

Although comprehensive analyses were performed to investigate the potential mechanisms of the prognostic biomarkers for CCA, this study has certain limitations. First, immune cell infiltration in bile duct cancer was assessed on the basis of publicly available TCGA, GSE107943, and NGDC RNA-seq data. Further biological experiments need to be implemented in vivo and in vitro to validate the immunotherapeutic value and interaction among immune cells and CRGs in CCA. Second, even if half of the genes included in the prognostic model have been reported to influence different cancers, biological experiments should be performed to elucidate the molecular mechanisms of CRGs in the pathogenesis and prognosis of CCA.

Conclusion

In summary, we identified and verified an accurate and sensitive 10-CRG signature to predict the OS of CCA patients. Mechanistically, cuproptosis plays a key role in CCA development and is closely related to immune cell infiltration, CAFs, drug sensitivity, and checkpoint and chemotherapy marker expression. These findings provide a novel prognostic signature to account for the OS of CCA patients and may benefit the timely diagnosis and individualized treatment of CCA.

Supplemental Material

sj-xlsx-1-tct-10.1177_15330338241239139 - Supplemental material for Identification of Novel Cuproptosis-Related Genes Mediating the Prognosis and Immune Microenvironment in Cholangiocarcinoma

Supplemental material, sj-xlsx-1-tct-10.1177_15330338241239139 for Identification of Novel Cuproptosis-Related Genes Mediating the Prognosis and Immune Microenvironment in Cholangiocarcinoma by Qiang Liu, Jianpeng Zhu, Zhicheng Huang, Xiaofeng Zhang and Jianfeng Yang in Technology in Cancer Research & Treatment