A novel cuproptosis-related prognostic signature and potential value in HCC immunotherapy

Introduction

Hepatocellular carcinoma (HCC) has high heterogeneity and poor prognosis (Dong et al., 2020) and led to the second-highest mortality rate among all cancers (Sung et al., 2021). The 5-year survival rate of HCC patients is less than 20% despite the continuous emergence of immunotherapy and targeted drugs (Llovet et al., 2021). Due to its high heterogeneity, the prognosis of HCC also varies greatly (Pinero et al., 2020). Alpha-fetoprotein (AFP), AFP-L3, glypican-3 (GPC3) and des-gamma-carboxy prothrombin (DCP) are important markers for the diagnosis and prognostication of HCC, but the sensitivity and specificity are limited (Cai et al., 2019; Choi et al., 2019). Therefore, finding novel regulatory mechanisms of HCC and identifying more reliable prognostic markers are crucial for improving the survival rate and promoting precision therapy.

Programmed cell death (PCD) is a process in which cells die after being stimulated by an external signal, which may be physiological or pathological (Hotchkiss et al., 2009). Apoptosis, autophagy and necrosis are three major forms of PCD that have been extensively studied (D'Arcy, 2019). In recent years, new PCD patterns including ferroptosis, pyroptosis, and necroptosis have been discovered, which play an important role in the occurrence and development of tumors (Frank and Vince, 2019; Hirschhorn and Stockwell, 2019; Tsvetkov et al., 2022). Copper ionophore-induced cell death, also called cuproptosis, is the latest proposed PCD pattern. The main mechanism is that copper binds to lipoylated components in the TCA cycle under overload conditions, thus inducing cell death (Tsvetkov et al., 2022). Copper metabolism plays an important role in the regulation of tumor microenvironment, however, the underlying mechanism of cuproptosis in HCC remains unclear. Therefore, it is crucially important to explore the association between cuproptosis and HCC, and to find new strategies for diagnosis and treatment.

In this study, we first preformed multi-genetic analysis of cuproptosis-related genes and further identified three cuproptosis-subgroups by nonnegative matrix factorization (NMF) clustering analysis. Differential pathway enrichment and immune cell infiltration analyses of these three subgroups were conducted to explore the biological mechanisms underlying the three subtypes. The cuproptosis-related prognostic signature (Cu-PS) was established based on the five identified prognostic core genes. The Cu-PS was significantly associated with TNM stage, tumor mutational burden (TMB), drug sensitivity, and immunotherapy response. In conclusion, our study explored the regulatory mechanism of cuproptosis in HCC and provides new targets and strategies for the clinical diagnosis and treatment of HCC patients.

Materials and methods

Data collection and preprocessing

A total of 716 HCC samples with complete survival and mRNA expression information from three HCC cohorts (TCGA-LIHC cohort, n = 370; GSE76427 cohort, n = 115; ICGC-LIRI cohort, n = 231) were included in this study. The clinical information and copy number variation data were downloaded from https://www.cancer.gov/. The mRNA expression and somatic mutation data of TCGA-LIHC were downloaded from https://xenabrowser.net. The clinical and mRNA expression data of the GSE76427 and ICGC-LIRI cohorts were downloaded from https://www.ncbi.nlm.nih.gov/geo/ and https://dcc.icgc.org/projects/, respectively.

The RNA data (FPKM or count format) were transformed into TPM format. The batch effects among the TCGA-LIHC, GSE76427 and ICGC-LIRI cohorts were eliminated via the SVA” package (Yang et al., 2020). The CNV diagram of 16 cuproptosis-related genes was generated by the “Rcircos” package.

Identification of HCC cuproptosis subtypes (Cu-clusters) based on NMF unsupervised clustering analysis

The 16 cuproptosis-related genes were obtained from published literatures, including 13 genes with positive relationships (FDX1, LIPT1, LIAS, DLD, DBT, GCSH, DLST, DLAT, PDHA1, PDHB, SLC31A1, ATP7A and ATP7B) and 3 genes with negative relationships (CDKN2A, MTF1 and GLS). Based on the expression of these 16 cuproptosis-related genes, we performed an unsupervised clustering analysis (NMF) and identified three different cuproptosis-related subtypes. The clustering analysis was performed with the “consensus-cluster plus” package and iterated 1,000 times (Wilkerson and Hayes, 2010).

Pathway enrichment analysis of hallmark and GO gene sets

The hallmark gene sets were downloaded from the MSigDB (http://www.gsea-msigdb.org). Single-sample enrichment scores for the hallmark gene set were estimated with the “GSVA” package (Hanzelmann et al., 2013). The package “clusterProfiler” was applied to annotate the GO functions of DEGs.

Estimation of immune cell infiltration with the CIBERSORT and MCPcounter algorithms

The “MCPcounter” package and CIBERSORT algorithm were employed to estimate the fractions of tumor-infiltrating immune cell subsets based on mRNA transcriptome profiles (Becht et al., 2016). The CIBERSORT algorithm (https://cibersort.stanford.edu/) can be used to evaluate the abundance of 22 immune cells, and the MCPcounter algorithm can be used to quantify the fractions of 8 tumor-infiltrating immune cell types, as well as endothelial cells and fibroblasts.

Identification of cuproptosis subtype-related DEGs

The “limma” package was applied to screen DEGs. A total of 140 DEGs were identified with significance criteria (adjusted p value < 0.001, |logFC>1|). The “heatmap” package was used to present the expression landscape of DEGs among three different cuproptosis subtypes.

Construction of a prognostic signature based on the cuproptosis-related DEGs

The “caret” package was used to divide the patients (from the TCGA-LIHC and GSE76427 cohorts) into training (70%) and testing (30%) cohorts. The independent cohort ICGC-LIRI was used for validation. The “survival” and “glmnet” packages were used for the univariate Cox and LASSO Cox analyses to identify the cuproptosis-related and prognostically significant hub genes. The cuproptosis-related prognostic signature (termed Cu-PS) was constructed by multivariate Cox regression. The formula is as follows: Cu-PS = ${\displaystyle \sum }_i^{n}coe{f}_i*mRN{A}_i$. The selection of the optimal cutoff value was based on the “surv_cutpoint” function in the “survival” package. The “survivalROC” and “survminer” packages were used to generate the receiver operating characteristic (ROC) and Kaplan–Meier (K−M) curves, respectively.

Quantitative real-time PCR (qRT–PCR) in cell lines

To validate the expression levels of the 5 hub genes in normal and HCC cell lines, we cultured two HCC cell lines (Huh7 and HLE) and one human hepatocellular cell line (MIHA). Total RNA was extracted from the above 3 cell lines, and cDNA was synthesized with a reverse transcription kit (TaKaRa, Japan). The qRT–PCR was performed using SYBR Green Mix (TaKaRa, Japan) and a C1000 system (Bio–Rad, Hercules, CA). The primer sequences of the 5 hub genes are listed in Table 1. The RNA quality was assessed, and RNA levels were normalized based on human GAPDH.

TABLE 1

TABLE 1. The primer sequences of the 5 hub genes.

Ability of the Cu-PS to predict immunotherapy response

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was employed to quantify immunosuppressive and dysfunctional factors in the tumor immune microenvironment (TIME) and to estimate the ability of cancer cells to escape antitumor immunity (Jiang et al., 2018). In addition, the IMvigor210 cohort (Mariathasan et al., 2018) (348 urothelial carcinoma patients) and Liu et al. cohorts (Liu et al., 2019) (121 melanoma patients) were employed to validate the relationship between the Cu-PS and immunotherapy response. The mRNA data of the two immunotherapy cohorts were transformed into TPM values before further analysis.

Analysis of the correlation of the Cu-PS with drug sensitivity

The “pRRophetic” package (Geeleher et al., 2014) was used to calculate the IC50 values of drugs in the Genomics of Drug Sensitivity in Cancer (GDSC) database. We analyzed the correlation between IC50 values and the Cu-PS with Spearman correlation analysis to explore the association of the Cu-PS with drug sensitivity. |Cor| > 0.2 and adjusted p < 0.05 were used as cutoffs for identifying significant correlations.

Statistical analyses

All statistical analyses were performed with the software of R-4.0.2. The Wilcoxon rank-sum test and Kruskal–Wallis tests were applied to identify differences between two and among three groups, respectively (Hazra and Gogtay, 2016). K−M analysis and the log-rank test were utilized to analyze differences between distinct Cu-clusters, Cg-clusters and Cu-PS subgroups. The mutation waterfall plot was generated with the “maftools” package (Mayakonda et al., 2018). The CNV circle graph of 16 cuproptosis-related genes in human chromosomes was generated with the “RCircos” package (Zhang et al., 2013). All tests were bilateral, p < 0.05 was considered significant, and the false discovery rate (FDR) was calculated for multiple hypothesis testing (Ferreira, 2007).

Results

Multiomics landscape of 16 cuproptosis-related genes and identification of cuproptosis subtypes in HCC

The study analysis flowchart was shown in Supplementary Figure S1. First, we analyzed the differences in the mRNA expression of 16 cuproptosis genes between HCC and normal liver tissues. The results showed that in HCC patients, 9 of 13 genes positively associated with cuproptosis were expressed at lower levels, 3 of 13 cuproptosis positive genes (LIPT1, DLAT and ATP7A) and the 3 cuproptosis negative genes (GLS, MTF1 and CDKN2A) were upregulated, indicating that the level of cuproptosis was lower in HCC than in normal liver tissues (Figure 1A). The mutation analysis showed that 18 of 361 (4.99%) HCC patients had mutations in cuproptosis-related genes (Figure 1B). The copy number alteration (CNA) frequency analysis showed that most of the cuproptosis genes had deletions affecting copy number. ATP7B had the highest frequency of deletion, while LIAS had the highest frequency of amplification (Figure 1C). The chromosomal locations of the 16 cuproptosis genes CNA are shown in Figure 1D.

FIGURE 1

FIGURE 1. Genetic variation landscape and Unsupervised clustering of 16 Cu-RGs in HCC. (A) Comparison of 16 Cuproptosis-related genes (Cu-RGs) between normal and HCC tissues in mRNA expression. (B) The mutation frequency and classification of 16 Cu-RGs in TCGA-LIHC cohort. (C) The CNV variation frequency of 16 Cu-RGs in TCGA-LIHC cohort. (D) The location of CNV alteration of 16 Cu-RGs on the chromosomes in TCGA-LIHC. (E) Co-expression and prognosis relationship of 16 Cu-RGs in TCGA-LIHC. The regulation of Cu-RGs to cuproptosis were depicted by circles (lift) in two colors: red, cuproptosis-down; gray, cuproptosis-up. The lines connecting 16 Cu-RGs represented their interaction with each other. The size of each circle represented the prognosis effect of each regulator and scaled by p-value. The color on the right half of the circle represents the effect of Cu-RGs on prognosis: green, protective factors; purple, risk factors. (F) Kaplan-Meier curves of overall survival (OS) for 485 HCC patients in TCGA-GEO cohort with different Cuproptosis-clusters (termed as C1, C2 and C3). The numbers of patients in C1, C2, and C3 are 115, 172, and 198, respectively.

TCGA-LIHC and GSE76427 samples with available survival information were employed to generate the comprehensive crosstalk network of the 16 cuproptosis genes (Figure 1E). We found that the three genes negatively associated with cuproptosis (CDKN2A, MTF1 and GLS) were risk factors for HCC overall survival (OS) (HR < 1, univariate Cox test), and most of the genes positively associated with cuproptosis were favorable factors (HR < 1, univariate Cox test), which indicated that HCC patients may benefit from cuproptosis. The NMF algorithm was used to classify 485 HCC samples into three distinct cuproptosis subgroups, termed Cu-cluster 1, Cu-cluster 2 and Cu-cluster 3, based on the expression of the 16 cuproptosis genes. The process of cluster analysis (rank = 2:7) is shown in Supplementary Figure S2. K−M survival analysis showed that the patients in Cu-cluster 3 had the best OS benefit, followed by those in Cu-cluster 2 and Cu-cluster 1 (p = 0.039, Figure 1F). We performed a multi-genetic analysis and identified three distinct cuproptosis-subgroups associated with different OS prognoses.

Analysis of enriched pathways and immune cell infiltration among the three Cu-clusters

We explore the cuproptosis level based on the expression of the 16 cuproptosis-related genes among the three Cu-clusters. The patients in Cu-cluster 1 had the lowest cuproptosis level, as indicated by lower expression of genes positively associated with cuproptosis and overexpression of genes negatively associated with cuproptosis. The patients in Cu-cluster 2 had a moderate cuproptosis level, and those in Cu-cluster 3 had the highest cuproptosis level (Figure 2A). To discover the underlying biological mechanisms behind the differences in survival between the three Cu-clusters, we performed gene set enrichment analysis (GSEA) and immune cell infiltration analysis. The GSEA results showed that Cu-cluster 1 was enriched in mTORC1 signaling, E2F targets and G2/M checkpoint pathways. Cu-cluster 3 was associated with metabolic pathway terms, such as bile acid metabolism, peroxisome, lipogenesis and fatty acid metabolism. Cu-cluster 2 was enriched in the terms Notch signaling, angiogenesis, epithelial-mesenchymal transition and TGF-β signaling (Figure 2B). We further analyzed immune cell infiltration with MCPcounter (Figure 2C) and the CIBERSORT algorithm (Figure 2D). The HCC patients in Cu-cluster 1 had the highest levels of inhibitory immune cells (myeloid dendritic cells (DCs), regulatory T cells (Tregs) and M0 macrophages), while those in Cu-cluster 2 had the highest levels of stromal cell subsets (endothelial cells and fibroblasts); those in Cu-cluster 3 had lower immune cell infiltration. Furthermore, the patients in cluster-C3 and normal samples were less infiltrated in T cell, B cell as well as endothelial cells and fibroblast (Supplementary Figure S3). These results indicated that cuproptosis subtype is associated with tumor microenvironment factors.

FIGURE 2

FIGURE 2. Biological characteristics and immune cell infiltration characteristics in distinct cuproptosis-clusters. (A) Expression heatmap of 16 Cu-RGs in the TCGA-GEO cohort. The patient annotations of Cuproptosis-cluster, OS-status, gender, TNM stage, age and cohorts were presented at the top of heatmap. (B) GSVA score heatmap of Hallmark pathways among three Cuproptosis-clusters. (C,D) Comparison of immune cell infiltration among three Cuproptosis-clusters, which was estimated by MCPcounter (C) and CIBERSORT (D) algorithm, respectively.

Identification of cuproptosis-related and prognosis-related hub genes in HCC

To explore potential biological behaviors, we identified 140 DEGs among the three subtypes (|logFC|>0.5, adjusted p < 0.01, Figure 3A) and further performed GO enrichment analysis. The results showed that the DEGs were mainly enriched in metabolic processes, such as steroid metabolic processes, xenobiotic metabolic processes, fatty acid biosynthetic processes, cellular hormone metabolic processes, estrogen metabolic processes, vitamin D metabolic processes and bile acid and bile salt transport (Figure 3B). In the training cohort, univariate, LASSO and multivariate Cox analyses were applied to identify the cuproptosis-related and prognosis-related hub genes. A total of 10 prognostic DEGs were selected through LASSO Cox analysis (Figures 3C,D), and 5 prognostic hub DEGs were identified via multivariate Cox analysis. The hazard ratios and p values of these selected genes are shown in the forest plot of the univariate Cox analysis results (Figures 3E,F). One hub DEG (C7) was the favorable factor for HCC prognosis, while the other hub DEGs (MAGEA6, HK2, CYP26B1 and EPO) were risk factors.

FIGURE 3

FIGURE 3. Screening of Cu-DEGs and construction of cuproptosis-related prognostic signatures (Cu-PS) (A) The Venn diagram illustrated the 140 Cu-DEGs among three cuproptosis-clusters. (B) GO enrichment analysis revealed the biological characteristics of the 140 DEGs. (C) 10 Cu-DEGs screened by LASSO regression. (D) 10-fold cross-validation plot of LASSO regression. (E) Univariate Cox forest-plot of the 10 selected Cu-DEGs. (F) Forest-plot of 5 hub Cu-DEGs in Cu-PS identified by multivariate Cox regression analysis.

Construction and validation of the cuproptosis-related prognostic signature in HCC

With the training cohort, we further constructed a cuproptosis-related prognostic signature termed the Cu-PS via multivariate Cox analysis. The formula used to calculate the Cu-PS is described in the Methods section. The Cu-PS the training cohort was calculated with the same formula. The patients with an increased Cu-PS had a high fraction of death and shortened survival time in both the training cohort and testing cohort (Figures 4A–D). As can be seen in the heatmaps, the risk genes (MAGEA6, HK2, CYP26B1 and EPO) were upregulated with increasing Cu-PS, while the protective gene C7 was downregulated (Figures 4E,F). HCC samples were divided into high- and low-risk groups based on the median Cu-PS. The HCC patients in the training cohort with low Cu-PS had longer OS, and a similar result was found in the testing cohort (Figures 4G,I; p < 0.001). The AUCs of Cu-PS for predicting 1-, 2-, and 3-year survival were 0.697, 0.704, and 0.682, respectively, in the training cohort and 0.713, 0.656, and 0.644, respectively, in the testing cohort (Figures 4H,J).

FIGURE 4

FIGURE 4. Construction and validation of Cu-PS in training and testing cohorts. (A,B) The range of Cu-PS in the training and testing cohorts, and patients were cutoff into high- and low-two subgroups by the and median value of Cu-PS. (C,D) The distribution plots showed the survival status of patients with increasing Cu-PS in the training and testing cohorts. (E,F) The heatmaps showed the expression of 5 hub Cu-DEGs between two Cu-PS subgroups in the training and testing cohorts. (G,H) The Kaplan-Meier curves of overall survival (OS) between the two subgroups in the training and testing cohorts. (I,J) The ROC curves of the Cu-PS in predicting 1-, 2-, and 3-year OS in the training and testing cohorts.

To further confirm the prognostic value of the 5 hub genes (C7, MAGEA6, HK2, CYP26B1 and EPO) in HCC, we performed K-M survival analysis in the training cohort, testing cohort and independent ICGC-LIRI cohort (Figures 5A–M). The optimal cutoff value was identified by the “cutpoint” function in the “survival” package. Consistent with the previous results (Figure 3F), the patients with high expression of MAGEA6, HK2, CYP26B1 and EPO had a worse prognosis in the training, testing and ICGC-LIRI cohorts (Figures 5D,E,G–M and G-M, all p < 0.05; Figure 5F, p = 0.105). Meanwhile, patients with high C7 expression lead a prognostic benefit in the training cohort and ICGC-LIRI cohort (Figures 5A,C; both p < 0.001), and the same trend was found in the testing cohort (Figure 5B, p = 0.322).

FIGURE 5

FIGURE 5. The prognosis value of 5 hub Cu-DEGs in three HCC cohorts. (A–M) The Kaplan-Meier curves of C7 in training (A), testing (B) and ICGC (C) cohorts, as well as CYP26B1(D–F), EPO (G–I), HK2 (H–J) and MAGEA6 (K–M).

We further validated the 5 hub genes in tissues (TCGA-LIHC dataset) and cultured normal (MIHA) and HCC cell lines (Huh7 and HLE). The results from the TCGA-LIHC dataset showed that C7, CYP26B1 and EPO were significantly downregulated in HCC, while MAGEA6 was significantly upregulated (Supplementary Figure S4A). Similar results were found in the cell lines, but EPO and CYP26B1 did not show the same differences (Supplementary Figure S4B). The expression of EPO in HCC cell lines tended to be higher (without a significant difference) than that in the normal cell line. The expression of CYP26B1 in MIHA and Huh7 cell lines was significantly higher than that in the HLE cell line.

Analysis of the correlation of the Cu-PS with clinical characteristics and TMB

Patients with TNM stage Ⅱ and Ⅲ/Ⅳ disease had a higher Cu-PS than those with stage Ⅰ disease (Figure 6A, both p < 0.01). In addition, patients who died had a higher Cu-PS than those who survived (Figure 6B, p < 0.001). Moreover, the TMB value in the high-Cu-PS group was significantly higher than that in the low-Cu-PS group (Figure 6C, p = 0.013). Waterfall plots of the top 20 frequently mutated genes are shown in Figure 6D. The patients with high Cu-PS had a significantly higher mutation frequency of TP53 but a lower frequency of AXIN1 mutation. As expected, the patients with TP53 mutation had a significant OS benefit compared with TP53 wild-type patients (Figure 6E, p = 0.006). We further performed GO enrichment analysis of the DEGs between the high- and low-Cu-PS groups. The results showed that the DEGs were mainly enriched in the cell cycle and metabolic processes, such as nuclear division, mitotic cell cycle phase transition, chromosome segregation, organic acid biosynthetic processes, carboxylic acid biosynthetic processes, and hormone metabolic processes (Figure 6F). These findings confirmed the prognosis values and revealed the underlying mechanism of Cu-PS in HCC.

FIGURE 6

FIGURE 6. Correlation between the Cu-PS and clinical features and the landscape of somatic mutation. (A,B) Comparison of Cu-PS among different TNM stages (A) and survival outcomes (B). Differences in risk scores between different survival outcome. (C) Comparison of TMB value between high- and low- Cu-PS subgroups. (D) The waterfall plots showed the somatic mutation spectrum of the high- and low-risk groups. (E) Kaplan-Meier curves between TP53 mutation (110 patients) and wild (248 patients) groups in TCGA-LIHC. (F) GO enrichment analysis of the DEGs between high- and low- Cu-PS subgroups.

Ability of the Cu-PS to predict drug sensitivity and immunotherapy efficacy

To further investigate the potential application value of the Cu-PS in HCC treatment, we explored the correlations with drug sensitivity and immunotherapy efficacy based on the GDSC database and two immunotherapy cohorts (the GSE78220 and IMvigor210 cohorts). We identified 16 drugs in the GDSC database that were significantly associated with the Cu-PS by Spearman correlation analysis (Figure 7A, p < 0.05). Among them, 10 drugs showed a correlation between drug sensitivity and the Cu-PS, including sorafenib, the Src inhibitor A.770,041, the RAF inhibitor AZ628, the Src/Abl dual-kinase inhibitor AZD.0530, and the JNK inhibitor AS601245 (all cor>0.2, all p < 0.05). Six drugs showed a correlation between the Cu-PS and drug resistance, including the PARP inhibitor ABT-888, all-trans retinoic acid (ATRA), the Bcl-2 inhibitor ABT.263, the Akt inhibitor A.443654 and the AMPK activator (all cor<0.2, all p < 0.05). Sorafenib is a first-line treatment for advanced HCC. As expected, the IC50 of sorafenib was significantly positively associated with the Cu-PS (Figure 7B, cor = 0.25, p < 0.001), which was in line with the sensitivity of HCC to sorafenib. These results imply that the Cu-PS is correlated with drug sensitivity and might be a potential biomarker for guiding drug treatment selection.

FIGURE 7

FIGURE 7. Efficacy prediction of chemotherapy drugs and immunotherapy. (A) The correlation of Cu-PS with IC50 of drugs in GDSC database. (B) The correlation scatter plot of Cu-PS and Sorafenib-IC50. (C) The proportion of immune cells in TME between the high- and low- Cu-PS subgroups. (D) Comparison of TIDE score between the two Cu-PS subgroups. (E) Kaplan-Meier curves for high- and low- Cu-PS subgroups in the anti-PD-L1 immunotherapy cohort (IMvigor210). (F) The proportion of immune response (CR, PR, SD and PD) between two Cu-PS subgroups in IMvigor210 cohort. (G) Kaplan-Meier curves for high- and low- Cu-PS subgroups in the anti-PD1 immunotherapy cohort (GSE78220 cohort) (H) The proportion of immune response between two Cu-PS subgroups in GSE78220 cohort.

Immunotherapy has made major breakthroughs in the treatment of liver cancer. The immune cell infiltration analysis showed that the patients with a high Cu-PS had a suppressive immune microenvironment, poor prognosis, and significantly enriched levels of Tregs, M0 macrophages, and neutrophils but decreased levels of B cells, CD4⁺ T cells and mast cells (Figure 7C). In addition, the patients in low-Cu-PS group had a high fraction of memory CD8 T cell, while the effector CD8 T cell was not significant (Supplementary Figure S5). The TIDE scores in the high-Cu-PS group were also significantly higher than those in the low-Cu-PS group, which was associated with resistance to immunotherapy (Figure 7D). We further confirmed the predictive ability of the Cu-PS in two immunotherapy cohorts (the IMvigor210 and Liu cohorts). The results showed that the patients in the Cu-PS-low group in the IMvigor210 cohort had better OS (Figure 7E, p < 0.001). The fractions of patients who achieved complete response (CR) and partial response (PR) in the Cu-PS-low group were higher than those in the Cu-PS-high group (Figure 7F). Similar results were found in the Liu cohorts (Figures 7G,H). In summary, these findings suggest that application of the Cu-PS might improve drug selection and immunotherapy response prediction in HCC.

Discussion

Cuproptosis is a newly discovered PCD process, and copper metabolism plays an important role in tumor development, invasion and metastasis. In this study, we systematically performed multiomics bioinformatics analyses to explore the association of cuproptosis-related genes with genomic and TIME characteristics, prognosis and immunotherapy response in HCC.

In this study, we assessed the relevance of cuproptosis and the immune microenvironment in HCC. We performed multiomics analysis of the 16 cuproptosis-related genes and found that the level of cuproptosis was significantly higher in normal liver tissues than in HCC tissues, which indicated that cuproptosis may suppress tumorigenesis to a certain extent. We further identified three distinct cuproptosis-related subgroups (Cu-clusters) associated with OS and with different enriched tumor hallmark pathways. Cu-cluster 1 was enriched in mTORC1 signaling, E2F targets and G2/M checkpoint pathways. Cu-cluster 3 was enriched in metabolism pathways, while Cu-cluster 2 was enriched in angiogenesis and epithelial-mesenchymal transition pathways. The HCC patients in Cu-cluster 1 had the highest levels of inhibitory immune cells (myeloid DCs, Tregs and M0-macrophages), while Cu-cluster 2 had the highest levels of stromal cell subsets (endothelial cells and fibroblasts), and Cu-cluster 3 had lower immune cell infiltration.

We further constructed a 5-gene (C7, MAGEA6, HK2, CYP26B1 and EPO) prognostic signature termed the Cu-PS via univariate Cox, LASSO and multivariate Cox regression analyses. Complement component 7 (C7) is an essential component of the complement system (CS) and is involved in membrane attack complex (MAC) formation. Moreover, the C7 peptide can inhibit Akt and Erk1/2 and further suppress HCC cell migration and invasion induced by HGF (Zhao et al., 2019).

Melanoma-associated antigen family A (MAGEA) antigens are expressed in a variety of malignancies. MAGEA6 can promote pancreatic (Tsang et al., 2020), lung cancer (Pineda et al., 2015) and colorectal cancer (Wu et al., 2018) carcinogenesis by inhibiting autophagy. In addition, MAGEA6 regulates stemness and self-renewal in HCC by activating the AMPK signaling pathway and is associated with poor prognosis (Guo et al., 2019).

Pericyte-hexokinase 2 (HK2) is the key rate-limiting enzyme of the glycolytic pathway, which is associated with pathological stage and prognosis (DeWaal et al., 2018). In the mouse model of HCC, HK2 silencing could synergistically enhance the sensitivity of HCC cells to sorafenib (Yu et al., 2022). CYP26 enzymes are the major enzymes responsible for the clearance of retinoic acid (RA), which is an important endogenous signaling molecule that regulates the cell cycle and maintains epithelial cells. CYP26 inhibitors can enhance endogenous RA activity in a cell-type-specific manner and might be new, attractive targets in cancer and skin disease treatment (Tsvetkov et al., 2022). Erythropoietin (EPO) is primarily synthesized in the kidney and can promote erythrocyte production. Under hypoxic conditions, EPO is upregulated to promote the production of red blood cells and enhance the oxygen-carrying capacity (Liu et al., 2020).

Although we conducted a comprehensive and systematic analysis, the identification of additional new cuproptosis-related genes will enrich our research. Due to the lack of an HCC immunotherapy cohort, we assessed the ability of a prognostic signature (Cu-PS) to predict the efficacy of immunotherapy in two cohorts of patients with different cancers treated with immunotherapy. In conclusion, we performed a comprehensive analysis of 16 cuproptosis-related genes in 716 HCC samples and identified a novel HCC prognostic signature (Cu-PS), providing a new strategy for predicting HCC prognosis and immunotherapy efficacy.

Conclusion

Compared with already known biomarkers such as AFP, GPC3 and DCP, our study constructed a novel cuproptosis-related prognostic signature (Cu-PS) that might be a useful biomarker for predicting immunotherapy response and enhancing diagnosis and treatment of HCC, which indicates that cuproptosis is associated with the TIME and HCC prognosis.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

H-MZ had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis. X-XW, L-HW and HJ contributed equally to this work. Concept and design, X-XW and H-MZ; Acquisition and analysis of data, X-XW, L-HW, HJ, Q-QL, S-ZD, and Q-YD; Drafting of the manuscript, X-XW and S-ZD; Statistical analysis, X-XW and Q-QL. All authors have read and agreed to the published version of the manuscript.

Funding

This study was supported by Natural Science Basic Research Program of Shaanxi (No. 2021JZ-35) and National Natural Science Foundation of China (No.81572699).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2022.1001788/full#supplementary-material

References

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17 (1), 218. doi:10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, J., Chen, L., Zhang, Z., Zhang, X., Lu, X., Liu, W., et al. (2019). Genome-wide mapping of 5-hydroxymethylcytosines in circulating cell-free DNA as a non-invasive approach for early detection of hepatocellular carcinoma. Gut 68 (12), 2195–2205. doi:10.1136/gutjnl-2019-318882

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, J., Kim, G. A., Han, S., Lee, W., Chun, S., and Lim, Y. S. (2019). Longitudinal assessment of three serum biomarkers to detect very early-stage hepatocellular carcinoma. Hepatology 69 (5), 1983–1994. doi:10.1002/hep.30233

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Arcy, M. S. (2019). Cell death: A review of the major forms of apoptosis, necrosis and autophagy. Cell Biol. Int. 43 (6), 582–592. doi:10.1002/cbin.11137

PubMed Abstract | CrossRef Full Text | Google Scholar

DeWaal, D., Nogueira, V., Terry, A. R., Patra, K. C., Jeon, S. M., Guzman, G., et al. (2018). Author Correction: Hexokinase-2 depletion inhibits glycolysis and induces oxidative phosphorylation in hepatocellular carcinoma and sensitizes to metformin. Nat. Commun. 9 (1), 2539. doi:10.1038/s41467-018-04182-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, L. Q., Peng, L. H., Ma, L. J., Liu, D. B., Zhang, S., Luo, S. Z., et al. (2020). Heterogeneous immunogenomic features and distinct escape mechanisms in multifocal hepatocellular carcinoma. J. Hepatol. 72 (5), 896–908. doi:10.1016/j.jhep.2019.12.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, J. A. (2007). The Benjamini-Hochberg method in the case of discrete test statistics. Int. J. Biostat. 3 (1), 11. doi:10.2202/1557-4679.1065

PubMed Abstract | CrossRef Full Text | Google Scholar

Frank, D., and Vince, J. E. (2019). Pyroptosis versus necroptosis: Similarities, differences, and crosstalk. Cell Death Differ. 26 (1), 99–114. doi:10.1038/s41418-018-0212-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9 (9), e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J. C., Yang, Y. J., Zhang, J. Q., Guo, M., Xiang, L., Yu, S. F., et al. (2019). microRNA-448 inhibits stemness maintenance and self-renewal of hepatocellular carcinoma stem cells through the MAGEA6-mediated AMPK signaling pathway. J. Cell. Physiol. 234 (12), 23461–23474. doi:10.1002/jcp.28915

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanzelmann, S., Castelo, R., and Guinney, J. (2013). Gsva: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hazra, A., and Gogtay, N. (2016). Biostatistics series module 3: Comparing groups: Numerical variables. Indian J. dermatol. 61 (3), 251–260. doi:10.4103/0019-5154.182416

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirschhorn, T., and Stockwell, B. R. (2019). The development of the concept of ferroptosis. Free Radic. Biol. Med. 133, 130–143. doi:10.1016/j.freeradbiomed.2018.09.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Hotchkiss, R. S., Strasser, A., McDunn, J. E., and Swanson, P. E. (2009). Cell death. N. Engl. J. Med. 361 (16), 1570–1583. doi:10.1056/NEJMra0901217

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, D., Schilling, B., Liu, D., Sucker, A., Livingstone, E., Jerby-Arnon, L., et al. (2019). Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat. Med. 25 (12), 1916–1927. doi:10.1038/s41591-019-0654-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, W., Varier, K. M., Sample, K. M., Zacksenhaus, E., Gajendran, B., and Ben-David, Y. (2020). Erythropoietin signaling in the microenvironment of tumors and healthy tissues. Adv. Exp. Med. Biol. 1223, 17–30. doi:10.1007/978-3-030-35582-1_2

PubMed Abstract | CrossRef Full Text | Google Scholar

Llovet, J. M., De Baere, T., Kulik, L., Haber, P. K., Greten, T. F., Meyer, T., et al. (2021). Locoregional therapies in the era of molecular and immune treatments for hepatocellular carcinoma. Nat. Rev. Gastroenterol. Hepatol. 18 (5), 293–313. doi:10.1038/s41575-020-00395-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 554 (7693), 544–548. doi:10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Pineda, C. T., Ramanathan, S., Fon Tacer, K., Weon, J. L., Potts, M. B., Ou, Y. H., et al. (2015). Degradation of AMPK by a cancer-specific ubiquitin ligase. Cell 160 (4), 715–728. doi:10.1016/j.cell.2015.01.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinero, F., Dirchwolf, M., and Pessoa, M. G. (2020). Biomarkers in hepatocellular carcinoma: Diagnosis, prognosis and treatment response assessment. Cells 9 (6), 1370. doi:10.3390/cells9061370

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Ca. Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsang, Y. H., Wang, Y., Kong, K., Grzeskowiak, C., Zagorodna, O., Dogruluk, T., et al. (2020). Differential expression of MAGEA6 toggles autophagy to promote pancreatic cancer progression. Elife 9, e48963. doi:10.7554/eLife.48963

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsvetkov, P., Coy, S., Petrova, B., Dreishpoon, M., Verma, A., Abdusamad, M., et al. (2022). Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 375 (6586), 1254–1261. doi:10.1126/science.abf0529

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking. Bioinformatics 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, F., Liu, F., Dong, L., Yang, H., He, X., Li, L., et al. (2018). miR-1273g silences MAGEA3/6 to inhibit human colorectal cancer cell growth via activation of AMPK signaling. Cancer Lett. 435, 1–9. doi:10.1016/j.canlet.2018.07.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Huang, X., Liu, Z., Qin, W., and Wang, C. (2020). Metabolism-associated molecular classification of hepatocellular carcinoma. Mol. Oncol. 14 (4), 896–913. doi:10.1002/1878-0261.12639

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Q., Dai, W., Ji, J., Wu, L., Feng, J., Li, J., et al. (2022). Sodium butyrate inhibits aerobic glycolysis of hepatocellular carcinoma cells via the c-myc/hexokinase 2 pathway. J. Cell. Mol. Med. 26 (10), 3031–3045. doi:10.1111/jcmm.17322

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Meltzer, P., and Davis, S. (2013). RCircos: an R package for Circos 2D track plots. BMC Bioinforma. 14, 244. doi:10.1186/1471-2105-14-244

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, M., Wang, Y., Liu, Y., Zhang, W., Liu, Y., Yang, X., et al. (2019). C7 peptide inhibits hepatocellular carcinoma metastasis by targeting the HGF/c-Met signaling pathway. Cancer Biol. Ther. 20 (12), 1430–1442. doi:10.1080/15384047.2019.1647051

PubMed Abstract | CrossRef Full Text | Google Scholar