On 15 July 2022, We downloaded the RNA sequencing (RNA-seq) data and the corresponding clinical information from 371 HCC patients using the Cancer Genome Atlas (TCGA) database, and the clinical information of HCC patients are shown in Supplementary Table S1. In addition, somatic cell datasets and Copy number variation (CNV) data for HCC have been downloaded from the TCGA database and the University of California Santa Cruz (UCSC) Xena websites. Data analysis was performed using R (version 4.0.3) and R Bioconductor packages. The expression data were normalized to transcripts per kilobase million (TPM) values before further analysis.

Based on previous studies (Polishchuk et al., 2019; Aubert et al., 2020; Dong et al., 2021; Ren et al., 2021; Bian et al., 2022; Chen, 2022; Kahlson and Dixon, 2022; Tsvetkov et al., 2022), a total of 19 CRGs (NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, DBT, GCSH, DLST) were obtained. Differences in CRGs expression between HCC and adjacent normal tissues were identified using the “ggplot2” package. We then constructed a protein-protein interaction (PPI) network (the minimum required interaction score of 0.9) for 19 CRGs using the Search Tool for Retrieving Interacting Genes (STRING).

Using the “maftools” package generated Mutation frequencies and oncoplot waterfall plots for 19 CRGs in HCC patients. Using the “RCircos” package in R, the locations of CNV alterations in 19 CRGs on 23 chromosomes were mapped.

By http://vip.sangerbox.com/home.html, 19 CRGs were mapped in the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis for enrichment analysis in HCC.

A prognostic overall survival (OS) analysis was performed to identify potential prognostic characteristics. Log-rank tests were performed to calculate p-values, hazard ratios (HRs), and 95% confidence intervals. Based on these prognostic CRGs, a prognostic model was constructed using LASSO Cox regression analysis. The resulting model risk score was formulated as follows: risk score = exp-gene1 * coef-gene1 + exp-gene2 * coef-gene2+ …+ exp-genei * coef-genei. According to the median risk score, TCGA-HCC patients were divided into low-risk and high-risk subgroups. The Kaplan-Meier analysis compared the two subgroups’ overall survival (OS) time. Predictive accuracy and risk scores were assessed for each gene by performing a temporal receiver operating characteristic (ROC) analysis. Predictive nomograms were developed to predict overall survival at 1, 3, and 5 years.

Then, we used the Tumor Immunity Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/), a web portal that comprehensively analyzes tumor-infiltrating immune cells to analyze the association between prognostic CRGs and immune infiltration. TIMER’s “Genes” module allows visualizing the correlation of gene expression with immune infiltrate levels in HCC. In the analysis of TMB and MSI, Spearman correlation analysis was used to calculate the correlation between gene expression and TMB and MSI scores, and p < 0.05 was considered statistically significant.

We analyzed the association of five prognostic genes with tumor grade and pTNM stage using the “ggplot2” package. HCC patients were divided into four phases for gene expression comparison and into early and late phases for comparison with adjacent normal tissues.

We selected CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT and SIGLEC15 as immune-checkpoint-relevant transcripts extracted the expression data of these eight genes to assess the expression of the immune checkpoints and co-expression-of-prognostic CRGs with these immune checkpoints.

We collected 6 pairs of HCC tissues and adjacent normal tissues from the Third Xiangya Hospital of Central South University. These tissues were used to detect CDKN2A expression levels by qRT-PCR and immunohistochemistry. The Ethics Committee of the Third Xiangya Hospital of Central South University has approved the study. The approval number is I-22287.

To lyse cells, 1 ml Trizol reagent (TaKaRa, Japan) was added to the sample tissue according to the actual manufacturer’s instructions and incubated on a shaker for 15 min at room temperature. Total RNA was subsequently extracted from target samples. One microgram of RNA was reverse transcribed into cDNA using Revert Aid First Strand cDNA Synthesis Kit (Thermo, United States). Quantitative RT-PCR was then performed with Pro Taq HS Premix Probe qPCR Kit (Accurate, Hunan, China). The GAPDH gene was used as an endogenous control gene for normalizing the expression of target genes. Primers used in this study included CDKN2A(forward 5′-CCG​TGG​ACC​TGG​CTG​AGG​AG-3′, reverse 5′-CGG​GGA​TGT​CTG​AGG​GAC​CTT​C-3′), GAPDH (5′-CAG​GAG​GCA​TTG​CTG​AT-3′, 5′-GAA​GGC​TGG​GGC​TCA​TTT-3′).

The tissues were sectioned and embedded in paraffin. Sections were incubated overnight at 4°C with anti-CDKN2A antibody (1:200 dilution; proteintech, China). Slides were washed with phosphate-buffered saline (PBS) and incubated with a goat anti-rabbit IgG secondary antibody conjugated with fluorescein isothiocyanate (ZSDB-BIO, China) for 30 min with washed slides. After washing with PBS, they were incubated with an antifade reagent (Invitrogen, United States). Staining was visualized to determine protein expression levels using an Olympus CX41 fluorescence microscope (Olympus, Japan). The analysis results were performed using ImageJ software.

We first explored the different expressions of 19 CRGs in HCC and the adjacent normal tissues by the UCSC Xena. In HCC, 15 CRGs were up-regulated or down-regulated (Figure 1). More specifically, compared to adjacent normal tissues, ATP7A, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A and DLST was increased, whereas the expression of NLRP3, SLC31A1, and DBT was decreased in HCC. A PPI analysis was constructed to detect the interaction of these CRGs, and the results showed that DLAT, DLST, PDHA1, DLD, and LIPT1 were central genes (Supplementary Figure S1A). Supplementary Figure S1B shows the correlation network containing all CRGs.

Then we summarize the incidence of copy number variants and somatic mutations in 19 CRGs in HCC. Figures 2A, B showed 38 of 368 (10.33%) HCC samples showed gene mutations. Missense mutations were the most ordinary variant classification (Figure 2A). SNPs were the most ordinary variant type, and C > A was ranked as the top SNV category. Among 19 CRGs, CDKN2A had the highest mutation frequency, followed by NFE2L2 and NLRP3 (Figure 2B). Figure 2C demonstrates the location of CNV alterations on chromosomes for these 19 CRGs. We also investigated the frequency of CNV alterations and showed that these 19 CRGs showed prevalent CNV alterations. Among 19 CRGs, NLRP3, DLD, and PBHB had copy number amplification, while CNV deletion frequencies were prevalent for ATP7B, PDHA1, MTF1, CDKN2A, and DLST (Figure 2D).

Using GO and KEGG databases, pathways were analyzed to reveal the function of CRGs. We found that GO was divided into three categories: biological pathway (BP), cytological component (CC), and molecular function (MF). In the BP, these 19 CRGs are mainly involved tricarboxylic acid cycle, acetyl COA metabolic process, aerobic respiration and immune response (Figure 3A). In the CC, these genes mainly participated in mitochondrion, oxidoreductase matrix, catalytic complex and pyruvate dehydrogenase complex (Figure 3B). In the MF, these genes mainly involved transferase activity transferring acyl groups, copper ion transmembrane transporter activity and DNA binding transcription factor activity (Figure 3C). In addition, KEGG pathway analysis showed that 19 CRGs have mainly participated in metabolic pathways, Glycolysis/Gluconeogenesis, Citrate cycle (TCA cycle), Platinum drug resistance and Human T-cell leukemia virus-1 infection (Figure 3D).

A univariate Cox regression analysis was performed to find CRGs with a prognostic value, and five CRGs were screened. As a result, the forest plot showed p-values for the 19 CRGs (Figure 4A). The results showed that the low expression group of HCC patients with CDKN2A (Figure 4B, p = 0.002), DLAT (Figure 4C, p = 0.002), DLST (Figure 4D, p = 0.014) and GLS (Figure 4E, p = 0.022) and PDHA1 (Figure 4F, p = 0.023) had better survival than the high expression group. A five-gene model was constructed according to the optimum λ value obtained by LASSO Cox regression analysis (Figures 5A, B). Risk score was calculated as Risk score = (0.1645) * CDKN2A + (0.3033) * DLAT + (-0.096) * DLST + (0.1196) * GLS + (0.0289)* PDHA1. Based on the median score, the distribution of risk scores and survival status, patients in the training set were randomly divided into high and low-risk groups, and the distribution of risk scores, survival status, and expression of these five genes are shown in Figure 5C, with patients having an increased risk of death and a decreased survival time as the risk score increased. Kaplan-Meier curves showed that high-risk HCC patients had a lower overall survival probability than low-risk patients (Figure 5D, p = 0.00155), and 1, 3, and 5 years ROC curves were showed AUCs of 0.735, 0.643, and 0.631, respectively (Figure 5E).

Later, we constructed a predictive nomogram to predict the survival probability. The five prognostic CRGs and the clinicopathological features were considered in the nomogram. CDKN2A expression and pT stage were found they be independent factors that could affect the prognosis of HCC patients in univariate and multivariate analysis (Figures 6A, B). Predicted nomograms indicated that 1-year, 3-year and 5-year OS curves could be predicted relatively well compared to the ideal model in the entire cohort (Figures 6C, D).

In developing the tumor’s immune microenvironment, cuproptosis plays a crucial role. Our study also elucidated the association between the expression of prognostic CRGs (DLAT, PDHA1, GLS, CDKN2A, DLST) and HCC immune infiltration using the TIMER database. The results showed that CDKN2A expression was positively associated with the expression of B cells (p = 1.08e-06), CD8⁺ T cells (p = 9.56e-07), CD4⁺ T cells (p = 5.26e-04), neutrophils (p = 1.33e-06), macrophages (p = 6.01e-06), and medullary dendritic cells (p = 1.40e-08) (Figure 7A). DLAT expression was positively associated with the expression of B cells (p = 6.19–04), CD8⁺ T cells (p = 1.64e-02), CD4⁺ T cells (p = 1.66e—05), neutrophils (p = 2.00e—13), macrophages (p = 1.03e-08), and medullary dendritic cells (p = 2.13e-07) (Figure 7B). DLST expression was positively associated with the expression of B cells (p = 1.91e-01), CD8⁺ T cells (p = 6.64e-01), CD4⁺ T cells (p = 1.90e-03), neutrophils (p = 5.24e-05), macrophages (p = 8.62e-02), and medullary dendritic cells (p = 1.91e-02) (Figure 7C). Moreover, GLS expression was also positively associated with the expression of B cells (p = 8.66e-05), CD8⁺ T cells (p = 2.90e-04), CD4⁺ T cells (p = 6.60e-22), neutrophils (p = 9.07e-20), macrophages (p = 6.77e-23), and medullary dendritic cells (p = 1.21e-07) (Figure 7D). PDHA1 expression was positively associated with the expression of B cells (p = 9.29e-04), CD8⁺ T cells (p = 1.28e-03), CD4⁺ T cells (p = 5.01e-03), neutrophils (p = 8.72e-13), macrophages (p = 2.81e-10), and medullary dendritic cells (p = 8.07e-04) (Figure 7E). This study suggests a significant association between CRGs and tumor-immune infiltration.

We analyzed the association of prognostic genes with tumor grade and pTNM stage. The expression of CDKN2A (Figures 8A, p = 2.9e-05), DLTA (Figure 8B, p = 0.027), GLS (Figure 8D, p = 0.0056) was significantly different at the pTNM stage in HCC. In addition, CDKN2A (Figure 8F, p = 2e-24), DLAT (Figure 8G, p = 8.8e-09), DLST(Figure 8H, p = 0.00047), GLS(Figure 8I, p = 2.5e-13) and PDHA1 (Figure 8J, p = 8.3e-14) were significantly different at the in HCC patients with high tumor grades, low tumor grades and adjacent normal tissues.

The above results showed that CRGs were significantly associated with tumor immune infiltration. Subsequently, we analyzed the relationship between CRGs and TMB and MSI, which could further elucidate whether CRGs could serve as a biomarker for drug screening and lay the foundation for subsequent studies. The results showed a positive correlation between MSI and PDHA1 (Figure 9H, p = 0.044) but no significant correlation between MSI and CDKN2A, DLAT, DLST and GLS (Figures 9A, C, E, F). Meanwhile, TMB was positively correlated with CDKN2A (Figure 9B, p = 1.64e-4) and negatively correlated with GLS (Figure 9G, p = 7.74e-05) but not significantly associated with DLAT, DLST, PDHA1 (Figures 9D, J, I).

Drug therapy is an important means of cancer treatment. In this study, drug sensitivity analysis and immune checkpoint analysis provide a certain basis for the mechanism of drug therapy. In drug sensitivity analysis, the expression of CDKN2A, DLAT, PDHA1, and GLS was negatively associated with some or most drugs in the Cancer Therapy Response Portal database, and the expression of DLST positively correlated (Figure 10A). HCC samples were divided into high-risk and low-risk groups and compared with adjacent normal tissues for immune checkpoint analysis. The expression of CD274, CTLA4, LAG3, PDCD1, PDCD1LG2 and SIGLEC15 was significantly different between adjacent normal tissues and tumor tissues (Figures 10B, C). Then, we performed a correlation analysis between five CRGs with prognostic value and immune checkpoint genes. We found that PDHA1 and CD274, PDCD1LG2, GLS were significantly correlated with CD274, CTLA4, PDCD1, PDCD1LG2, DLST and CD274, PDCD1LG2, SIGLEC15, DLAT and CD274, PDCD1LG2, CDKN2A was significantly correlated with CD274, CTLA4, LAG3 and PDCD1, while DLTA was negatively correlated with LAG3 (Figure 10D). Based on the positive correlation between CD274 and all five prognosis-related genes, we performed another analysis of CD274 in HCC and showed a significant difference between high-risk and low-risk HCC and adjacent normal tissues (Figure 10E, p = 0.0033).

qRT-PCR was performed to validate CDKN2A expression levels in HCC and adjacent normal tissues. The results showed the expression level of CDKN2A in HCC was significantly higher than that in adjacent normal tissues (Figure 11A), which was consistent with the results of gene difference analysis. The protein expression level of CDKN2A was subsequently determined by immunohistochemical staining, and the results were consistent with the expression level of mRNA (Figure 11B).