Comprehensive analysis of the relationship between cuproptosis-related genes and esophageal cancer prognosis

Abstract

BACKGROUND

Esophageal cancer is one of the most common malignant tumors of the digestive system, with a 5-year survival rate of 15% to 50%. Cuproptosis, a unique kind of cell death driven by protein lipoylation, is strongly connected to mitochondrial metabolism. The clinical implications of cuproptosis-related genes in esophageal cancer, however, are mainly unknown.

AIM

To identify cuprotosis-related genes that are differentially expressed in esophageal cancer and investigate their prognostic significance.

METHODS

With |log fold change| > 1 and false discovery rate < 0.05 as criteria, the Wilcox test was used to evaluate the differentially expressed genes between 151 tumor tissues and 151 normal esophageal tissues. Cuproptosis-related genes were selected to be linked with prognosis using univariate Cox regression analysis. Genes were separated into high- and low- expression groups based on their cutoff value of gene expression, and the correlation between the two groups and overall survival or progression-free survival was investigated using the log-rank test. The C-index, calibration curve, and receiver operator characteristic (ROC) curve were used to assess a nomogram containing clinicopathological characteristics and cuproptosis-related genes.

RESULTS

Pyruvate dehydrogenase A1 (PDHA1) was found to be highly correlated with prognosis in univariate Cox regression analysis (hazard ratio = 22.96, 95% confidence interval = 3.09-170.73; P = 0.002). According to Kaplan-Meier survival curves, low expression of PDHA1 was associated with a better prognosis (log-rank P = 0.0007). There was no significant correlation between PDHA1 expression and 22 different types of immune cells. Tumor necrosis factor superfamily member 15 (TNFSF15) (P = 3.2 × 10^-6; r = 0.37), TNFRSF14 (P = 8.1 × 10^-8; r = 0.42), H long terminal repeat-associating 2 (P = 6.0 × 10^-8; r = 0.42) and galectin 9 (P = 3.1 × 10^-6; r = 0.37) were all found to be considerably greater in the high PDHA1 expression group, according to an analysis of genes related to 47 immunological checkpoints. The low PDHA1 expression group had significantly lower levels of cluster of differentiation 44 (CD44) (P = 0.00028; R = -0.29), TNFRSF18 (P = 1.2 × 10^-5; R = -0.35), programmed cell death 1 ligand 2 (P = 0.0032; R = -0.24), CD86 (P = 0.018; R = -0.19), and CD40 (P = 0.0047; R = -0.23), and the differences were statistically significant. We constructed a prognostic nomogram incorporating pathological type, tumor-node-metastasis stage, and PDHA1 expression, and the C-index, calibration curve, and ROC curve revealed that the nomogram’s predictive performance was good.

CONCLUSION

Cuproptosis-related genes can be used as a prognostic predictor for esophageal cancer patients, providing novel insights into cancer treatment.

Key Words: Esophageal cancer; Cuproptosis; Pyruvate dehydrogenase A1; Overall survival; Nomogram

Core Tip: Esophageal carcinoma has a poor prognosis and is one of the major causes of cancer-related deaths worldwide. Despite recent advancements in the surgical and pharmacological treatment of esophageal cancer, the prognosis remains poor. Copper toxicity has been linked to the incidence and progression of esophageal cancer in numerous studies. At the gene level, however, the probable biochemical mechanism is unknown. We included 19 cuproptosis-related genes and screened a gene that could successfully predict the prognosis of esophageal cancer by statistical analysis to further elucidate the role of cuproptosis-related genes in impacting the prognosis of esophageal cancer.

INTRODUCTION

Esophageal cancer is a major global health issue, and its incidence is rapidly increasing[1]. Esophageal cancer is classified into two types: esophageal squamous cell carcinoma (ESCC), which accounts for 90% of all occurrences, and esophageal adenocarcinoma (EAC). In recent years, epidemiological research has revealed that the incidence of EAC has grown 3-4 times, with the proportion increasing[2,3]. Despite significant advances in the diagnosis and treatment of esophageal cancer, the mortality rate varies from 15% to 20%, placing it fourth among all cancer-related deaths[1].

Esophageal cancer is the result of a complex process involving various causes and polygene alterations. Using high-throughput sequencing technologies, a comprehensive mutation catalog was evaluated, and substantial genetic alterations have been discovered in the malignancies. Gene alterations are often linked to aberrant expression and are becoming more essential for the early diagnosis and prognosis of esophageal cancer[4]. Currently, several gene expression products are employed as indicators for esophageal cancer diagnosis and prognosis[5,6]. Somatic mutations in tumor protein p53 (TP53) have been found in more than 83% of ESCCs. Adenocarcinoma and squamous cell carcinoma both have TP53 point mutations[7,8]. In addition, numerous cell cycle-controlling genes are overexpressed in ESCC. For example, cyclin-dependent kinase 4/cyclin-dependent kinase 6 accounted for 23.6%, murine double minute 25.7%, and cyclin D1 46.4%, showing that the above components are implicated in the incidence and development of ESCC[9]. As a result, there is a pressing need to uncover genetic anomalies in esophageal cancer and understand their molecular basis to enhance early diagnosis and minimize esophageal cancer mortality.

Previous studies have reported several types of precisely programmed cell death, including apoptosis, pyroptosis, necrosis, and ferroptosis[10]. Similar to iron, copper is a trace metal in cells that plays an integral role in maintaining protein functions. Excessive copper can cause cytotoxicity, but the exact mechanism is unclear[11]. Tsvetkov et al[12] discovered that the copper carrier elesclomol, which was originally used to treat cancer, kills cells in excess of copper. Elesclomol does not trigger cell death on its own, suggesting that copper toxicity is to blame. Dar et al[13] found that patients with esophageal cancer had considerably higher mean blood copper levels than controls, with a mean copper concentration of 169 g/dL in the cancer group and 149 g/dL in the control group. Therefore, we hypothesized that copper shortage or excess is linked to the occurrence and progression of esophageal cancer. Copper, however, has been linked to the development and progression of esophageal cancer in few studies.

There has not been any research on the role of cuproptosis-related genes in the genesis, progression, and metastasis of esophageal cancer to date. Cuproptosis-related genes and their processes need to be better understood in order to improve the prognosis of malignant tumors and uncover novel treatment targets. Bioinformatics analysis was used to evaluate the expression profile of cuproptosis-related genes and its predictive significance in esophageal cancer in this study.

MATERIALS AND METHODS

The study did not include any human participants, data, or tissue, nor did it include any animals. All of the information was gathered from a public database.

Data collection

The Cancer Genome Atlas (https://portal.gdc.cancer.gov/) provided gene expression data and clinical information for 151 esophageal cancer samples, whereas the Genotype-Tissue Expression database (https://xenabrowser.net/) provided gene expression levels for 151 healthy tissue samples. The "limma" package in R was used to conduct matrix normalization. We selected 19 genes from the scientific literature that have been linked to cuproptosis in prior studies[12].

Screening of cuproptosis-related prognosis genes

In 151 patients with esophageal cancer and 151 healthy controls, a total of 55185 genes were acquired. The log2 (x+1) scale was used to standardize all of the expression data. Differentially expressed genes (DEGs) were screened between tumor tissue and normal esophageal tissue using the "limma" package and Wilcox-test, with the conditions of |log fold change [FC]| > 1 and false discovery rate (FDR) < 0.05. Using the "limma" package, we evaluated those genes that were differently expressed between esophageal cancer and healthy control tissues based on 19 cuproptosis-related genes. The FDR less than 0.05 was used as a criterion for further investigation. The researcher employed a univariate Cox regression analysis to find predictive markers linked to overall survival (OS) and progression-free survival (PFS). The optimum cutoff value for a prognostic gene was identified using a receiver operating characteristic (ROC) curve. Patients were classified into high-risk (less than the cutoff) and low-risk (greater than the cutoff) esophageal cancer groups based on gene expression cutoff values.

Function enrichment analysis

The R package "ClusterProfiler" was used to analyze gene ontology (GO) and pathway enrichment [Kyoto Encyclopedia of Genes and Genomes (KEGG)]. It was deemed substantially enriched when the P value and adjusted P value were both less than 0.05. To identify the functional role of the genes, GO analysis was done on the significantly expressed genes, and the expression levels of the genes were displayed in GO circle plots using the R package "GOplot". To determine tumor-infiltrating immune-cell fractions in esophageal cancer patients, we employed the CIBERSORT algorithm[14]. On samples with a CIBERSORT result of P value less than 0.05, further analysis was performed. Then for reference, we collected 547 gene expression profiles from the CIBERSORT website (http://cibersort.stanford.edu/). The Pearson’s test and the "corrplot" program were used to correlate infiltrated immune cells.

Development of a nomogram for predicting OS

A cuproptosis-related gene and numerous clinicopathological variables were also subjected to univariate and multivariate Cox regression analyses. For predicting the OS, we utilized the R "rms" package to generate a nomogram including pathological categories, American Joint Committee on Cancer (AJCC)-tumor-node-metastasis (TNM) stages, and differential cuproptosis-related genes. The uniformity of the nomogram was assessed using calibration curves to anticipate various OS results. Harrell’s concordance index was used to generate the C-index. The area under the receiver operating characteristic (AUC) curve and ROC were both utilized to evaluate how predictive our nomogram was.

Statistical analyses

All statistical analyses were performed using R software, version 4.1.0 (http://www.rproject.org/). The Students t-test was used to compare variables such as age at diagnosis, sex, AJCC-TNM stage, mutational status of v-Ki-ras2 Kirsten rat sarcoma viral oncogene homolog (KRAS), v-Raf murine sarcoma viral oncogene homolog B1 genes (BRAF), epidermal growth factor receptor (EGFR), phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha (PIK3CA, unpaired t-test). The Cox proportional hazards model was used to find genes and clinicopathological variables linked to survival. Based on the optimal cutoff of differential cuproptosis-related genes, the Youden index method determined that esophageal cancer patients should be divided into high- and low-risk groups. For assessing survival in the high- and low-risk groups, Kaplan-Meier survival curves were employed, and log-rank tests were performed to compare survival rates. For two-sided P values, statistical significance was considered as a value less than 0.05.

RESULTS

Demographics of the 151 esophagus cancer patients

Our study comprised 151 esophageal cancer patients and 151 healthy controls. Table 1 shows the characteristics of patients with EAC (n = 74) and ESCC (n = 77). All esophageal cancer patients had a median survival duration of 13.4 mo. EAC had a median survival time of 14.8 mo, whereas ESCC had a median survival time of 13 mo.

Table 1 Demographic characteristics and clinicopathological characteristics.

Variables	Total, n = 151	Esophagus cancer, n = 151
		Adenocarcinoma, n = 74	Squamous cell carcinoma, n = 77
Age, yr
Median	60	68.5	57
Interquartile range	(53-72)	(57-77)	(51-63.5)
Sex
Female	22 (14.6%)	11 (14.9%)	11 (14.3%)
Male	129 (85.4%)	63 (85.1%)	66 (85.7%)
AJCC-TNM Stage
Ⅰ	18 (11.9%)	11 (14.9%)	7 (9.1%)
Ⅱ	70 (46.4%)	24 (32.4%)	46 (59.7%)
Ⅲ	51 (33.8%)	31 (41.9%)	20 (26.0%)
Ⅳ	12 (7.9%)	8 (10.8%)	4 (5.2%)
EGFR status
Mutant	5 (3.3%)	3 (4.1%)	2 (2.6%)
Wild-type	146 (96.7%)	71 (95.9%)	75 (97.4%)
KRAS status
Mutant	2 (1.3%)	2 (2.7%)	0 (0)
Wild-type	149 (98.7%)	72 (97.3%)	77 (100%)
BRAF status
Mutant	1 (0.7%)	0 (0)	1 (1.3%)
Wild-type	150 (99.3%)	74 (100%)	76 (98.7)
PIK3CA status
Mutant	14 (9.3%)	4 (5.4%)	10 (13.0%)
Wild-type	137 (90.7%)	70 (94.6%)	67 (87.0%)
OS event
Event	58 (38.4%)	36 (48.6%)	22 (28.6%)
Non-event	93 (61.6%)	38 (51.4%)	55 (71.4%)
OS months
Median	13.4	14.8	13.0
Range	(7.8-22.9)	(7.3-27.5)	(11.0-18.6)
PFS event
Event	73 (48.3%)	38 (51.4%)	35 (45.5%)
Non-event	78 (51.7%)	36 (48.6%)	42 (54.5%)
PFS months
Median	10.7	10.3	10.7
Range	(5.1-18.9)	(5.4-24.5)	(3.7-15.8)

AJCC: American Joint Committee on Cancer; BRAF: v-Raf murine sarcoma viral oncogene homolog B1 genes; EGFR: Epidermal growth factor receptor; KRAS: v-Ki-ras2 Kirsten rat sarcoma viral oncogene homolog; OS: Overall survival; PFS: Progression-free survival; PIK3CA: Phosphatidylino-sitol-4,5-bisphosphate 3-kinase, catalytic subunit alpha; TNM: Tumor-node-metastasis.

Relationship between clinical parameters and PFS/OS

Table 2 summarizes the relationships between clinical pathological features and these individuals’ OS or PFS. When compared to ESCC, EAC was related to a lower OS (P = 0.011). AJCC-TNM stages III and IV were linked to OS (P = 0.021) and PFS (P = 0.013). When compared to females, males had a lower PFS (P = 0.009). Age, EGFR, BRAF, KRAS, or PIK3CA status, among other clinicopathological variables, were not linked to OS or PFS.

Table 2 Relationship between the clinicopathological characteristics and survival outcome.

Parameters	Overall survival			Progression-free survival
	Non-event	Event	P value	Non-event	Event	P value
Age, yr			0.758			0.628
≤ 65	57	37		50	44
> 65	36	21		28	29
Sex			0.102			0.009
Female	17	5		17	5
Male	76	53		61	68
AJCC-TNM stage			0.021			0.013
I and II	61	27		53	35
III and IV	32	31		25	38
Pathological type			0.011			0.469
Adenocarcinoma	38	36		36	38
SCC	55	22		42	35
EGFR status						0.596
Wild-type	3	2	0.941	2	3
Mutant	90	56		76	70
BRAF status			0.428			0.332
Wild-type	1	0		1	0
Mutant	92	58		77	73
KRAS status			0.734			0.962
Wild-type	1	1		1	1
Mutant	92	57		77	72
PIK3CA status			0.828			0.896
Wild-type	9	5		7	7
Mutant	84	53		71	66

AJCC: American Joint Committee on Cancer; BRAF: v-Raf murine sarcoma viral oncogene homolog B1 genes; EGFR: Epidermal growth factor receptor; KRAS: v-Ki-ras2 Kirsten rat sarcoma viral oncogene homolog; PIK3CA: Phosphatidylino-sitol-4,5-bisphosphate 3-kinase, catalytic subunit alpha; SCC: Squamous cell carcinoma; TNM: Tumor-node-metastasis.

Identification of prognostic cuproptosis-related genes

With |log FC| > 1 and FDR < 0.05, a total of 7055 DEGs were identified in the healthy control group, comprising 3494 upregulated genes and 3561 downregulated genes. Then, using the FDR < 0.05 threshold, gene expression analysis was done to identify cuproptosis-related genes and 18 genes fulfilled our requirements (Supplementary Table 1). In the intersection of two gene sets, there were five genes: high-affinity copper uptake protein 1 (SLC31A1), ferredoxin, lipoyl(octanoyl) transferase 2, pyruvate dehydrogenase A1 (PDHA1), and programmed cell death 6 interacting protein (DEGs and cuproptosis-related genes). Furthermore, using univariate Cox regression analysis, we discovered that one gene (PDHA1) was strongly linked to prognosis [hazard ratio (HR) = 22.96, 95% confidence interval (CI) = 3.09-170.73, P = 0.002] (Table 3).

Table 3 The results of the univariate Cox regression analysis.

Gene name	HR	HR.95L	HR.95H	P value
PDHA1	22.96	3.09	170.73	0.002
ATP7A	3.83	0.96	15.21	0.057
CDKN2A	0.84	0.64	1.10	0.204
PDHB	1.76	0.37	8.50	0.479
GLS	1.50	0.48	4.67	0.484
ATP7B	1.15	0.77	1.71	0.496
FDX1	1.52	0.41	5.59	0.527
DLD	1.82	0.27	12.11	0.534
LIPT2	0.87	0.46	1.64	0.662
DLST	0.58	0.05	6.70	0.662
NFE2L2	0.73	0.18	3.04	0.666
DLAT	1.43	0.23	8.89	0.702
DBT	0.83	0.27	2.54	0.744
GCSH	0.87	0.32	2.36	0.777
MTF1	0.85	0.26	2.80	0.785
LIAS	0.86	0.25	2.93	0.806
NLRP3	0.91	0.44	1.90	0.809
LIPT1	1.14	0.29	4.41	0.853
SLC31A1	1.08	0.19	6.15	0.928

HR: Hazard ratio; HR.95L: Low 95% confidence interval of HR; HR.95H: High 95% confidence interval of HR.

Functional enrichment analysis

Using GO enrichment and KEGG pathway analyses, we explored the biological processes, cellular components, and molecular activities of DEGs. Energy metabolism and glycolysis/gluconeogenesis signaling pathways were highly enriched among DEGs with cuproptosis relevance. The phrases "tricarboxylic acid cycle," "acetylCoA metabolic process," and "acetylCoA biosynthetic process from pyruvate" are considerably enriched in Figure 1A. "Carbon metabolism," "tricarboxylic acid cycle," "pyruvate metabolism," "glycolysis/gluconeogenesis," "platinum drug resistance," "biosynthesis of cofactors," "mineral absorption," and "central carbon metabolism in cancer" were the most substantially enriched pathways (Figure 1B).

Figure 1 Esophageal cancer functional enrichment analysis of cuproptosis-related genes. A: The expression of cuproptosis-related genes is represented by the outer circle in each enriched Gene Ontology (GO) terms; red dots on each GO term shows upregulated cuproptosis-related genes, and the inner circle indicates the importance of GO terms (log10 adjusted P value). The downregulated cuproptosis-related genes are shown by blue dots, and the chart on the right shows the distribution of cuproptosis-related genes in important GO terms; B: The top 16 Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis findings are shown. The GeneRatio denotes the number of differentially expressed genes (DEGs) detected in one GO pathway in proportion to the total number of DEGs.

Prognostic value of PDHA1

The 151 esophageal cancer patients were divided into two groups based on the cutoff value: Those with high PDHA1 expression and those with low PDHA1 expression. In the analysis of OS, high PDHA1 expression was linked to considerably lower OS rates than low expression (P = 0.007; Figure 2A). However, there was no link between PDHA1 expression and PFS (Figure 2B). Within a multivariate context, we investigated the link between PDHA1 expression and survival result (Table 4). When we controlled for clinical prognostic indicators that were significant (P < 0.05) in univariate Cox regression models, we discovered that overexpression of the PDHA1 gene might predict poor clinical outcomes. The expression of PDHA1 (HR: 1.67, 95 CI: 1.03-2.73; P = 0.0386) and AJCC-TNM stage (HR: 2.30, 95%CI: 1.58-3.35; P < 0.001) were independent risk factors for OS, but not for PFS, according to multivariate Cox regression analysis.

Figure 2 The Kaplan-Meier survival of patients with esophageal cancer who were classified according to pyruvate dehydrogenase A1 expression. A and B: In esophageal cancer, patients with low pyruvate dehydrogenase A1 (PDHA1) expression levels had better overall survival (A) and progression-free survival (B) than patients with high PDHA1 expression levels. The significance of the difference between high and low expression was calculated using the log-rank test.

Table 4 Univariate and multivariate Cox regression analyses of pyruvate dehydrogenase A1 expression and the clinicopathological factors in the esophagus cancer patients.

Variables	Overall survival		Progression-free survival
	HR (95%CI of HR)	P value	HR (95%CI of HR)	P value
PDHA1 low	1.00 (Reference)		1.00 (Reference)
UVA PDHA1 high	2.34 (1.42-3.87)	0.0009	1.51 (0.95-2.38)	0.08
MVA PDHA1 high	1.67 (1.03-2.73)	0.0386
Age ≤ 65	1.00 (Reference)		1.00 (Reference)
Age > 65	0.81 (0.47-1.41)	0.46	0.97 (0.60-1.55)	0.89
Sex, Female	1.00 (Reference)		1.00 (Reference)
UVA Male	2.12 (0.84-5.34)	0.11	2.97 (1.19-7.39)	0.02
MVA Male			2.27 (0.90- 5.73)	0.08
Pathological types, adenocarcinoma	1.00 (Reference)		1.00 (Reference)
SCC	0.81 (0.47-1.39)	0.45	1.11 (0.70-1.77)	0.66
AJCC stage I, II	1.00 (Reference)		1.00 (Reference)
UVA AJCC stage III, IV	2.53 (1.76-3.64)	5.86 × 10-7	1.98 (1.46-2.67)	9.48 × 10-6
MVA AJCC stage III, IV	2.30 (1.58-3.35)	1.46 × 10-5	1.86 (1.37-2.52)	6.98 × 10-5

CI: Confidence interval; HR: Hazard ratio; PDHA1: Pyruvate Dehydrogenase A1; UVA: Univariate; MVA: Multivariate; SCC: Squamous cell carcinoma.

Different infiltration levels of immune cells between low and high PDHA1 groups

Our study used the CIBERSORT algorithm to examine the ratio of tumor-infiltrating immune cells (TICs) in esophageal cancer to further verify the correlation between PDHA1 and immune cells (Figure 3A). P < 0.05 was used to classify samples as statistically different. The ratio of TICs in the low-risk group was represented by the first 108 of 151 esophageal cancer patients, whereas the ratio of TICs in the high-risk group was represented by the last 43 samples. The correlation between 22 immune cells was represented using a heatmap (Figure 3B). Neutrophils and activated mast cells, macrophages M1 and activated natural killer cells, plasma cells and naive B cells, activated mast cells and activated dendritic cells, and macrophages M1 and naive CD4⁺ T cells were the top five results with a positive correlation. CD8⁺ T cells and resting memory CD4⁺ T cells, on the other hand, were the immune cells that were most negatively associated. Furthermore, the bar graph revealed that the immune infiltration of 22 immune cell types in high and low PDHA1 esophageal cancer patients did not statistically significantly differ (Figure 3C).

Figure 3 Tumor-infiltrating immune cell profile. A: A barplot shows the ratio of 22 different tumor-infiltrating immune cells in esophageal cancer patients. Column names (X-axis): sample ID; B: A heatmap shows the association between 22 different types of tumor-infiltrating immune cells. Each dot represents the P value of the correlation between two different cell types, and Pearson’s correlation coefficient was used to determine significance; C: A bar graph shows the difference between esophageal cancer patients with high or low expression of pyruvate dehydrogenase A1 in terms of 22 different types of tumor-infiltrating immune cells. The Wilcoxon rank-sum test was used for the significance test.

PDHA1 expression and immune checkpoint correlation analysis

We also investigated the relationships between PDHA1 and 47 genes associated with immunological checkpoints, which have been reported in the literature[15]. The results showed that tumor necrosis factor superfamily member 15 (TNFSF15) (P = 3.2 × 10^-6; r = 0.37), TNFRSF14 (P = 8.1 × 10^-8; r = 0.42), H long terminal repeat-associating 2 (HHLA2) (P = 6.0 × 10^-8; r = 0.42) and galectin 9 (LGALS9) (P = 3.1 × 10^-6; r = 0.37) were significantly higher than that of the low expression group of PDHA1. Cluster of differentiation 44 (CD44) (P = 0.00028; R = -0.29), TNFRSF18 (P = 1.2 × 10^-5; R = -0.35), programmed cell death 1 ligand 2 (PDCD1LG2) (P = 0.0032; R = -0.24), CD86 (P = 0.018; R = -0.19), and CD40 (P = 0.0047; R = -0.23) were significantly lower than that of the low expression group of PDHA1 (Figure 4A). The expression profiles of 47 immune checkpoint genes involved in cuproptosis were investigated further. In addition, we categorized the above genes into high and low expression groups based on their cutoff values. Kaplan-Meier analysis and log-rank testing were used to assess OS. The findings revealed that low HHLA2, TNFRSF18, and CD44 overexpression was linked to a significantly shorter OS and a worse prognosis (Figure 4B-D).

Figure 4 Correlation between pyruvate dehydrogenase A1 expression and immune checkpoint-related genes. A: Barplot reveals expression of 47 immune checkpoint-related genes in high- and low-pyruvate dehydrogenase A1 (PDHA1) esophageal cancer patients; B-D: H long terminal repeat-associating 2 (B), tumor necrosis factor superfamily member 18 (C), and cluster of differentiation 44 (D) Kaplan-Meier survival curves of high and low expression of immune checkpoint-related genes. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

Constructing and evaluating a predictive nomogram

The 1-, 2-, and 3-year OS probability was calculated using a nomogram that included the pathological type, AJCC-TNM stage, and PDHA1 expression (Figure 5A). The actual observed vs anticipated rates of the 1-, 2-, and 3-year OS demonstrated near to the ideal 45° oblique line on the calibration curve (Figure 5B). Furthermore, the AUC of ROC curves for 1-, 2-, and 3-year survival were 0.725, 0.776, 0.619, and 0.810, respectively (Figure 5C). With a C-index of 0.703 for OS, the nomogram showed promising discrimination.

Figure 5 Developing a predictive nomogram. A: A nomogram that predicts patients with esophageal cancer’s 1-, 2-, and 3-year survival; B: A calibration curve of 1-, 2-, and 3-year survival in the nomogram and ideal model; C: Receiver operating characteristic curves and area under the curve for the nomogram’s 1, 2, and 3-year survival. AUC: Area under the curve; PDHA1: Pyruvate dehydrogenase A1.

DISCUSSION

In the present study, we investigated the expression signature of 19 cuproptosis-related genes in esophageal cancer tissues and explored their relationships with OS and PFS. A prognostic nomogram involving gene expression and clinicopathological parameters was constructed for the first time. Functional analysis showed that DEGs were enriched in energy metabolism, especially in pathways related to the tricarboxylic acid cycle. Cuproptosis-related genes were also confirmed to be associated with immune checkpoint genes.

In our study, 151 patients with esophageal cancer and 151 healthy controls were included. In univariate Cox regression analysis, PDHA1 was selected from 19 cuproptosis-related genes to be associated with the prognosis of esophageal cancer. At the same time, PDHA1 expression was also different in tumor tissues and healthy tissues. The results of the present study suggested that the PDHA1 expression was relatively low in cancer tissues, which was consistent with previous studies[16,17]. Li et al[17] also reported that high expression of PDHA1 in ovarian cancer cells was significantly correlated with better OS and PFS. Our results suggested that the high expression of PDHA1 expression was associated with poor OS, but not with PFS.

Li et al[18] showed that PDHA1 knockout inhibited glucose entering the tricarboxylic acid cycle, resulting in the reconnection of glutamine metabolism by increasing glutaminase 1 and glutamate dehydrogenase 1 expression, thus increasing the survival rate of glutamine-dependent cells. Consistent with our results of functional enrichment, these DEGs were significantly enriched in tricarboxylic acid cyclic-related pathways. We hypothesized that PDHA1 gene knockdown or low expression causes mitochondrial malfunction, resulting in aberrant generation of intracellular reactive oxygen species (ROS) and adenosine triphosphate, which is compatible with the cuproptosis mechanism.

TNM staging is linked to OS and PFS in esophageal cancer, whereas pathological types are associated with OS in esophageal cancer, according to our findings. Through univariate and multivariate Cox analyses, Zhang et al[19] also claimed that pathological stage is an independent risk factor for OS and PFS in patients with operable esophageal cancer. Therefore, pathological types and TNM stages were included in the prognostic nomogram. Men and women had different rates of esophageal cancer, with men having nearly twice as many cases as women[20]. We did not include sex in our nomogram because, while our analysis indicated that sex impacted PFS, the OS of esophageal cancer was not connected to sex.

CIBERSORT and ssGSEA techniques were used to analyze the composition of tumor-infiltrating immune cells (TICs) in each sample. Regulatory T cells and CD8⁺ T cells have been shown to play an important role in anti-tumor immunity in previous studies[21-23]. Studies have demonstrated that PDHA1 mediates metabolic reprogramming in macrophages[24,25]. However, there was no significant difference between the 22 types of immune cells in the high and low PDHA1 expression groups, according to our findings.

Immune checkpoint inhibitors have recently been studied in a variety of cancers, and they provide a novel therapy option[26]. However, no evidence of a link between cuproptosis and immune checkpoint genes has been found. TNFSF15, TNFRSF14, HHLA2, LGALS9, CD44, TNFRSF18, PDCD1LG2, CD86, and CD40 all have a strong relationship with PDHA1 expression in our recent research. The link between immune checkpoint-related genes expression and prognosis in certain esophageal cancer patients is still debated[27,28]. The immune checkpoint-related genes associated with PDHA1 expression were analyzed using the log-rank test, and the results revealed that HHLA2, TNFRSF18, and CD44 were substantially correlated with prognosis.

Our study had a number of limitations. First, the sample size must be increased to analyze EAC and ESCC individually. Second, if the sample size is high enough, the treatment procedures and stages of the esophagus must be unified. Finally, given that prognostic characteristics were generated and analyzed using data from public databases, further biological evidence, in addition to the statistical evidence we present, was required.

CONCLUSION

This study analyzed the association between cuproprosis-related genes and the prognosis of esophageal cancer in a systematic method. Cuproprosis-related genes, especially PDHA1, can be used as prognostic predictors in esophageal cancer patients, providing additional information on how to treat the disease.

ARTICLE HIGHLIGHTS

Research background

Despite many breakthroughs in treatment, the general prognosis for esophageal cancer, one of the least responsive malignancies to cancer therapy, remains poor. As a result, identifying biomarkers and understanding the molecular mechanisms of esophageal cancer were critical for improving patient outcomes.

Research motivation

A nomogram for predicting the prognosis of esophageal cancer would be developed by evaluating cuproprosis-related genes features and their correlation with prognosis in order to predict the prognosis of esophageal cancer.

Research objectives

Considering cuproprosis-related genes expression was linked to patient prognosis, we intended to develop a nomogram to predict prognosis based on cuproprosis-related genes characteristics and evaluate its prediction performance.

Research methods

Cuproprosis-related genes were found to be linked with esophageal cancer prognosis using univariate COX regression analysis on 151 esophageal cancer samples. The C-index, calibration curve, and receiver operator characteristic (ROC) curve were used to evaluate the prediction ability of a prognostic nomogram created by combining clinicopathological variables and cuproprosis-related genes.

Research results

Univariate COX regression analysis of 19 Cuproprosis-related genes revealed that the expression of pyruvate dehydrogenase A1 (PDHA1) was associated with the prognosis of esophageal cancer. The low PDHA1 expression group had a better prognosis of esophageal cancer, according to the log-rank test. There was no statistical correlation between PDHA1 expression and 22 immune cells; however there was a correlation between PDHA1 expression and several immune checkpoint genes. The C-index, calibration curve, and ROC curve were used to confirm the predictive ability of the esophageal cancer prognostic nomogram, which was developed by combining pathological type, tumor-node-metastasis stage, and PDHA1 expression.

Research conclusions

Cuproprosis-related genes are correlated to esophageal cancer prognosis, and a deep understanding of its molecular mechanism might contribute to novel cancer treatments in the future.

Research perspectives

To enhance the overall survival of esophageal cancer patients, researchers must investigate cuproprosis biomarkers and anticipate possible therapy targets.

ACKNOWLEDGEMENTS

We appreciate the analytical data provided by The Cancer Genome Atlas database (https://cancergenome.nih.gov/) and Genotype-Tissue Expression (https://xenabrowser.net/).