Prognostic Implication of a Novel Metabolism-Related Gene Signature in Hepatocellular Carcinoma

Background: Hepatocellular carcinoma (HCC) is one of the main causes of cancer-associated deaths globally, accounts for 90% of primary liver cancers. However, further studies are needed to confirm the metabolism-related gene signature related to the prognosis of patients with HCC.

Methods: Using the “limma” R package and univariate Cox analysis, combined with LASSO regression analysis, a metabolism-related gene signature was established. The relationship between the gene signature and overall survival (OS) of HCC patients was analyzed. RT-qPCR was used to evaluate the expression of metabolism-related genes in clinical samples. GSEA and ssGSEA algorithms were used to evaluate differences in metabolism and immune status, respectively. Simultaneously, data downloaded from ICGC were used as an external verification set.

Results: From a total of 1,382 metabolism-related genes, a novel six-gene signature (G6PD, AKR1B15, HMMR, CSPG5, ELOVL3, FABP6) was constructed based on data from TCGA. Patients were divided into two risk groups based on risk scores calculated for these six genes. Survival analysis showed a significant correlation between high-risk patients and poor prognosis. ROC analysis demonstrated that the gene signature had good predictive capability, and the mRNA expression levels of the six genes were upregulated in HCC tissues than those in adjacent normal liver tissues. Independent prognosis analysis confirmed that the risk score and tumor grade were independent risk factors for HCC. Furthermore, a nomogram of the risk score combined with tumor stage was constructed. The calibration graph results demonstrated that the OS probability predicted by the nomogram had almost no deviation from the actual OS probability, especially for 3-year OS. Both the C-index and DCA curve indicated that the nomogram provides higher reliability than the tumor stage and risk scores. Moreover, the metabolic and immune infiltration statuses of the two risk groups were significantly different. In the high-risk group, the expression levels of immune checkpoints, TGF-β, and C-ECM genes, whose functions are related to immune escape and immunotherapy failure, were also upregulated.

Conclusions: In summary, we developed a novel metabolism-related gene signature to provide more powerful prognostic evaluation information with potential ability to predict the immunotherapy efficiency and guide early treatment for HCC.

Introduction

As the sixth most common cancer worldwide, liver cancer has become the fourth most prevalent cause of tumor-associated deaths globally (1). Hepatocellular carcinoma (HCC) is considered the most common pathological type of primary malignant liver cancer and has become a significant global health concern (2). At the time of diagnosis, the vast majority of HCC patients are in the middle-to-late stages of the disease. Although surgery and multi-mode treatment have been improved, the overall 5-year survival rate of HCC patients remains low, approximately 10–20%, due to its invasiveness, metastatic potential, and recurrence rate (3). Therefore, identifying effective biomarkers for HCC prognosis is important for the evaluation and treatment of HCC. Common prognostic indicators of HCC include clinicopathological features, such as alpha-fetoprotein and microvessel infiltration (4, 5). However, the specificity and sensitivity of existing prognostic markers do not provide meaningful prognosis patterns.

Recent studies have indicated that metabolic changes associated with cancer represent a novel concern in managing this disease. Many studies have confirmed that changes in cell metabolism play a critical role in the initiation and progression of cancer (6–8). The liver is an important metabolic organ, playing a critical role in many important metabolic pathways including glycolysis, the tricarboxylic acid cycle, and nucleotide biosynthesis (9). Therefore, elucidating the relationship between metabolism-related genes and HCC is essential for understanding its pathogenesis (10). However, the value of a metabolism-related gene signature in HCC prognosis evaluation remains unclear.

In this study, we first downloaded transcriptome profiling data and the corresponding clinical information of liver hepatocellular carcinoma (LIHC) from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) databases, respectively. We then used the least absolute contraction and selection operator (LASSO) regression analysis to construct a six-gene signature related to metabolism in the TCGA set, and verified it in the ICGC set. Using the six-gene signature, we accurately predicted the overall survival (OS) of patients with HCC. Finally, we identified six differentially expressed genes (DEGs) with prognostic value to construct a signature, which represents a promising predictive indicator for HCC patients.

Materials and Methods

Data Collection

We downloaded transcriptome profiling data (fragments per kilobase million, FPKM), up to and including October 20, 2020, of 374 HCC patients and 50 normal tissues, along with the corresponding clinical information of 377 patients with HCC from the TCGA database (https://portal.gdc.cancer.gov/repository), as a training set. The clinical information included survival time, survival status, age, sex, histological grade, tumor stage, and TNM stage. Clinical characteristics are shown in Table 1. Subsequently, these data were matched according to the sample names, and the date of 343 HCC patients who survived more than 30 days was used for follow-up analyses to rule out non-neoplastic causes of death, including heart failure, infection, and bleeding. We then evaluated the other 239 HCC samples by using complete expression profiling data downloaded from the ICGC database (https://dcc.icgc.org/) as a validation set. These samples included liver cancer tissues from Japan with a background of hepatitis B virus (HBV) and hepatitis C virus (HCV) infection. The clinical information of 259 HCC patients, including survival time, survival status, age, sex, and tumor stage of the corresponding patients, was also obtained and extracted (Table 2). After matching these data according to sample names, 231 patients with HCC were evaluated. Homo sapiens GTF files in the Ensembl database were used to annotate the probes.

TABLE 1

Table 1 Clinical characteristics of the HCC patients in the TCGA database.

TABLE 2

Table 2 Clinical characteristics of the HCC patients in the ICGC database.

Establishment of a Prognostic Gene Signature That Related to Metabolism

Metabolism-related genes were collected from the Molecular Signatures Database (MSigDB) (REACTOME_METABOLISM_OF_AMINO_ACIDS_AND_DERIVAT IVES M727, REACTOME_METABOLISM_OF_CARBOHYDRATES M16864, and REACTOME_METABOLISM_OF_LIPIDS M27451). The “Limma” R package was used to screen metabolism-related DEGs between tumor and paired adjacent non-tumor tissues (FDR <0.05, |log₂FC| >3). We then applied univariate Cox analysis to identify the prognostic metabolism-related DEGs. Based on these overlapping prognostic DEGs, an interactive network was generated using the STRING dataset. To avoid overfitting, the “glmnet” R software package was used for LASSO analysis to construct a prognostic signature. The following formula was applied to the definition of the risk score for each patient (11):

$$\text{Risk score}=\beta \mathrm{gene}(1)\times \text{E gene}(1)+\beta \text{gene}(2)\times \text{E gene}(2)+\dots +\beta \text{gene}(\text{n})\times \text{E gene}(\text{n})$$

E denotes the normalized expression level of the gene, and β denotes the corresponding regression coefficient. Patients obtained from the TCGA and ICGC databases were divided into two risk groups respectively, according to the median risk score of the TCGA dataset.

Verification of the Prognostic Evaluation Efficiency of the Gene Signature

To identify the correlations between gene signature and the clinical outcomes of HCC patients, we conducted Kaplan–Meier survival analysis to assess the relationship between risk grouping and OS. The “survivalROC” R package was utilized to execute receiver operating characteristic (ROC) analysis to estimate the specificity and sensitivity of prognostic evaluation of gene signature prognosis assessment. Principal components analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were used to investigate the distribution patterns of the two groups. Moreover, we conducted an independent prognostic analysis to evaluate the independent risk predictors for the prognosis of HCC patients.

Verification of the mRNA and Protein Expression of Six Metabolism-Related Genes Between HCC and Adjacent Normal Liver Tissues

Three paired HCC samples and adjacent normal liver tissue samples were collected. Total RNA was extracted from the tissues by using RNAiso Plus (Takara, Otsu, Japan) according to the manufacturer’s instructions. Next, we reverse-transcribed RNA into cDNA by HiScript^® III RT SuperMix for qPCR Kit (Vazyme). The ChamQ SYBR qPCR Master Mix Kit (Vazyme) was used to conduct real-time reverse transcription polymerase chain reaction (RT-PCR) using Bio‐Rad CFX96™ (Bio‐Rad, Hercules, California, USA). Gene expression was standardized as the beta-actin gene ACTB. The sequences of the primers for the prognostic genes and ACTB are presented in Table S1. Each sample was repeated three times. The 2^−ΔΔCt method was used to compare the relative expression levels of metabolism-related genes in HCC and adjacent normal liver tissues. Representative immunohistochemical staining images of the six metabolism-related genes in HCC and adjacent normal liver tissues were downloaded from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/).

Construction of Nomogram

Nomogram analysis was carried out with the “rms” R package, by fitting the predictors significantly related with OS of HCC in multivariate analysis into research to predict the OS of patients with HCC. A calibration plot was constructed to evaluate the discrimination between the 1-, 2-, and 3-year OS predicted by the nomogram, and the actual values. Furthermore, the reliability of the nomogram was estimated using the concordance index (C-index) and decision curve analysis (DCA).

Functional Enrichment Analysis

We further conducted a gene set enrichment analysis (GSEA) for functional annotation. Moreover, single-sample gene set enrichment analysis (ssGSEA) in the “GSVA” R package was performed to quantify the difference in tumor-infiltrating immune cell scores and immune-related pathway activity between the two risk groups.

Statistical Analysis

R software and GraphPad Prism 7 were used for all statistical analyses. Differences between the two groups were assessed using the Student’s t-test. Statistical significance was set at P <0.05.

Results

Acquisition of Prognostic Metabolism-Related DEGs

A flowchart of this process is displayed in Scheme 1. A total of 1,382 metabolism-related genes were collected from the MSigDB (Table S2). Among them, 40 metabolism-related genes were differentially expressed between HCC and adjacent paired normal tissues (FDR <0.05, |log₂FC| >3, Table S3). Meanwhile, 162 genes related to HCC prognosis were screened out using univariate Cox analysis (P <0.001, Table S4). Subsequently, six prognostic DEGs (G6PD, AKR1B15, HMMR, CSPG5, ELOVL3 and FABP6) were identified for further analysis (Figure 1A). Figures 1B, C present the heatmap and forest plots of the six prognostic DEGs. All six genes were defined as deleterious effectors with all-hazard ratio (HR) values >1. The correlation network between the six metabolism-related genes is shown in Figure 1D.

SCHEME 1

Scheme 1 Flow chart of data collection and analysis.

FIGURE 1

Figure 1 Identification of the candidate metabolism-related genes in the TCGA set. (A) Venn diagram to identify prognostic DEGs between tumor and adjacent normal tissue. (B) The heatmap of the six genes between tumor and adjacent normal tissue. (C) The forest plot between gene expression and OS. (D) The correlation network of six metabolism-related genes.

Establishment of a Prognostic Metabolism-Related Gene Signature

LASSO analysis was performed to identify a six-gene signature according to the optimal value of λ. The risk score of patients was determined as follows: Risk score = 0.209817253236067 × E (G6PD) + 0.123704944063381 × E (AKR1B15) + 0.230463210311424 × E (HMMR) + 0.21119779243939 ×E (CSPG5) + 0.271858044998993 × E (ELOVL3) + 0.0373144717068065 × E (FABP6). According to the median risk scores, 343 HCC patients were divided into a high-risk group (171 patients) and a low-risk group (172 patients). The clinical characteristics of the patients are described in Table 3. As depicted in Figures 2A, B, an increase in risk score was related to poor OS. The Kaplan–Meier curve analysis showed that patients in the high-risk group exhibited reduced OS (P <0.001, Figure 2C). Time-dependent ROC analysis was used to evaluate the prognostic evaluation ability of a six-gene signature (Figure 2D). The area under the curve (AUC) at 1-, 2-, and 3-year OS were 0.803, 0.731, and 0.699, respectively. Moreover, PCA and t-SNE analyses were used to identify the different distributions between the two risk groups. As shown in Figures 2E, F, the distribution patterns of patients in the two groups were different.

TABLE 3

Table 3 Baseline characteristics of the patients in different risk groups.

FIGURE 2

Figure 2 Prognostic analysis of the six-gene signature in the TCGA set. (A, B) The distribution and median value of the risk scores in the TCGA set. (C) Kaplan–Meier curves for the OS of patients in the high-risk group and low-risk group. (D) AUC of time-dependent ROC curves in the TCGA set. (E) PCA and (F) t-SNE analysis of the ICGC set.

Verification of the Metabolism-Related Gene Signature

According to the median risk score of the TCGA set, 231 patients from the ICGC set were divided into a high-risk group (175 patients) and a low-risk group (56 patients). Similar to the TCGA set, patients in the low-risk group exhibited better OS (Figures 3A, B). The Kaplan–Meier curve also confirmed that the OS of patients in the high-risk group was worse than that of patients in the low-risk group (P <0.001, Figure 3C). Similarly, the AUC of the six-gene signature was 0.712, 0.697, and 0.715 at the 1-, 2-, and 3-year timepoints, respectively (Figure 3D). Additionally, PCA and t-SNE also showed that the two groups were presented in two different patterns (Figures 3E, F).

FIGURE 3

Figure 3 Validation of the six-gene signature in the ICGC set. (A, B) The distribution and median value of the risk scores in the ICGC set. (C) Kaplan–Meier curves for the OS of patients in the high-risk group and low-risk group. (D) AUC of time-dependent ROC curves in the ICGC set. (E) PCA and (F) t-SNE analysis of the ICGC set.

Performance Comparison of the Metabolism-Related Gene Signature With Other Gene Signatures in Prognosis Evaluation

We further compared the prediction performance of the metabolism-related gene signature with four other published gene signatures obtained from Huo’s (12), Chen’s (13), Jiang’s (14) and Xu’s (15) studies, in the TCGA database. Figure 4A reveals that the AUC of the metabolism-related gene signature for 1-year OS was 0.803, which is significantly larger than that of Huo’s (0.789), Chen’s (0.769), Jiang’s (0.746) and Xu’s (0.688) gene signatures. The results demonstrate that the metabolism-related gene signature offers superior prognosis evaluation performance than the four previously published gene signatures for HCC.

FIGURE 4

Figure 4 The expression level of six genes in HCC. (A) The ROC analysis at 1 year of overall survival for the MyGeneSig, HuoGeneSig, ChenGeneSig, JiangGeneSig, and XuGeneSig. (B) The mRNA expression analysis by qRT-RCR. (C) The immunohistochemistry staining images of G6PD, AKR1B15, HMMR and CSPG5 from the HPA. *P<0.05; **P<0.01.

Verification of the Expression of Six Metabolism-Related Genes Between HCC and Adjacent Normal Liver Tissues

To evaluate the expression of the six metabolism-related genes (G6PD, AKR1B15, HMMR, CSPG5, ELOVL3, and FABP6) between HCC and adjacent normal liver tissues, qRT-PCR was performed to quantify mRNA expression levels. The results of qRT-PCR revealed that in comparison with adjacent normal liver tissues, the six metabolism-related genes were all upregulated in HCC tissues (Figure 4B, P <0.05). To further validate the differences in protein expression of the six genes, representative immunohistochemical images from the HPA database were obtained. G6PD, AKR1B15, HMMR, and CSPG5 expression in HCC tissues was significantly higher than in the adjacent normal liver tissues (Figure 4C). However, the expression of ELOVL 3 and FABP 6 was negative in both HCC and adjacent normal liver tissues (Figure S1).

Correlation Analysis Between Clinicopathological Characteristics and the Metabolism-Related Gene Signature

We conducted a series of correlation analyses to study the relationship between the metabolism-related gene signature and clinicopathological characteristics. The results of the heat map (Figure 5A) and scatter diagrams (Figures 5B–E) showed that tumor grade, clinical stage, T stage, and survival status were significantly related to the risk score. Next, we conducted univariate and multivariate Cox analyses to evaluate and validate whether the gene signature represents an independent risk factor for OS (Figures 5F, G). Univariate Cox regression analysis indicated that the risk score (training set: HR = 3.328, 95% CI = 2.423–4.571, P <0.001; validation set, HR = 2.518, 95% CI = 1.612–3.933, P <0.001) and tumor stage (training set: HR = 2.836, 95% CI = 1.934–4.158, P <0.001; validation set, HR = 2.492, 95% CI = 1.351–4.599, P <0.01) were closely correlated with OS. Multivariate Cox regression analysis further confirmed that the risk score was an independent risk predictor of OS (training set: HR = 2.962, 95% CI = 2.149–4.084, P <0.001; validation set: HR = 2.264, 95% CI = 1.434–3.574, P <0.001). These results indicate that the six-gene signature is an independent risk indicator for the prognostic evaluation of HCC.

FIGURE 5

Figure 5 Associations between the signature risk scores and clinicopathological characteristics. (A) Heatmap of the clinicopathological characteristics and the risk score; (B–E) The risk score in different groups divided by clinical characteristics. Univariate Cox regression analyses and multivariate Cox regression analyses regarding OS in the TCGA (F) and the ICGC (G) set. **P<0.01; ***P<0.001.

Construction and Validation of a Predictive Nomogram

Since independent prognostic analysis confirmed that tumor stage and risk score were independent risk factors for HCC, we constructed a nomogram to estimate the probability of 1-, 2-, and 3-year OS (Figure 6A). The calibration chart showed that the OS probability predicted by the nomogram approximated the actual OS probability well, especially for the 3-year OS probability (Figures 6B–D). Moreover, the C-index of the nomogram (0.735, 95% CI: 0.710–0.760) was higher than that of both stage (0.637, 95% CI: 0.609–0.665) and risk score (0.716, 95% CI: 0.692–0.741). We further constructed a DCA curve to predict the reliability of the nomogram (Figures 6E–G). The results confirmed that the nomogram provides the highest reliability, in comparison to single tumor stage and risk score, especially for predicting 3-year OS.

FIGURE 6

Figure 6 Nomograms to predict OS in hepatocellular carcinoma patients. (A) Nomograms using clinical traits shared between the TCGA set. The calibration and DCA curve for determining the reliability of the nomogram to predict the 1-year (B, E), 2-year (C, F) and 3-year (D, G) OS.

Functional Analyses

To elucidate the molecular mechanism associated with the six-gene signature, GSEA was applied to the training set. As we have seen, the metabolism-related KEGG pathways were significantly enriched in the high-risk group (Figures 7A, B). These KEGG pathways include amino sugar and nucleotide sugar metabolism (FDR = 0.013, P = 0.002), glutathione metabolism (FDR = 0.311, P = 0.024), inositol phosphate metabolism (FDR = 0.003, P <0.000), purine metabolism (FDR = 0.002, P <0.000), pyrimidine metabolism (FDR = 0.002, P <0.000), and selenoamino acid metabolism (FDR = 0.028, P = 0.004). Interestingly, many immune-related responses in the high-risk group were also significantly enriched, including antigen processing and presentation (FDR = 0.032, P = 0.040), B cell receptor signaling pathway (FDR = 0.040, P = 0.026), chemokine signaling pathway (FDR = 0.029, P = 0.016), leukocyte transendothelial migration (FDR = 0.038, P = 0.018), natural killer cell mediated cytotoxicity (FDR = 0.017, P = 0.010), and T cell receptor signaling pathway (FDR = 0.027, P = 0.010).

FIGURE 7

Figure 7 Functional analyses of the low-risk and high-risk groups. The differences of the metabolism-related (A) and immune-related (B) pathways between high-risk and low-risk groups. The scores of 16 immune cells (C) and 13 immune-related functions (D) between the two risk groups. Adjusted P values were showed as: ns, not significant; *P <0.05; **P <0.01; ***P <0.001.

An increasing number of studies have indicated that metabolic reprogramming of tumors is related to immune infiltration and immune responses (16). Therefore, we further investigated whether there were different subgroups of immune cells in the two risk groups and scored the related functions by using ssGSEA. The scores of the immune cells (including aDCs, DCs, iDCs, macrophages, mast cells, NK cells, Tfh cells, Th1 cells, Th2 cells, and Treg cells) were significantly different between the two risk groups (P <0.05, Figure 7C). Additionally, the high-risk group showed higher scores for APC_co_stimulation, APC_co_inhibition, CCR, Check-point, HLA, MHC_class_I, Parainflammation, and T_cell_co-inhibition, while the score for Type_II_IFN_Response exhibited the opposite trend (P <0.05). The ICGC set study confirmed the differences between the two risk groups in aDC cells, DC cells, iDC cells, macrophages cells, NK cells, TH2 cells, Treg cells, APC_co_stimulation, check-point, and HLA (P <0.05, Figure 7D).

Relationship Between Metabolism-Related Gene Signature and Expression of Immune Checkpoints

Immune-based therapies have become a systematic treatment method for improving the prognosis of advanced cancers (17). Therefore, we further analyzed the correlation between the risk groups and the expression levels of immune checkpoint proteins. As shown in Figures 8A–F, in comparison with the low-risk group, the high-risk group had significantly higher expression levels of programmed death-1 (PD-1), programmed death ligand 1 (PD-L1), cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), lymphocyte activation gene-3 (LAG3), T-cell immunoglobulin and mucin-domain containing-3 (TIM-3), and T cell immunoreceptor with Ig and ITIM domains (TIGIT). Moreover, recent studies have indicated that TGF-β signaling and a cohort of 30 extracellular matrix genes (C-ECM) are significantly associated with cancer immunosuppression and poor prognosis. Figures 8G, H describe the increased expression of these genes in the high-risk group.

FIGURE 8

Figure 8 Expression of immune checkpoints (A–F), TGF-β–encoding (G) and C-ECM (H) signature genes.

Discussion

HCC has become a major health concern and threat to global mortality, especially in China. From 1991 to 2011, the mortality rate associated with liver cancer has continuously increased (18). According to statistics, in 2015, there were approximately 370,000 new cases of liver cancer, and 326,000 related deaths (19). Progress in radical surgery plays a vital role in alleviating the global burden of HCC. However, the long-term outcome of HCC remains unclear. Therefore, it is essential to establish an effective prognostic signature for the evaluation and treatment of HCC.

In this study, we used the transcriptome profiling and the corresponding clinical information of HCC patients, obtained from the TCGA database, to identify a six metabolism-related gene signature (G6PD, AKR1B15, HMMR, CSPG5, ELOVL3, FABP6) for the prognostic evaluation of HCC, and verified them in the ICGC database. G6PD is a rate-limiting enzyme that catalyzes the pentose phosphate pathway and participates in regulating the redox homeostasis of cells exposed to oxidative damage (20). Abnormal behavior of G6PD is related to a variety of pathological processes and diseases, including inflammation (21), diabetes (22), and tumors (23). High levels of G6PD expression were observed in human livers infected with HBV and HBV-related cancers, which provides further confirmation of the findings presented here (24). AKR1B15, a newly discovered Aldo-keto reductase (AKR), shares 92% amino acid sequence identity with AKR1B10. Studies have shown that AKR1B10 can induce a variety of cancers, such as liver cancer (25), non-small cell lung cancer (26), and pancreatic cancer (27), and represents a promising potential cancer target. Presently, a few studies have indicated that mutation of the AKR1B15 allele is related to mitochondrial oxidative phosphorylation and serous ovarian cancer, which has attracted increasing attention from researchers (28, 29). HMMR is a protein that regulates cell growth and is involved in maintaining homeostasis and the directional regulation of mitotic and meiotic spindles (30). It has been shown that upregulation of HMMR expression is related to the aggressive growth and low survival rate of various cancers, such as breast cancer (31), colorectal cancer (32), and gastric cancer (33). Several other types of cancer, such as breast cancer (34) and malignant peripheral nerve sheath tumors (35), have been found to exhibit poor patient survival, which has been correlated with low expression of HMMR. ELOVL3 is an ultra-long-chain fatty acid elongase that participates in the synthesis of fatty acids (36). A small number of studies have shown that ELOVL3 is highly expressed in colorectal and prostate cancers (37, 38). FABP6 is a protein involved in bile acid digestion, absorption, metabolism and enterohepatic circulation (39). Recent studies have indicated that elevated levels of FABP6 are associated with the initiation and development of colorectal cancer by regulating the NF-κB pathway (40). The relationship between these six genes and the progression of HCC in patients remains to be clarified, as very few studies have examined these genes. In this study, these six genes were all found to be upregulated in HCC tissues and were significantly correlated with a short OS (P <0.001).

Furthermore, according to the risk scores calculated based on the expression levels and regression coefficient values of the six metabolism-related genes, we divided the patients into two risk groups. The low-risk group exhibited better OS than the high-risk group. Thus, the six-gene signature also represents an independent prognostic risk indicator for HCC. Based on multivariate Cox analysis, we proposed a nomogram based on tumor stage and risk score. The C-index and DCA curve proved that the nomogram was superior to both tumor stage and single risk score in predicting tumor prognosis, especially for predicting OS in 3 years.

In addition, metabolism-related pathways were identified in the high-risk group. Interestingly, patients in the high-risk group were also more closely associated with changes in immune-related pathways and upregulated expression of immune checkpoint proteins (PD-1, PD-L1, TIM3, TIGIT), TGF-β, and C-ECM genes. Immune checkpoint blocking therapy is promising in the treatment of HCC (41). It has been suggested that HCC patients with higher expression of PD-1 or PD-L1 might be more likely to respond to immune checkpoint blocking therapy (42). Therefore, maybe we could identify high-risk patients early based on our gene signature and conduct immune checkpoint blocking treatment or other applicable therapy to improve the prognosis of patients, which needs further verification in the future. Previous researches have indicated that an increase in tumor-related macrophages (43), NK cells (44), and Treg cells (45) in the tumor microenvironment is related to reduced OS in HCC patients, and strategies for immunotherapy have been proposed. However, it is worth noting that the relationship between metabolism-related genes and immunity remains to be clarified.

This study has several limitations. Our risk signature significantly stratifies HCC patients, which makes the prognosis and immunotherapy response evaluation accurate and reliable. However, our research is mainly based on a public database, it is retrospective, and belongs to a small sample study. Additional studies using large-scale, prospective, and multicenter clinical trials are needed to validate the robustness and reproducibility of the gene signature. Although we have performed a comprehensive bioinformatics analysis to construct and verify the prognostic metabolism-related gene signature in HCC, it may not be as accurate for different types of HCC. Furthermore, the underlying molecular mechanism of how these six genes affect the development of HCC remains unclear, and further experiments are needed. The correlation between risk score and immune status has not been verified by basic experiments, which is an important issue worthy of further study.

Conclusion

In conclusion, this study developed a prognostic signature comprised of six metabolism-related genes for HCC. The signature was confirmed to be an independent risk indicator related to the OS of HCC in the TCGA and ICGC databases. Moreover, we constructed a nomogram composed of risk score and tumor grade, which provided higher accuracy in predicting the OS of patients with HCC. However, longitudinal clinical experiments should be conducted to verify this hypothesis.

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 authors.

Author Contributions

CY and MY conceived the study, performed the bioinformatics analyses and wrote the manuscript. MC and JO performed the experiments. WT and FD downloaded and organized the clinical and gene expression data. DY, SL and YZ performed the statistical analyses. CZ and YC critically revised the article for essential intellectual content and administrative support. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Nature Science Foundation of China (82071655, 81860276), Scientific Research Project Foundation of Hubei Provincial Health Commission (WJ2019M179).

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.

Acknowledgments

We thank the TCGA and ICGC databases for their generous sharing of large amounts of data. We apologize for any omission of citations and references due to space limitations.

Supplementary Material

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

References

1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2018) 68:394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Jemal A, Bray F, Center MM, Ferlay J, Ward E, Forman D. Global Cancer Statistics. CA Cancer J Clin (2011) 61:69–90. doi: 10.3322/caac.20107

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Siegel R, Naishadham D, Jemal A. Cancer Statistics, 2013. CA Cancer J Clin (2013) 63:11–30. doi: 10.3322/caac.21166

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Marrero JA, Feng Z, Wang Y, Nguyen MH, Befeler AS, Roberts LR, et al. Alpha-Fetoprotein, Des-Gamma Carboxyprothrombin, and Lectin-Bound Alpha-Fetoprotein in Early Hepatocellular Carcinoma. Gastroenterology (2009) 137:110–8. doi: 10.1053/j.gastro.2009.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Hanazaki K, Kajikawa S, Koide N, Adachi W, Amano J. Prognostic Factors After Hepatic Resection for Hepatocellular Carcinoma With Hepatitis C Viral Infection: Univariate and Multivariate Analysis. Am J Gastroenterol (2001) 96:1243–50. doi: 10.1111/j.1572-0241.2001.03634.x

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Vander Heiden MG, DeBerardinis RJ. Understanding the Intersections Between Metabolism and Cancer Biology. Cell (2017) 168:657–69. doi: 10.1016/j.cell.2016.12.039

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Kreuzaler P, Panina Y, Segal J, Yuneva M. Adapt and Conquer: Metabolic Flexibility in Cancer Growth, Invasion and Evasion. Mol Metab (2020) 33:83–101. doi: 10.1016/j.molmet.2019.08.021

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lacroix M, Riscal R, Arena G, Linares LK, Le Cam L. Metabolic Functions of the Tumor Suppressor P53: Implications in Normal Physiology, Metabolic Disorders, and Cancer. Mol Metab (2020) 33:2–22. doi: 10.1016/j.molmet.2019.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Tarasenko TN, McGuire PJ. The Liver is a Metabolic and Immunologic Organ: A Reconsideration of Metabolic Decompensation Due to Infection in Inborn Errors of Metabolism (IEM). Mol Genet Metab (2017) 121:283–8. doi: 10.1016/j.ymgme.2017.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ward PS, Thompson CB. Metabolic Reprogramming: A Cancer Hallmark Even Warburg did Not Anticipate. Cancer Cell (2012) 21:297–308. doi: 10.1016/j.ccr.2012.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Chen W, Tang D, Ou M, Dai Y. Mining Prognostic Biomarkers of Hepatocellular Carcinoma Based on Immune-Associated Genes. DNA Cell Biol (2020) 39:499–512. doi: 10.1089/dna.2019.5099

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Huo J, Wu L, Zang Y, Dong H, Liu X, He F, et al. Eight-Gene Metabolic Signature Related With Tumor-Associated Macrophages Predicting Overall Survival for Hepatocellular Carcinoma. BMC Cancer (2021) 21:31. doi: 10.1186/s12885-020-07734-z

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Chen Q, Li F, Gao Y, Xu G, Liang L, Xu J. Identification of Energy Metabolism Genes for the Prediction of Survival in Hepatocellular Carcinoma. Front Oncol (2020) 10:1210. doi: 10.3389/fonc.2020.01210

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Jiang N, Zhang X, Qin D, Yang J, Wu A, Wang L, et al. Identification of Core Genes Related to Progression and Prognosis of Hepatocellular Carcinoma and Small-Molecule Drug Predication. Front Genet (2021) 12:608017. doi: 10.3389/fgene.2021.608017

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xu W, Chen Z, Liu G, Dai Y, Xu X, Ma D, et al. Identification of a Potential Ppar-Related Multigene Signature Predicting Prognosis of Patients With Hepatocellular Carcinoma. PPAR Res (2021) 2021:6642939. doi: 10.1155/2021/6642939

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yang M, Chen G, Gao K, Wang Y. Tumor Immunometabolism Characterization in Ovarian Cancer With Prognostic and Therapeutic Implications. Front Oncol (2021) 11:622752. doi: 10.3389/fonc.2021.622752

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Ruf B, Heinrich B, Greten TF. Immunobiology and Immunotherapy of HCC: Spotlight on Innate and Innate-Like Immune Cells. Cell Mol Immunol (2021) 18:112–27. doi: 10.1038/s41423-020-00572-w

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Shao HY, Peng MJ, Chen G, Yan H, Li S, Treatment. Age-Period-Cohort Modeling Study on the Mortality of Hepatic Carcinoma From 1991 to 2011. Chin J Cancer Prevention Treatment (2016) 23:1465–9. doi: 10.16073/j.cnki.cjcpt.2016.22.001

CrossRef Full Text | Google Scholar

19. Cao M, Li H, Sun D, Chen W. Cancer Burden of Major Cancers in China: A Need for Sustainable Actions. Cancer Commun (Lond) (2020) 40:205–10. doi: 10.1002/cac2.12025

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yang HC, Wu YH, Yen WC, Liu HY, Hwang TL, Stern A, et al. The Redox Role of G6PD in Cell Growth, Cell Death, and Cancer. Cells (2019) 8:1055. doi: 10.3390/cells8091055

CrossRef Full Text | Google Scholar

21. Peiro C, Romacho T, Azcutia V, Villalobos L, Fernandez E, Bolanos JP, et al. Inflammation, Glucose, and Vascular Cell Damage: The Role of the Pentose Phosphate Pathway. Cardiovasc Diabetol (2016) 15:82. doi: 10.1186/s12933-016-0397-2

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Ham M, Choe SS, Shin KC, Choi G, Kim JW, Noh JR, et al. Glucose-6-Phosphate Dehydrogenase Deficiency Improves Insulin Resistance With Reduced Adipose Tissue Inflammation in Obesity. Diabetes (2016) 65:2624–38. doi: 10.2337/db16-0060

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Cai T, Kuang Y, Zhang C, Zhang Z, Chen L, Li B, et al. Glucose-6-phosphate Dehydrogenase and NADPH Oxidase 4 Control STAT3 Activity in Melanoma Cells Through a Pathway Involving Reactive Oxygen Species, c-SRC and SHP2. Am J Cancer Res (2015) 5:1610–20.

PubMed Abstract | Google Scholar

24. Liu B, Fang M, He Z, Cui D, Jia S, Lin X, et al. Hepatitis B Virus Stimulates G6PD Expression Through HBx-mediated Nrf2 Activation. Cell Death Dis (2015) 6:e1980. doi: 10.1038/cddis.2015.322

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Tsuzura H, Genda T, Sato S, Murata A, Kanemitsu Y, Narita Y, et al. Expression of Aldo-Keto Reductase Family 1 Member b10 in the Early Stages of Human Hepatocarcinogenesis. Int J Mol Sci (2014) 15:6556–68. doi: 10.3390/ijms15046556

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Fukumoto S, Yamauchi N, Moriguchi H, Hippo Y, Watanabe A, Shibahara J, et al. Overexpression of the Aldo-Keto Reductase Family Protein AKR1B10 is Highly Correlated With Smokers’ non-Small Cell Lung Carcinomas. Clin Cancer Res (2005) 11:1776–85. doi: 10.1158/1078-0432.CCR-04-1238

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Chung YT, Matkowskyj KA, Li H, Bai H, Zhang W, Tsao MS, et al. Overexpression and Oncogenic Function of Aldo-Keto Reductase Family 1B10 (AKR1B10) in Pancreatic Carcinoma. Mod Pathol (2012) 25:758–66. doi: 10.1038/modpathol.2011.191

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Calvo SE, Compton AG, Hershman SG, Lim SC, Lieber DS, Tucker EJ, et al. Molecular Diagnosis of Infantile Mitochondrial Disease With Targeted Next-Generation Sequencing. Sci Transl Med (2012) 4:118ra10. doi: 10.1126/scitranslmed.3003310

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Yashin AI, Wu D, Arbeev KG, Ukraintseva SV. Joint Influence of Small-Effect Genetic Variants on Human Longevity. Aging (Albany NY) (2010) 2:612–20. doi: 10.18632/aging.100191

PubMed Abstract | CrossRef Full Text | Google Scholar

30. He Z, Mei L, Connell M, Maxwell CA. Hyaluronan Mediated Motility Receptor (Hmmr) Encodes an Evolutionarily Conserved Homeostasis, Mitosis, and Meiosis Regulator Rather Than a Hyaluronan Receptor. Cells (2020) 9:819. doi: 10.3390/cells9040819

CrossRef Full Text | Google Scholar

31. Assmann V, Gillett CE, Poulsom R, Ryder K, Hart IR, Hanby AM. The Pattern of Expression of the Microtubule-Binding Protein RHAMM/IHABP in Mammary Carcinoma Suggests a Role in the Invasive Behaviour of Tumour Cells. J Pathol (2001) 195:191–6. doi: 10.1002/path.941

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Zlobec I, Baker K, Terracciano LM, Lugli A. Rhamm, p21 Combined Phenotype Identifies Microsatellite Instability-High Colorectal Cancers With a Highly Adverse Prognosis. Clin Cancer Res (2008) 14:3798–806. doi: 10.1158/1078-0432.CCR-07-5103

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Li H, Guo L, Li JW, Liu N, Qi R, Liu J. Expression of Hyaluronan Receptors CD44 and RHAMM in Stomach Cancers: Relevance With Tumor Progression. Int J Oncol (2000) 17:927–32. doi: 10.3892/ijo.17.5.927

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Blanco I, Kuchenbaecker K, Cuadras D, Wang X, Barrowdale D, de Garibay GR, et al. Assessing Associations Between the AURKA-HMMR-TPX2-TUBG1 Functional Module and Breast Cancer Risk in BRCA1/2 Mutation Carriers. PloS One (2015) 10:e0120020. doi: 10.1371/journal.pone.0120020

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Mantripragada KK, Spurlock G, Kluwe L, Chuzhanova N, Ferner RE, Frayling IM, et al. High-Resolution DNA Copy Number Profiling of Malignant Peripheral Nerve Sheath Tumors Using Targeted Microarray-Based Comparative Genomic Hybridization. Clin Cancer Res (2008) 14:1015–24. doi: 10.1158/1078-0432.CCR-07-1305

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Guillou H, Zadravec D, Martin PG, Jacobsson A. The Key Roles of Elongases and Desaturases in Mammalian Fatty Acid Metabolism: Insights From Transgenic Mice. Prog Lipid Res (2010) 49:186–99. doi: 10.1016/j.plipres.2009.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Alnabulsi A, Cash B, Hu Y, Silina L, Alnabulsi A, Murray GI. The Expression of Brown Fat-Associated Proteins in Colorectal Cancer and the Relationship of Uncoupling Protein 1 With Prognosis. Int J Cancer (2019) 145:1138–47. doi: 10.1002/ijc.32198

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Yang Y, Liu L, Li M, Cheng X, Fang M, Zeng Q, et al. The Chromatin Remodeling Protein BRG1 Links ELOVL3 Trans-Activation to Prostate Cancer Metastasis. Biochim Biophys Acta Gene Regul Mech (2019) 1862:834–45. doi: 10.1016/j.bbagrm.2019.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zwicker BL, Agellon LB. Transport and Biological Activities of Bile Acids. Int J Biochem Cell Biol (2013) 45:1389–98. doi: 10.1016/j.biocel.2013.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Fang C, Dean J, Smith JW. A Novel Variant of Ileal Bile Acid Binding Protein is Up-Regulated Through Nuclear Factor-Kappab Activation in Colorectal Adenocarcinoma. Cancer Res (2007) 67:9039–46. doi: 10.1158/0008-5472.CAN-06-3690

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Cheng AL, Hsu C, Chan SL, Choo SP, Kudo M. Challenges of Combination Therapy With Immune Checkpoint Inhibitors for Hepatocellular Carcinoma. J Hepatol (2020) 72:307–19. doi: 10.1016/j.jhep.2019.09.025

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Macek Jilkova Z, Aspord C, Decaens T. Predictive Factors for Response to PD-1/PD-L1 Checkpoint Inhibition in the Field of Hepatocellular Carcinoma: Current Status and Challenges. Cancers (Basel) (2019) 11:1554. doi: 10.3390/cancers11101554

CrossRef Full Text | Google Scholar

43. Zhang Q, He Y, Luo N, Patel SJ, Han Y, Gao R, et al. Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma. Cell (2019) 179:829–45.e20. doi: 10.1016/j.cell.2019.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Cai L, Zhang Z, Zhou L, Wang H, Fu J, Zhang S, et al. Functional Impairment in Circulating and Intrahepatic NK Cells and Relative Mechanism in Hepatocellular Carcinoma Patients. Clin Immunol (2008) 129:428–37. doi: 10.1016/j.clim.2008.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Zhou SL, Zhou ZJ, Hu ZQ, Huang XW, Wang Z, Chen EB, et al. Tumor-Associated Neutrophils Recruit Macrophages and T-Regulatory Cells to Promote Progression of Hepatocellular Carcinoma and Resistance to Sorafenib. Gastroenterology (2016) 150:1646–58.e17. doi: 10.1053/j.gastro.2016.02.040

PubMed Abstract | CrossRef Full Text | Google Scholar