The RNA-sequencing data (FPKM value) of 374 HCC patients and the corresponding clinical information were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/, accessed in July 2022). The mRNA expression data and clinical information of ICGC-LIRI-JP dataset were downloaded from International Cancer Genomics Consortium (ICGC) website. Patients without survival information were excluded for further analysis. A total of 368 HCC patients and 50 adjunct non-tumor samples in TCGA and 240 HCC patients in ICGC were involved in the present study. The clinical information of HCC samples was listed in Supplementary Table S1. For data normalization, the values of FPKM data in two cohorts were transformed into transcripts per million kilobase (TPM) values. Furthermore, somatic mutations data and copy number variation (CNV) data were also obtained from TCGA database. This study was based on public databases, the approval of the local ethics committee was not required. Data analysis was performed with the R (version 4.2.1), R-studio and Bioconductor packages.

Human normal liver cells (LO2) and Hepatocellular Carcinoma cells (Huh7 and MHCC97H) were purchased from Shanghai Cell Bank of Chinese Academy of Sciences (China), which were maintained in DMEM plus 10% FBS (Gibco, United states) and 1% penicillin-streptomycin. All these cell lines were incubated at 37 °C under the condition of 5% CO2.

Male C57BL/6 mice of 4-week-old were separated into control and HCC group. For HCC group, 20 μg pT3-EF1a-HA-myr-AKT (Addgene), 20 μg pT/Caggs-NRASV12 (Addgene) along with 2.85 μg SB transposase plasmid (Addgene) were diluted in 2.0 mL saline, and then injected into mice by tail vein injection. For control group, equal volume of saline was injected into mice by tail vein injection. All mice were sacrificed to collect liver and tumor tissue 4–6 weeks later.

A total of 10 Cuproptosis-Related Genes (CRGs) were obtained from prior papers (Tsvetkov et al., 2022) and listed in Supplementary Table S2. The differential expression of CRGs between HCC and normal liver tissues were identified using the “limma” package (Ritchie et al., 2015).

The mutation frequency and waterfall plot of 10 CRGs in HCC patients were generated by the “maftools” package (Mayakonda et al., 2018). The location of 10 CRGs on 23 chromosomes was drawn by the “circlize” package (Gu et al., 2014). The CNV values were set by .3 as a threshold. The Cleveland dot plot was visualized the frequency of CNV by the “ggpubr” package.

Gene Ontology (GO) of cuproptosis-related genes was analyzed via “clusterProfiler” package (Yu et al., 2012). Similarly, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis by the same package. We applied the Benjamini−Hochberg method for the multiple correction, and terms with adjusted p < .05 were significantly enriched. The results of GO and KEGG were shown in Supplementary Tables S3,S4.

We performed the univariate and multivariate Cox regression analysis for overall survival (OS) in HCC clinical data. p-values and hazard ratio (HR) with 95% confidence interval (CI) were generated by Cox proportional hazard regression. CRGs with a significant prognostic value were selected for further analysis. The forest plot of each variables was generated through the “forestplot” R package (Dettori et al., 2021).

We established a predictive scoring model using LASSO-Cox analysis. CRGs-score = (coefficient A) * gene A+ (coefficient B) * gene B+ …… + (coefficient N) * gene N. where coefficient N, gene N represented the coefficient index, and the gene expression level, respectively. Subsequently, we divided patients into the high-risk score and low-risk score subgroups according to the median value of CRGs-score, and the overall survival (OS), disease-free interval (DFI), disease-specific survival (DSS), and progression-free interval (PFI) were compared between the two subgroups via Kaplan-Meier analysis. The predictive accuracy of this model was evaluated by performing time-dependent receiver operating characteristic (ROC) analysis (Kamarudin et al., 2017).

Independent prognostic variables were integrated to construct the nomogram in predicting 1-, 3-, and 5-year overall survival time with stepwise Cox regression analysis. The calibration curves of the CRGs-score and other clinical indicators were estimated using 500 bootstrapping to determine bias-corrected estimates of predicted versus actual values. The curves were plotted to compare nomogram-predicted and observed 1-, 3-, and 5-year survival time.

We performed the CIBERSORT algorithm to quantify the proportion of immune cells between CRGs-score subgroups. We calculated the normalized gene expression data by LM22 signature and 500 permutations. Then we performed pathway enrichment analysis by the “GSVA” package. We downloaded immune associated signatures to carry out GSVA analysis. Moreover, we examined seven important immune checkpoint molecules including CTLA4, HAVCR2, LAG3, PD-L1, PD-1, PDCD1LG2, TIGIT between two subtypes. The TIMER database (cistrome.shinyapps.io/timer) (Li et al., 2017) was used to investigate the correlation between the expression of CRGs and the abundance of six immune cells (B Cells, CD8⁺ T-cell, CD4⁺ T-cell, macrophages, neutrophils, and dendritic cells).

Total RNA was extracted by Trizol reagent (Invitrogen, United states) according to the manufacturer’s protocol. Then, the RNA was reverse transcribed using PrimeScript RT reagent Kit with gDNA Eraser (Takara, Japan). The cDNAs were subjected to SYBR Green dye qPCR analysis. The primers used in real-time qPCR assays were listed in Supplementary Table S5.

The “pRRophetic” package was used to analyze the drug sensitive prediction based on the Genomics of Drug Sensitivity in Cancer (GDSC) database. The half-maximal inhibitory concentrations (IC50) of the samples were calculated by ridge regression algorithm (Geeleher et al., 2014) and the difference of anticancer drug sensitivity were compared between the high CRGs and low CRGs-score group.

The log-rank test and “ggsurvplot” package were utilized to perform the survival analysis. The cut-off value was 2.309 which was the median value of CRGs-score. The risk factors and drug sensitivity analysis between the two groups were performed using the Wilcox test. QPCR data of validation experiments was calculated in GraphPad Prism eight software. The independent Student’s t-test was applied to compare data between two subgroups. All p < .05 represents statistical significance.

The prevalence of Cuproptosis-Related Genes (CRGs) in HCC was determined focusing on somatic mutations and copy number variations (CNVs). At the genetic level, mutation of CRGs was observed in 18 of 364 (4.95%) HCC samples (Figure 1A). Among different classification of mutations, we found that missense mutation was most frequently observed. And single nucleotide polymorphism (SNP) was the most prevalent variant type. C > T and C > A were both ranked as the top single nucleotide variants (SNV) class. Of these, we also found that CDKN2A (3%), MTF1 (1%) and DLD (1%) showed higher mutation frequencies than other genes (Figure 1B). In addition, the location of CNV alterations of these 10 CRGs on chromosomes were shown in Figure 1C. However, the CNV alterations were not universally prevalent among these CRGs. Specifically, LIAS, GLS, and DLD showed distinctive copy number amplification, while CDKN2A, MTF1, DLAT, FDX1, and PDHB exhibited significant CNV deletion (frequency >3%, Figure 1D), indicating that changes in CNVs of CRGs may be important factors leading to abnormal gene expression.

To further clarify the function of CRGs, we further performed enrichment analysis using GO and KEGG database. The 10 CRGs were mainly involved in acetyl-CoA biosynthetic and metabolic process from pyruvate, TCA cycle, thioester biosynthetic process, nucleoside and ribonucleoside bisphosphate biosynthetic process, pyruvate metabolic process, glucose metabolic process, lipoate metabolic process, protein lipoylation, dicarboxylic acid metabolic process, alpha-amino acid catabolic process, dihydrolipoamide metabolic process, senescence-associated heterochromatin focus assembly, cartilage homeostasis, positive regulation of macrophage apoptotic process, protein succinylation in GO analysis (Figure 2A). KEGG pathway analysis also suggested that these CRGs were mainly involved in TCA cycle, glycolysis, pyruvate and carbon metabolism, lipoic acid metabolism, central carbon metabolism in cancer, biosynthesis of cofactors, glucagon signaling pathway, HIF-1 signaling pathway, D-Amino acid metabolism (Figure 2B, data was shown in Supplementary Tables S3,S4 ).

Next, we explored the expression of the 10 CRGs in 368 HCC samples (6 samples without clinical information were excluded for this analysis) and normal liver tissue (50 samples) using the TCGA-LIHC dataset. We showed that these 10 CRGs were differentially expressed (either upregulated or downregulated) in HCC compared to normal liver (Figure 3A). Specifically, the expression of CDKN2A, DLD, DLAT, LIAS, GLS, LIPT1, MTF1, PDHA1, and PDHB was upregulated, while the expression of FDX1 was significantly downregulated in HCC. Moreover, we analyzed the prognostic values of CRGs by univariate Cox regression analysis. A total of five genes (CDKN2A, GLS, MTF1, DLAT, LIPT1) with a significant p-value were identified (p < .05). These five genes were associated with increased risk of HCC with HRs >1 and recognized as prognostic CRGs (Figure 3B).

To establish a prognostic gene scoring model, LASSO-Cox regression analysis was performed based on the above five prognostic CRGs, and finally all these genes were selected according to the minimum criteria (Figures 4A,B). CRGs-score was calculated by CRGs and the corresponding coefficients (Table 1). Next, we divided the HCC patients into low and high CRGs-score subgroups with the median value of 2.309. The Kaplan-Meier (KM) curve revealed that HCC patients with high CRGs-score had a worse overall survival than those with low CRGs-score (median time = 41.8 vs. 81.7 months, p = .00018, HR = 1.94, Figure 4C). With the CRGs-score increasing, the survival status of patients, and genes expression pattern of these five genes are presented in Figure 4D. The results showed that increased CRGs-score was accompanied by decreased survival and increased risk of death. The expression of CDKN2A, DLAT, and GLS was higher, while MTF1 was lower in high CRGs-score compared to low CRGs-score (Figure 4D). Furthermore, we applied a weighted CRGs-score incorporating all related genes to estimate 1-, 3- and 5-year overall survival. The prediction accuracy evaluated by AUCs was .72, .66, and .62 in the 1-year, 3-year and 5-year ROC curves, respectively (Figure 4E).

Considering the importance of the CRGs-score in evaluating the prognosis of HCC patients, we next investigated its applied value in clinical diagnosis. Univariate and multivariate COX analysis revealed that stage and CRGs-score were independent factors affecting HCC patient’s prognosis (Figures 5A,B). We also constructed a nomogram featuring age, gender, stage, grade and CRGs-score to predict the survival rate of HCC patients at 1, 3, and 5 years (Figure 5C). Intriguingly, the nomogram exhibited good predictive value with a C-index of .67. The calibration of the 1-year, and 3-year overall survival rates were relatively well compared with an ideal model in the entire cohort (Figure 5D). In the KM curves, the disease-free interval (DFI), disease-specific survival (DSS) and progression-free interval (PFI) in the high CRGs-score group were all significantly lower than those in the low CRGs-score group, with the HRs of 1.68, 1.97, and 1.77 respectively (Figures 6A–C). Moreover, patients with advanced stage or higher grade had a higher CRGs-score significantly (Figures 6D,E). These results indicate that the CRGs-score is an excellent predictive prognostic trait.

Next, we further assessed whether the CRGs signature affected tumor immunity. Immune infiltration analysis was performed using CIBERSOFT algorithm. We found that the proportion of CD8⁺ T-cell, gamma delta T-cell, nature Killer cells, M1 macrophage and resting Mast cells were relatively higher in low CRGs-score group, while the proportion of follicular helper T-cell, regulatory T-cell, M0 macrophage, and neutrophils cells were lower in low CRGs-score group compared to high CRGs-score group (p < .05, Figure 7A). Additionally, we also analyzed thirteen immune related pathways between high and low CRGs-score subgroups by ssGSEA method. The high CRGs-score group exhibited downregulation of cytolytic activity, inflammation promoting, type Ⅰ and type Ⅱ IFN response and upregulation of MHC class Ⅰ pathway (p < .05, Figure 7B). We next compared the expression of seven immune checkpoint genes between two groups and found that the expression of immune checkpoint molecules including CTLA4, HAVCR2 (also known as TIM-3), LAG3, PD-L1, PD-1, PDCD1LG2 (PD-L2) and TIGIT. Interestingly, we found that CTLA4, HAVCR2, PD-L1, PD-1 and TIGIT were significantly higher in the high CRGs-score group, indicating that HCC in cuproptosis status might be sensitive to immune checkpoint blockade therapies (p < .05, Figure 7C). Moreover, we performed a correlation analysis between CDKN2A, DLAT, GLS, LIPT1, MTF1 and different immune cell types. Notably, the expression level of CDKN2A, DLAT, GLS, LIPT1, and MTF1 were positively associated with the immune infiltration level of B-cell, CD8⁺ T-cell, CD4⁺ T-cell, Macrophage, Neutrophil, and Dendritic Cell (all p < .05, Figures 8A–E).

The CRG-score in the ICGC validation cohort (n = 240) was calculated with the weighted formula derived from the training cohort, and the HCC patients were also divided into low and high CRGs-score subgroups with the cut-off value (2.309) consistent with the training cohort. The OS of patients with high CRGs-score was obviously decreased (p = .015, Figure 9A). The AUC values of OS at 1, 3, and 5 years predicted by the CRGs-score were .65, .68, and .68 respectively (Figure 9B). Consistently, in both cohorts, high CRGs-score correlated with poor survival. The expression of CDKN2A, DLAT, and GLS was increased, while expression of MTF1 was decreased accompanied with raised CRGs-score (Figure 9C).

We further validated the expression of the five prognostic CRGs (CDKN2A, DLAT, GLS, LIPT1, and MTF1, which constructed the CRGs-score model) in HCC cell lines and mouse HCC model. We analyzed the mRNA level of these five genes in the normal liver cell (LO2) and two HCC cells (Huh7, MHCC97H) by RT-qPCR. As shown in Figure 10, the transcription of CDKN2A, DLAT, GLS, LIPT1, and MTF1 were significantly enhanced in Huh7 and MHCC97H cells compared to LO2 cells (Figures 10A–E). We also established spontaneous HCC mouse model by hydrodynamic tail vein delivery of the NRASV12 and myr-AKT in C57BL/6 mice (Figure 11A). Mice injected with NRASV12 and myr-AKT plasmids successfully developed HCC tumors and exhibited higher liver/body ratio than normal group (Figure 11B). Additionally, increased mRNA level of CDKN2A, DLAT, GLS, LIPT1, and MTF1 in tumor were observed compared to normal group (Figures 11C–G). Furthermore, to investigate tumor immunity, we collected the HCC and normal tissues and detected CD8 expression by Immunohistochemistry (IHC) staining. Notably, CD8⁺T-cell were significantly decreased in HCC compared to normal group (Figures 11H–J).

Elesclomol is an anticancer drug that targets mitochondrial metabolism and plays an important role in transportation copper ions that has been shown to induce cuproptosis. Therefore, we tested the sensitivity of Elesclomol between the low and high CRGs-score subtypes in HCC patients. The elesclomol IC50 of low CRGs-score was significantly lower than high CRGs-score subtype, indicating increased sensitivity (p = 4.5e-06, Figure 12A). To clarify the roles of CRGs on sensitivity of chemotherapeutics and molecular targeted drugs, we also revealed the sensitivity for nine commonly used drugs in HCC between low and high CRGs-score subtypes. The results showed that the IC50 of 5-Fluorouracil, Sunitinib, Gemcitabine, and Bleomycin was lower in high CRGs-score subtype. However, the IC50 of Temsirolimus, Pazopanib, Erlotinib, and Axitinib was higher, compared to low CRGs-score subtype in HCC (Figures 12B–J).