Introduction

Hepatocellular carcinoma (HCC) is a common malignant tumor with a high mortality rate globally. Even though early screening and multidisciplinary diagnosis of HCC have greatly improved, the survival rate of postoperative HCC patients remains unsatisfactory1. Recently, immune checkpoint blockade (ICB) therapy, specifically targeting PD-1/PD-L1, has brought revolutionary progress in the advanced HCC treatment2. Still, only a minority of HCC patients benefit from current ICB therapies. The overall response rate (ORR) of ICB monotherapy is approximately 20%, while the ORR after combination treatment increase to almost 30%3. Immunotherapy based on the application of immune checkpoint inhibitors (ICIs), as single agents or combined with chemotherapeutic drugs, anti-angiogenic agents, and kinase inhibitors, offers great promise for HCC treatment4. Therefore, it is of great significance to focus on the expression pattern of immune checkpoint genes (ICGs) and develop more relevant predictive biomarkers to predict the therapeutic response and prognosis of HCC.

Copper (Cu) is an important metallic trace element in the human body and is engaged in various biochemical functions. Disorder of Cu homeostasis is associated with diverse pathologies, like Menkes disease, hepatic Wilson disease5, and even tumorigenesis6. Several studies have found that Cu levels are significantly upregulated in tumor tissue and serum of cancer patients7,8,9,10,11,12,13. Moreover, alteration of intracellular Cu levels may be involved in cancer development14. Recently, a novel mechanism of Cu-dependent cell death, also termed cuproptosis, was reported by Tsvetkov et al.15,16,17. It was proposed that cuproptosis is attributed to intracellular Cu accumulation inducing aggregation of lipoylated tricarboxylic acid (TCA) cycle components and resulting in proteotoxic stress and cell death16.

Cu homeostasis is essential for maintaining normal immune function18,19,20. Recent studies found that elevated Cu levels in tumor cells contribute to immune escape by enhancing PD-L1 expression21. The relationship between cuproptosis and ICGs has not been elucidated. Therefore, studying the prognostic signature of cuproptosis-related immune checkpoint genes (CRICGs) is of great interest. Based on The Cancer Genome Atlas (TCGA) and Gene expression omnibus (GEO) databases, we identified and validated the prognostic CRICG model. In the meantime, we evaluated the potential correlation between this signature and clinicopathological characteristics, tumor immunity, and somatic mutation.

Methods

Data acquisition and identification of cuproptosis-related ICGs

The transcriptome RNA‐seq dataset (HT Seq-TPM) and the corresponding clinical information of 424 HCC samples (including 374 tumor samples and 50 normal samples) were available from TCGA database (https://portal.gdc.cancer.gov). The progression-free survival (PFS) information of TCGA-HCC samples was acquired from the UCSC Xena platform (https://xenabrowser.net/). According to the research of Tsvetkov et al.16, we focused on 13 cuproptosis-related genes (CRGs). These 13 CRGs included 7 cuproptosis resistant genes (DLD, DLAT, FDX1, LIAS, LIPT1, PDHA1, and PDHB), as well as 3 other key enzymes in the TCA cycle (DBT, GCSH, and DLST) and 3 Cu transport-related factors (SLC31A1, ATP7A and ATP7B). The profile of 79 immune checkpoint genes was obtained as previously reported22. R package limma was applied to extract the gene expression profile of 13 CRGs and 79 ICGs from TCGA database. The CRICGs co-expressed with 13 CRGs were determined with the criteria of Pearson correlation coefficient > 0.2 and p-value < 0.001. The predictive CRGs-CRICGs co-expression network was established using Cytoscape software (https://cytoscape.org).

Research strategy

A flowchart of our study was carried out (Fig. S1). In brief, the gene expression profile of 13 CRGs and 79 ICGs was extracted from the TCGA database. Then, the predictive CRICGs were determined according to Pearson correlation analysis. Subsequently, the overall HCC cohort (n = 370) was randomly classified into the training (n = 185) and testing (n = 185) cohorts. LASSO-COX analysis was performed on the training group to construct the prognostic risk model. Subsequently, internal prognostic validation was performed for the testing and the overall cohorts. Furthermore, clinical correlation, pathway prediction, immune landscape, and TMB correlation analyses were applied to study the potential biological functions of this CRICG signature. Finally, external validation and experimental validation were adopted to verify this signature.

Establishment of CRICG signature

Based on the relevant clinical information, the overall HCC cohort (n = 370) was randomly classified into training and testing cohorts with a ratio of 1:1. The Chi-square test were utilized to verify the effectiveness of random assignment. First, univariate Cox regression analysis was performed on the training cohort for filtering prognostic CRICGs. Next, LASSO regression analysis was adopted to screen CRICGs with the minimum ten-fold cross‐validation using R package glmnet23. Further, multivariate Cox regression analysis was employed to screen CRICGs and determine the corresponding risk regression coefficients. Thus, a prognostic risk model based on 5-ICGs was established. The following equation was utilized to compute the risk score:

$${\text{Risk}}\,\,{\text{score }} = \mathop \sum \limits_{i = 1}^{n} coefi \times expi$$

The expi and coefi indicate the expression levels of ICG and its corresponding coefficients, respectively.

Internal validation of CRICGs signature

The training, testing, and overall cohort of HCC patients were separated into high- and low-risk populations based on the median value of the training risk score. Kaplan–Meier (K–M) survival analyses were carried out using R package survival to evaluate overall survival (OS) values. K–M survival analyses were also performed on subgroups stratified by different clinicopathological characteristics. The independent variables of the risk model were evaluated sequentially by Cox regression analysis. Subsequently, the ROC curves were utilized to assess the validity of this prognostic model by R package timeROC24. The concordance index (C-index) was applied to evaluate the predictive capability of this prognostic model.

Nomogram establishment

Based on the risk group and clinicopathological characteristics, the nomogram of 1-, 3- and 5-years OS was established via R package rms. In addition, the corresponding calibration curves were applied to illustrate the consistency level of the prediction results.

Clinical correlation of prognostic CRICGs signature

The Chi‐square test and Wilcox signed-rank test were applied to study the correlation between risk groups and clinicopathologic factors. The heatmap displayed the clinical correlation findings of the Chi-square test. The boxplot depicted the clinical correlation data of the Wilcox signed-rank test.

Functional enrichment analyses

The three-dimensional (3D) distribution patterns of risk sample classification in different models were presented by principal component analysis (PCA). With the criteria of |log fold change|> 1 and false discovery rate (FDR) < 0.05, the differentially expressed genes (DEGs) among the low- and high-risk populations were determined through R package limma. In order to study the relevant biological processes engaged in this prognostic signature, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)25 analyses were performed on the identified DEGs using R package clusterProfiler (p < 0.05). Furtherly, Gene Set Enrichment Analysis (GSEA) was used to determine the biological functions and pathways significantly enriched in high‐ and low‐risk populations (FDR < 0.05).

Assessment of immune cell infiltration

The tumor microenvironment (TME) scores of HCC samples were evaluated by R package ESTIMATE26,27. TME scores, including estimate scores, immune scores, and stromal scores, could infer tumor purity in tumor tissue. The immune infiltration scores of TCGA samples quantified by multiple algorithms (MCP-counter, quanTIseq, EPIC, CIBERSORT-abs, CIBERSORT, xCell, and TIMER) were obtained from the TIME platform (http://timer.cistrome.org)28. Using the single-sample gene set enrichment analysis (ssGSEA) method with R package GSVA, the degree of 16 immune cells infiltration and the activity level of 13 immune-related functions were evaluated in HCC samples29. The differential expression between risk groups and ICG expression was analyzed via R package limma. Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) platform was applied to forecast the biomarker sensitivity to ICB therapies in tumors like melanoma and non-small cell lung cancer30. The TCGA-HCC expression profiles were imported into the TIDE platform to obtain the TIDE scores of HCC samples. The poor ICB response is associated with high TIDE and dysfunction prediction scores in tumors, as well as low exclusion scores.

Drug susceptibility prediction

To further assess the significance of the prognostic signature in predicting clinical HCC therapy response, the half-maximal inhibitory concentration (IC50) values for chemotherapeutic or targeted drugs in high- and low-risk populations were evaluated via the R package pRRophetic. The Wilcoxon signed-rank test was applied to compare the differences in IC50 values between low- and high-risk populations (p < 0.01).

Somatic mutation analysis and TMB correlation analysis

The simple nucleotide variation information of HCC samples was acquired from the TCGA database. Then, the gene mutation type, frequency, and tumor mutation burden (TMB) scores of HCC samples were calculated using Perl software. The mutation genes among high- and low-risk populations were identified and visualized by R package maftools31.

External validation of CRICGs signature

GSE14520 dataset based on GPL3921 platform was available from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE14520 dataset, containing 221 HCC samples with microarray gene expression and survival information, was adopted to validate the prognostic significance of this CRICG signature externally. K–M survival analyses of OS and PFS were conducted with the optimal cut-off value. The ROC curves of 1-year OS were carried out in the GSE14520 dataset.

Cell culture and qRT-PCR assay

The human normal liver cell (MIHA) and HCC cell (SMMC-7721) were gifts from the College of Basic Medical Sciences, China Medical University, Shenyang, China. Human HCC cell line HEPG2 was maintained in our laboratory. These cell lines were cultured in high-glucose Dulbecco's Modified Eagle Medium (DMEM, Shanghai Basal Media Technologies, China) containing 10% fetal bovine serum (FBS, Biological Industrie, Israel) at 37 °C in a humidified incubator with 5% CO2. As instructed by the corresponding manufacturer, RNAiso Plus (Code No. 9108, Takara) was utilized to extract the total RNA from cell lines, and the PrimeScript RT reagent kits (RR047A, Takara) were used for cDNA synthesis. TB green Master reagents (RR820A, Takara) were used for amplificated through Applied Biosystems 7500 Real-Time PCR System. The relative gene expression level was computed through a 2−ΔΔCt method. All primers were synthesized by Shanghai Invitrogen Biotechnology (Shanghai, China). The primer sequences utilized are displayed in Table 1.

Table 1 The Primer sequences list for qRT-PCR.
Full size table

Statistics analysis

The Student's t-test was used to compare the gene expression differences of qRT-PCR. All statistical analyses were carried out using R software (version 4.1.0). Unless indicated otherwise, p-value < 0.05 was deemed to demonstrate statistically significant.

Results

Identification of CRICGs

The expression heatmap of 13 interested CRGs was shown in Fig. 1A. Based on the Pearson correlation analysis, 32 CRICGs co-expressed with these CRGs were obtained. The predictive CRGs-CRICGs co-expression network was constructed as shown in Fig. 1B.

Figure 1
figure 1

Identification of cuproptosis-related ICGs in HCC. (A) The expression heatmap of 13 cuproptosis‐related genes. (B) The co-expression network between cuproptosis-related genes and ICGs.

Full size image

Establishment and internal validation of prognostic CRICG model

The overall cohort of 370 HCC patients was randomly classified into the training and testing cohorts with a ratio of 1:1. Further, the Chi-square test verified the effectiveness of random assignment (Table 2). The training cohort was then applied to develop the prognostic model. In the first step, 7 ICGs were screened from 32 CRICGs through the univariate Cox regression analysis (Fig. 2A,B). Subsequently, 6 CRICGs were identified by LASSO regression analysis (Fig. 2C,D), 5 of which (TNFRSF4, TNFRSF14, TNFSF4, TNFSF9, and CD226) were screened by means of the multivariate Cox regression analysis (Fig. 2E,F). Then, the risk score was computed using the following equation:

$$ \begin{aligned} Risk\,\,score = & \,TNFRSF{4} \times \left( {0.{4827}} \right) \\ & + \,TNFRSF{14} \times \left( {0.{4436}} \right) \\ & + \,TNFSF{4} \times \left( {0.{28}00} \right) \\ & + \,TNFSF{9} \times \left( {0.{2}0{89}} \right) \\ & + \,CD{226} \times \left( { - {1}.{9}0{82}} \right) \\ \end{aligned} $$
Table 2 The detailed clinicopathologic factors information between the training and testing groups.
Full size table
Figure 2
figure 2

Construction of cuproptosis‐related ICGs. (A, B) The expression profiles of 7 ICGs preliminarily screened by univariate Cox regression analysis. (C, D) The cvfit and lambda curves performed by the LASSO regression model. (E, F) The Sankey diagram and the correlation heatmap of prognostic cuproptosis‐related ICGs screened by the multivariate Cox regression analysis.

Full size image

Additionally, K–M survival analyses were applied to internally validate the prognostic significance of this 5-CRICGs model. In the training, testing, and overall HCC cohorts, the median OS periods of the high-risk populations were remarkably shorter than those of the low-risk populations (Fig. 3A–F). Similar results were observed for PFS (p = 0.008, Fig. S2A). Moreover, subgroups stratified by different clinicopathological characteristics suggest that the high-risk populations are strongly related to unfavorable outcomes (Fig. 3G). The hazard ratio (HR) of this risk score in multivariate Cox regression analysis was 1.213 (p < 0.001, Fig. 4A,B), revealing that this signature could independently predict the HCC outcomes. The ROC curves analyses of 1‐year (AUC = 0.715), 3 years (AUC = 0.671), and 5 years (AUC = 0.684) OS show the high accuracy of this prognostic signature (Fig. 4C). The ROC curves and C-index analyses of this signature were considerably superior to other clinicopathological features, indicating the strong predictive ability of this risk model signature (Fig. 4D,E).

Figure 3
figure 3

Prognosis value evaluation of cuproptosis‐related ICGs model. (A–C) The overall survival value of the training, testing, and overall cohorts. (D–F) The distributions of survival status, risk scores, and expression heatmap in the training, testing, and overall cohorts. (G) Overall survival value stratified by clinicopathologic features.

Full size image
Figure 4
figure 4

Internal validation of prognosis risk model and Nomogram demonstration. (A, B) Cox regression analyses of clinicopathologic features and risk scores. (C) The ROC curves for 1-, 3- and 5-years OS. (D, E) The ROC curves and C-index of risk score and clinicopathologic features. (F, G) Nomogram and calibration curves of the prognosis model.

Full size image

Construction of nomogram

Thus, we constructed a nomogram to forecast 1-, 3- and 5-years OS incidences for HCC patients (Fig. 4F). The calibration curves for each year closely correspond to the diagonal line (Fig. 4G), indicating the reliability of this prognostic model for predicting HCC patient survival.

Clinical value of prognostic signature

Then, we applied the Chi-square test and Wilcox signed-rank test to investigate the correlation among this prognostic signature and clinicopathological features. As illustrated in the heatmap and boxplot, this model signature was markedly related to the grade, Tumor-node-metastasis (TNM) stage, and T stage of HCC (p < 0.001, Fig. 5A,B). In addition, this signature was dramatically associated with the M stage of HCC through Wilcox signed-rank test (p = 0.034, Fig. 5B).

Figure 5
figure 5

Clinical correlation analysis between the prognostic signature and clinicopathologic features by (A) the Chi-square test and (B) Wilcox signed-rank test. ***p < 0.001.

Full size image

Functional enrichment analyses

The 3D PCA analysis demonstrated the feasibility and availability of the CRICGs model in risk sample classification (Fig. 6A). Next, the DEGs between the low‐ and high‐risk populations were identified by R package limma. The GO analysis illustrated that these DEGs were enriched in various biological processes and molecular functions associated with cell replication, including chromosome segregation, nuclear division, and microtubule binding (Fig. 6B). The KEGG analysis suggests that these DEGs are enriched in cell cycle, drug metabolism-cytochrome p450, metabolism of xenobiotics by cytochrome p450, and retinol metabolism (Fig. 6C). The GSEA analysis indicated that high-risk populations were involved in several tumor-related pathways, such as cell cycle, DNA replication, p53 signaling pathway, and DNA repair pathways32(base excision repair, nucleotide excision repair, mismatch repair, and homologous recombination) (Fig. 6D). While the low-risk populations were involved in many kinds of immune or metabolic pathways, including complement and coagulation cascades, retinol metabolism, fatty acid metabolism, and PPAR signaling pathway (Fig. 6E).

Figure 6
figure 6

Functional analyses of cuproptosis‐related ICGs prognostic signature. (A) The PCA scatterplot of sample 3D distribution. (B) GO and (C) KEGG enrichment analysis of risk differential genes. (D, E) GSEA analysis of prognostic signature.

Full size image

Immune-related analysis of prognostic signature

We next employed the ESTIMATE method to study the correlation with the immune microenvironment. As illustrated in Fig. 7A, the stromal score is relatively lower in high-risk populations (p < 0.01), suggesting low tumor purity in high-risk populations. In addition, immune infiltration scores quantified by multiple algorithms reveal that this 5-CRICG signature was strongly associated with immune cell infiltration (Fig. 7B). Next, we adopted the ssGSEA method to investigate the immune-related functions of this prognostic signature. Compared with the low-risk populations, the high-risk populations had a relatively increased degree of aDCs and macrophage cell infiltration and a moderately decreased degree of Neutrophils, T helper cells, and natural killer (NK) cells infiltration (Fig. 7C). In addition, the high-risk populations had rather attenuated activity in several immune functions in terms of type II IFN response and cytolytic activity (Fig. 7D). Besides, the other immune checkpoint gene expression levels among these two populations were compared. As shown in Fig. 7E, 28 immune checkpoint genes (including these 5 ICGs) were differentially expressed. Various immune checkpoint genes as effective immunotherapy targets are highly expressed in high-risk populations, including CTLA4, TIGIT, and CD276 (B7-H3). Notably, the TIDE scores (p = 0.0031) and dysfunction scores (p = 0.0028) were relatively lower in the high-risk populations, while the exclusion scores (p < 0.001) were significantly higher, indicating that the high-risk populations might correlate with better ICB responses (Fig. 7F).

Figure 7
figure 7

Immune-related analyses of cuproptosis‐related ICGs prognostic signature. (A) TME score prediction. (B) The immune cell infiltration by multiple methods. The correlation between the predictive signature and (C) 16 immune cells and (D) 13 immune-related functions. (E) The relationship between prognostic signature and immune checkpoints. (F)TIDE score prediction. *p < 0.05; **p < 0.01; ***p < 0.001.

Full size image

Drug susceptibility analysis

After that, we screened 10 chemotherapeutic or targeted drugs through drug sensitivity comparison (Fig. S3). The drug sensitivity results showed that the high-risk population was more resistant to Elesclomol (cuproptosis inducer), verifying the consistency of this model signature. In addition, the high-risk populations were more resistant to Camptothecin (DNA topoisomerase I inhibitor), Axitinib (multi-targeted tyrosine kinase inhibitor), Erlotinib (EGFR tyrosine kinase inhibitor), and Gefitinib (EGFR tyrosine kinase inhibitor); and more sensitive to 5-Fluorouracil (DNA/RNA Synthesis inhibitor), Doxorubicin (DNA Topoisomerase inhibitor), Mitomycin C (DNA cross-linking agent), Paclitaxel (microtubule-associated inhibitor), and Sunitinib (multi-targeted receptor tyrosine kinase inhibitor).

Somatic mutation and TMB analysis

Subsequently, we evaluated the somatic mutation status among the high- and low-risk populations. The frequency top 10 mutated genes were TP53, CTNNB1, TTN, MUC16, PCLO, ALB, RYR2, APOB, CSMD3, and LRP1B (Fig. 8A,B). Among them, the frequencies of gene mutations in TP53, MUC16, PCLO, CSMD3, and LRP1B were increased in the high-risk populations. Besides, a significant difference was shown in the survival analysis of the risk groups combined with TMB (p < 0.001, Fig. 8C). The high-risk populations combined with high TMB scores had the worst OS.

Figure 8
figure 8

Somatic mutation analysis and external validation of cuproptosis‐related ICGs model. (A, B) The waterfall plot showing the difference in somatic mutation features between the low- and high-risk populations. (C) Kaplan–Meier survival curves of the prognostic signature combined with TMB. (D) OS in the GSE14520 dataset. (E) The ROC curves for 1-year OS in the GSE14520 dataset. (F) Expression levels of TNFRSF4, TNFRSF14, TNFSF4, TNFSF9, and CD226 in HCC cell lines. *p < 0.05; ***p < 0.001.

Full size image

External validation of prognostic signature

We externally validated this prognostic signature with the GSE14520 dataset. K–M survival curves for OS indicate that high-risk populations strongly correlate with unfavorable outcomes (p = 0.002, Fig. 8D). The AUC of 1-year OS was 0.601 (Fig. 8E). The PFS analysis yielded identical results (p = 0.005, Fig. S2B). Thus, we further evaluated the expression levels of TNFRSF4, TNFRSF14, TNFSF4, TNFSF9, and CD226 in the HCC cell lines. We observed significant downregulation of CD226 and TNFRSF4 in HCC cell lines, as well as prominent upregulation of TNFRSF14, TNFSF4, and TNFSF9 (Fig. 8F).

Discussion

Cu is an important metallic trace element in the human body and is engaged in various biochemical functions such as mitochondrial respiration, iron absorption, autophagy, and antioxidant defense33. The liver is the main organ responsible for removing Cu from the human body34,35. In hepatocytes, the Cu excretion is primarily mediated by the copper exporter ATP7A/B36,37,38,39. The disorders of Cu metabolism are closely related to liver diseases, like Wilson’s disease and even hepatocarcinoma. Elevated levels of Cu are closely correlated with various malignant cancers. The Cu imbalance plays a prominent role in cancer development40. In particular, the Cu level may increase the incidence of hepatocarcinoma in Wilson’s disease patients41.

A growing number of studies have explored predictive models for ICGs as potential predictors in multiple cancers. Li et al.42 investigated a signature of 7 ICGs to predict survival outcomes and therapeutic responses for endometrial carcinoma patients. Liu et al.43 identified a signature of 3 ICGs to predict prognosis and immune status in triple-negative breast cancer. The latest research study by Tsvetkov et al. revealed a novel form of Cu-dependent programmed death named cuproptosis. Elevated Cu levels in tumor cells contribute to immune escape by enhancing PD-L1 expression. Thus, studying the connection of cuproptosis with other ICGs is of great interest. Here, we developed a novel prognostic CRICG model, providing a new theoretical foundation for CRICG targets as the future hepatocarcinoma therapy.

Firstly, we focused on 13 CRICGs based on Tsvetkov et al.16 research. These genes included 7 cuproptosis resistant genes identified by whole-genome CRISPR-Cas9 positive selection screening. In addition, 3 Cu exporters or importers (ATP7A, ATP7B, and SLC31A1) could regulate cuproptosis sensitivity by affecting intracellular Cu levels16. The other 3 critical enzymes in the TCA cycle (DBT, GCSH, and DLST) were also included. Alterations in the expression and catalytic activity of essential enzymes of the TCA cycle are associated with the development of HCC44,45.

According to the LASSO-COX regression analyses, we developed a prognostic signature based on 5-CRICGs (TNFRSF4, TNFRSF14, TNFSF4, TNFSF9, and CD226). This prognostic signature was strongly related to the TNM stage, T stage, M stage, and grade of HCC. GO and KEGG analyses showed that the DEGs between these two populations enriched in several biological functions like cell division, DNA replication, and drug response. GSEA results reveal significant differences between these two populations enriched in various tumor-associated signaling pathways, including cell cycle, DNA replication, p53 signaling pathway, and DNA repair pathways. The above results suggest that this signature is remarkably correlated with the tumorigenesis and progression of HCC.

Again, we investigated the potential connection between the prognosis model and the immune landscape. According to the multiple immune infiltration quantification algorithms, this 5-CRICGs signature was strongly related to the immune cell infiltration. Analyzed by the ssGSEA method, this prognosis signature was attenuated in several immune cell infiltration (including neutrophils, NK cells, and T helper cells) and multiple immune functions in terms of cytolytic activity and type II IFN response. It has been evidenced that type II IFN (also known as IFNγ) could suppress tumor growth by affecting immune cells and tumor cells46. The loss of IFNγ is also pinpointed as related to ICB response47. These results suggest that the high-risk populations might be in an immunosuppressive state by modulating the IFN response. High-risk populations had prominently lower TIDE scores, indicating that the high-risk populations might correlate with preferable ICB responses. In addition, the drug sensitivity prediction of this CRICG prognostic model in response to common chemotherapeutic or targeted drugs can provide the possibility of personalized treatment for HCC patients.

Finally, the GSE14520 dataset externally validated the prognostic value of this CRICG signature. Furthermore, the expression levels of these 5 CRICGs in the HCC cell lines were detected. These 5 CRICGs were characterized as co-stimulatory molecules in the immune response48. CD226 and TNFRSF4 were expressed relatively lower in HCC cell lines, while TNFRSF14, TNFSF4, and TNFSF9 were the opposite trend. Except for TNFSF4, the HCC cell lines expression differences of CD226, TNFSF9, TNFRSF4, and TNFRSF14 correspond to Fig. 2B above. CD226, famous for competing with TIGIT, is an important co-stimulator that regulates T cell and NK cell function49. It was demonstrated that CD226 could induce IFNγ generation50 and promote HCC cell lysis51,52. Recent studies have found that CD226 expression is essential for PD-(L)1 or TIGIT immunotherapy response53,54. TNFSF9, also known as CD137L or 4-1BBL, was reported to be upregulated in primary biliary cirrhosis55. High TNFRSF14 expression was considered associated with colorectal liver metastasis56. Recently, TNFRSF4 (OX40) has emerged as a potential target for cancer immunotherapy57. TNFRSF4 binds its ligand TNFSF4 (OX40L) and transmits co-stimulatory signals58. The up-regulation of TNFRSF459,60,61 and down-regulation of TNFSF461 in our HCC cell lines experiment were consistent with the previous studies.

Here, we constructed and validated a prognostic risk model based on 5 CRICGs. This risk model is feasible and effective. However, our research has some limitations. First, the detailed information on treatment options and therapeutic responses of HCC patients in the TCGA database is unavailable. In addition, we only applied the cell lines for expression level validation due to the acquisitiveness constraints of tissue samples and commercial microarrays. Finally, the biological functions and potential mechanisms of these ICGs in hepatocellular carcinoma have not been explored in depth. Furtherly, we will continue to collect clinical tissue samples and design more in vitro experiments to investigate the function of these prognostic ICGs.

Conclusions

Our prognostic signature based on 5 CRICGs could validly forecast the outcomes and immune response of HCC patients, and is promising as biomarkers to be developed for anticancer therapy.