A novel metabolism-related gene signature in patients with hepatocellular carcinoma

Introduction

The incidence of liver cancer is increasing worldwide (Llovet et al., 2021; Villanueva, 2019). Liver cancer is the sixth most common malignancy and the third leading cause of cancer-related death globally. According to the latest statistics, there are 41,210 new cases of liver cancer and 29,380 deaths due to liver cancer in the United States in 2022 (Siegel et al., 2023). Hepatocellular carcinoma (HCC) patients account for more than 90% of liver cancer cases, and with the implementation of targeted and immunotherapy, the life expectancy of patients with hepatocellular carcinoma has increased (Llovet et al., 2018; Zucman-Rossi et al., 2015; Schulze et al., 2015). These treatments have a good therapeutic effect in patients with early-stage hepatocellular carcinoma. Due to chemotherapy tolerance and insensitivity to radiotherapy, the choice of effective treatments for patients with advanced hepatocellular carcinoma is seriously affected, resulting in poor prognosis (Llovet, Burroughs & Bruix, 2003). In the United States, the 5-year survival of patients with HCC is only 18% (Jemal et al., 2017). Traditional prognosis prediction systems for cancer patients, such as TNM staging, are increasingly difficult to cover the diversity of clinical features of HCC patients (Shao et al., 2016). The development of new prognostic models at multiple levels will help to better distinguish between different types of HCC patients (Shibata, 2021). Recent studies have demonstrated favorable outcomes in discriminating between HCC and predicting HCC prognosis, through the creation of multiple test-based indicators (Luo et al., 2022; Xie et al., 2022). This facilitates better extraction of features from different types of patients to facilitate more effective treatment of different types of patients.

Abnormal metabolism is one of the characteristics of HCC (Hanahan & Weinberg, 2011). During the malignant transformation of cells, metabolism usually undergoes drastic changes. Tumor cells regulate metabolism, promoting energy production and accumulation (Cheng et al., 2018). Hepatocellular carcinoma cells have elevated levels of metabolism to maintain their high metabolic levels for plasma membrane synthesis and energy production (Sangineto et al., 2020). Oncogenic signaling pathways, including B-Raf kinase (BRAF) and epidermal growth factor receptor (EGFR), drive dysregulation of fatty acid metabolism, affecting membrane composition and saturation to regulate tolerance to reactive oxygen species (ROS) and cancer cell survival (Talebi et al., 2018; Gimple et al., 2019; Bi et al., 2019). In addition, recent studies have shown that alterations in metabolism in tumor cells are able to promote tumor invasion and metastasis, regulate oxidative stress, and provide energy under various cellular stress situations (Pascual et al., 2017; Rohrig & Schulze, 2016; Ladanyi et al., 2018). Given that metabolism abnormalities play an important role in tumor progression, it has not been explored whether metabolism genes can be analyzed to build a prognostic model for HCC patients.

The role of metabolism in the regulation of immune cells has recently attracted widespread attention. Research evidence in several types of solid tumors shows the importance of metabolic reprogramming of immune cells in tumors, suggesting a new strategy for the treatment of HCC patients (Zhang et al., 2018). The function of immune cells in the HCC tumor microenvironment is closely related to abnormal metabolism (Hu et al., 2020). However, whether the expression level of MRGs in HCC patients affects immune cell infiltration in the tumor microenvironment remains to be explored.

In our research, we developed a novel prognostic prediction model based on MRGs, which is named “Metabolism-Related Risk Score” (MRRS), to predict the prognosis and survival of HCC patients. The immune-related components of different risk groups were analyzed, and the results showed that there were significant differences in the tumor immune microenvironment between the two groups. The evaluation of our predictive typing model shows great potential to guide the classification and treatment of HCC patients, demonstrating its clinical value and significance. Currently, a variety of genes and proteins that regulate metabolism have been identified, including GOT2, which has been reported to be involved in metabolism (Liu et al., 2022). However, the effect of this metabolism-related gene on the development and progression of hepatocellular carcinoma is unclear. Our study preliminarily explored the relationship between GOT2 and the prognosis of patients with hepatocellular carcinoma and preliminarily explored its effect on the migration ability of hepatocellular carcinoma.

Methods

Collection and processing of gene expression data and clinical information in HCC patients

Gene expression data and corresponding clinical information of 374 tumor tissues and 50 normal tissues of 374 patients with hepatocellular carcinoma as of March 03, 2022, were obtained from TCGA (https://portal.gdc.cancer.gov/repository). Of these, 44 patients were removed due to incomplete clinical data. Therefore, tissue gene expression data and complete follow-up information from the remaining patients (n = 330) were incorporated into our training dataset for further analysis. The test dataset used for validation was gene expression data and clinical information from an additional 231 tumor samples obtained from the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP). These samples were mainly from Japan (Fujimoto et al., 2016). Since this data is all available online and comes with permission to use it, no additional ethical approval is required. The current study follows the TCGA and ICGC data access policies and publication guidelines. The flowchart of this research is shown in Fig. 1.

Result

Analysis of the differential expression and prognostic relationship of DEGs in HCC based on TCGA database

By using univariate COX regression analysis, we found that 79 DEGs in tumor tissues of HCC patients were associated with the prognosis of HCC patients, of which 21 genes with high expression indicated a better prognosis and the other 58 genes were associated with a poor prognosis (Fig. 3A). Of these, 75 genes were MRGs (Fig. 3B). Among the 75 MRGs associated with the prognosis of HCC patients, 55 genes were highly expressed in tumor tissues of HCC patients, and 20 genes were poorly expressed. We used heatmaps for illustration (Fig. 3C). Then, we evaluated the protein-protein interactions of these genes, further revealing the strong interaction activity of these molecules at the protein level, as shown in Fig. 2D. Similarly, correlated networks built based on mRNA expression levels in TCGA showed a negative (blue) and positive (red) correlation between these prognostically relevant MRGs, as shown in Fig. 3E.

Development and evaluation of metabolism-related risk score

The expression profiles of the above 75 genes were used to establish a prognostic model by LASSO Cox regression analysis. Based on the optimal value of λ, 14 genes that contribute the most to the prognosis of HCC patients were identified (Table 2). Among these genes, except for ACADS, GOT2, ADH4, and CSAD, high expression of the other 10 genes was associated with a poor prognosis in HCC patients (P < 0.05). The risk score is calculated as follows: risk score = (0.0875 * DLAT expression) + (0.2953 * SEPHS1 expression) + (-0.1116 * ACADS expression) + (0.1978 * UCK2 expression) + (-0.0143 * GOT2 expression) + (-0.0295 * ADH4 expression) + (-0.3244 * LDHA expression) + (0.0520 * ME1 expression) + (0.0105 * TXNRD1 expression) + (0.0433 * B4GALT2 expression) + (0.1975* AK2 expression) + (0.1783* PTDSS2 expression) + (-0.023 * CSAD expression) + (0.0207 * AMD1 expression). Based on the median of the risk score, the sample was rated as low and high risk (Fig. 4A). Of these, 164 patients were in the high-risk group and 165 patients were in the low-risk group. Principal component analysis (PCA) and t-distribution random neighborhood embedding (t-SNE) analysis showed that patients in different risk groups were distributed in two directions (Figs. 4B–4C). This suggests that the MRRS model can distinguish well between high- and low-risk patients. It is shown that the patients in the high-risk group had a higher probability of death than those in the low-risk group (Fig. 4D). In addition, Kaplan–Meier survival analysis also confirmed that the survival time of the high-risk group was significantly reduced compared to the low-risk group (Fig. 4E, P < 0.001). The predictive performance of the model was evaluated using a time-varying ROC curve, with the area under the curve (AUC) reaching 0.829 at 1 year, 0.770 at 2 years, and 0.760 at 3 years (Fig. 4F). Overall, these findings suggest that MRRS was able to distinguish the prognosis of HCC patients in the TCGA database.

Table 2:

The genes that contribute to the prognosis of HCC patients.

Based on the optimal value of λ, 14 genes that contribute the most to the prognosis of HCC patients were identified.

Gene	Coef
DLAT	0.087527704
SEPHS1	0.295255334
ACADS	−0.111567502
UCK2	0.197793632
GOT2	−0.014390167
ADH4	−0.029554667
LDHA	0.324459855
ME1	0.052008356
TXNRD1	0.010457298
B4GALT2	0.043281636
AK2	0.19745112
PTDSS2	0.178350248
CSAD	−0.023302247
AMD1	0.02060711

DOI: 10.7717/peerj.16335/table-2

Validation of MRRS in the ICGC cohort

To test the robustness of the model developed by the TCGA cohort, we validated the model with HCC patients’ data from the ICGC database. We calculated the risk score for each patient in the ICGC cohort using the same formula. Based on the above median, we divided HCC patients from the ICGC cohort into high-risk (n = 214) and low-risk (n = 13) (Fig. 5A). In the ICGC cohort, the results of PCA and t-SNE analyses were similar to those of TCGA (Figs. 5B–5C). The high-risk group had a worse prognosis than the low-risk group (Fig. 5D). The survival curve showed that MRRS was effective in differentiating patients with different prognostic conditions, and patients in the high-risk group tended to have a higher risk of death and significantly shorter survival times than those in the low-risk group (Fig. 5E, P < 0.05). ROC analysis showed that MRRS had good predictive power on the prognosis of patients in the ICGC cohort (AUC = 0.757, 0.738, and 0.725; at 1, 2, and 3 years, respectively in ICGC, Fig. 5F). Taken together, this evidence shows that the model was able to predict the prognosis of HCC patients well.

Independent prognostic value of MRRS in patients with HCC

To further determine whether MRRS can be used as an independent predictor of prognosis in HCC patients, we first performed univariate and multivariate Cox regression analysis in the TCGA cohort. As shown in Figs. 6A–6B, MRRS was found to be significantly correlated with the prognosis of HCC patients in both univariate Cox regression analyses (Hazard ratio (HR) = 3.757, 95% CI [2.819–5.007], p < 0.001) and multivariate Cox regression analysis (HR = 3.695, 95% CI [2.695–5.066], p < 0.001). Univariate and multivariate Cox regression analyses were then performed in the ICGC cohort for validation. The results showed that MRRS remained an independent predictor of prognosis in the ICGC cohort (HR = 3.023, 95% CI [1.958–4.667], p < 0.001, and HR = 2.574, 95% CI [1.602–4.135], p < 0.001; Figs. 6C–6D). In addition, the multivariate analysis to remove confounding factors, age, sex, and stage could not be used as independent prognostic indicators for HCC patients, showing the superiority of MRRS in the prognosis assessment of HCC patients. Overall, our study suggests that MRRS can serve as an independent prognostic factor for patients with HCC.

GO and KEGG analysis of DEGs in high and low-risk groups in the TCGA and ICGC cohorts

To elucidate the biological functions and pathways associated with MRRS risk scores, we first used GO enrichment analysis and KEGG pathway analysis for DEGs between high and low-risk groups in the TCGA cohort. The results of GO enrichment analysis showed that in terms of Biological Processes, signal transduction, cell division, and apoptotic process were significantly enriched. In terms of cellular components, cytosol, cytoplasm, and plasma membrane were significantly enriched; In terms of Molecular Function, protein binding, identical protein binding, and ATP binding were significantly enriched (Fig. 7A). KEGG enrichment results showed significant enrichment of metabolic pathways such as drug metabolism-cytochrome P450, retinol metabolism, central carbon metabolism, and tyrosine metabolism. Cell cycle, glycolysis/gluconeogenesis, ECM-receptor interaction, and signaling pathways such as the HIF-1 signaling pathway and PPAR signaling pathway were also enriched (Fig. 7B). In addition, immune-related signaling pathways such as the IL-17 signaling pathway and the TNF signaling pathway were significantly affected (Fig. 7B). Then we performed GO enrichment analysis and KEGG enrichment analysis for DEGs between high and low-risk groups in the ICGC cohort. GO enrichment analysis showed significant enrichment of cell adhesion, cell division, and positive regulation of T cell proliferation and activation pathways (Fig. 7C). KEGG enrichment results showed that in addition to the significant enrichment of metabolic signaling pathways, a variety of immune-related signaling pathways were also significantly enriched, such as Th17 cell differentiation, human T-cell leukemia virus 1 infection, Th1 and Th2 cell differentiation, leukocyte transendothelial migration, antigen processing and presentation, Fc gamma R-mediated phagocytosis, and B cell receptor signaling pathway (Fig. 7D). In addition, consistent with the results of the analysis in the TCGA cohort, the PPAR signaling pathway and PI3K-Akt signaling pathway were also significantly affected (Fig. 7D). Taken together, these findings highlight that HCC patients in the high and low-risk groups may lead differences in tumor malignant characterization through these pathways, ultimately leading to different prognoses of HCC patients.

Infiltration of immune cells in TME in high- and low-risk populations

GO and KEGG enrichment analysis showed differences in function and pathways between risk groups in ECM-receptor interaction pathways and multiple immune-related signaling pathways. To further explore the correlation between different risk scores and immune status, we further explored the content of immune matrix components in the tumor microenvironment. According to the ESTIMATE algorithm, B cells, CD4+ T cells, CD8+ T cells, DC cells, macrophages, and neutrophil subsets in tumor tissues of high-risk groups were significantly upregulated significantly higher in the high-risk group (P < 0.05, Fig. 8), indicating that there are more immune cell components in TME in high-risk HCC patients.

Downregulation of GOT2 promotes the migration capacity of hepatocellular carcinoma

By analyzing the relationship between prognosis and GOT2 mRNA expression levels in patients with hepatocellular carcinoma in the TCGA database, we found that patients with low GOT2 expression had a worse prognosis in the TCGA database (Fig. 9A). To further explore the biological significance of GOT2 in the progression of hepatocellular carcinoma, we first transfected siRNAs targeting GOT2 (siGOT2) or negative control siRNA (siNC) in HEK293 cells and verified their knockdown efficiency by RT-qPCR (Fig. 9B). The wound healing assay was used to assess the effect of GOT2 on the migration of hepatocellular carcinoma cells (Huh7 and MHCC97H). As shown in Figs. 9C–9F, down-regulation of GOT2 significantly inhibited migration in both cell lines compared to the control group. Together, our findings provide new insights into the role of GOT2 in influencing malignant phenotypes by regulating migration in hepatocellular carcinoma cells.

Discussion

HCC is an extremely common malignancy worldwide, it is important to develop reliable prognostic indicators for HCC patients. In this study, we established 14 genetic biomarkers as new prognostic models of MRRS and analyzed their ability to predict prognosis in HCC patients. The prognostic performance of the model is verified by the survival curve and ROC curve, and the results show that the model has good predictive performance. Its prediction efficiency was verified in the ICGC cohort. MRRS is a good independent indicator for predicting the prognosis of HCC patients in both TCGA and ICGC cohorts. Our prognostic model can help predict the prognosis of HCC patients clinically, thereby recommending better treatment measures for high-risk HCC patients.

GO and KEGG analysis showed that cell proliferation signals were significantly altered in patients in the high-risk group, suggesting that metabolism significantly affects cell proliferation, and cell proliferation may be the cause of poor prognosis in HCC patients. The mechanisms involved may be related to the PPAR signaling pathway, and PI3K-Akt signaling pathway, so PPAR inhibitors and PI3K-Akt inhibitors can be tried in high-risk patients to have a better prognosis for HCC patients. Functional analysis showed that immune signals were widely involved in tumor processes in high-risk patients, and more immune cells were infiltrated in the tumor microenvironment of high-risk patients. Our results are consistent with previous studies, and that increased expression of tumor blasts such as regulatory T cells, tumor-associated macrophages, tumor-associated neutrophils, and myeloid-derived suppressor cells in most cancers generally predicts worse outcomes (Senovilla et al., 2012). TME is generally divided into three categories (Chen & Mellman, 2017; Lanitis et al., 2017): (1) immune-inflamed: immune cells exist near tumor cells, (2) immune-excluded: immune cells exist around the stroma but do not penetrate the tumor, (3) immune-desert: lack of immune cell infiltration. In the current study, we reasonably speculate that tumors in patients with high-risk HCC may be Immune-excluded tumors (IETs). In this case, although TME shows abundant immune cell infiltration, cytotoxic T lymphocytes (CTL) cannot effectively infiltrate the tumor and exert a killing effect. This suggests that immunotherapy alone in high-risk patients may be less effective in patients. At present, the possibility that high-risk patients are immune-inflamed cannot be ruled out, and further evaluation of immunotherapy responses between different patients is needed to obtain better treatment outcomes. High-risk patients with immune inflammation should be aggressively treated with immunotherapy to achieve better clinical outcomes. This strategy may also be used in other cancer patients to seek better treatment targets and carry out precise treatment of cancer.

Although the mechanism of metabolism in tumors has become an increasingly popular area of research in recent years, the potential regulation between tumor immunity and metabolism still needs more research. Different states of tumor-associated macrophages (TAMs) are able to adapt to the tumor microenvironment by altering metabolism. It can inhibit the differentiation of M2-type TAMs by inhibiting the metabolism pathway (Van den Bossche et al., 2016; Divakaruni et al., 2018). Cancer-related MDSCs, whose main energy supply mode is converted from glycolysis to FAO. Increased fatty acid uptake and higher expression of key enzymes, which in turn upregulated the FAO rate necessary to produce immunosuppressive ARG1 (Hossain et al., 2015). Cytokines that drive the expansion of MDSCs. Fatty acid transporters (FATPs), as long-chain fatty acid transporters, upregulate and control the inhibitory activity of MDSCs on tumors in tumor-derived MDSCs (Veglia et al., 2019). This evidence suggests that fatty acid metabolism targeting MDSCs and M2 TAMs may be an important means of enhancing the efficacy of cancer immunotherapy. In this study, antigen presentation signaling pathways differed significantly between different risk groups. One possible hypothesis is that differences in mediators in tumor tissues attract antigen-presenting cells (APCs) to the site of tumor cells (Friedmann Angeli, Krysko & Conrad, 2019). Our study shows that many immune-related biological processes and pathways are enriched among HCC patients with different metabolism risk groups. Therefore, it is reasonable to assume that metabolism in tumor tissues in HCC patients is closely related to tumor immunity.

Recent studies have found that amino acid metabolism (Bertero et al., 2019; Ericksen et al., 2019) and tricarboxylic acid metabolism (Nie et al., 2020) play an important role in tumor development and cancer drug resistance (Yoo & Han, 2022; Ma et al., 2018). As a member of glutamate-oxaloacetate aminotransferase, a recent study showed that GOT2 expression levels affect metabolism (Wang et al., 2018), suggesting that GOT2 may be an important regulator of cellular metabolism. In our study, the prognostic curve confirms that GOT2 plays an important role in predicting survival in patients with hepatocellular carcinoma. To gain a deeper understanding of the mechanism by which GOT2 affects tumor progression, wound healing assays are performed. These findings suggest that GOT2 has a significant activating effect on the migration ability of hepatocellular carcinoma.

In summary, our study defines a new prognostic model MRRS constructed from 14 MRGs. The model was shown to be independently correlated with OS in the derivation and validation cohorts, providing insights for predicting HCC prognosis. The study had several limitations. First, our prognostic model was constructed and validated with retrospective data from a public database. More forward-looking, real-world data are needed to validate its clinical utility. Secondly, since the model focuses on the evaluation ability of metabolism in HCC patients, a variety of indicators should be included in the comprehensive judgment of clinical decision-making. For example, it should be emphasized that the responsiveness of HCC patients in the high-risk score group to immunotherapy needs more research to elucidate.

Supplemental Information