Identification and Validation of an 6-Metabolism-Related Gene Signature and Its Correlation With Immune Checkpoint in Hepatocellular Carcinoma

Hepatocellular carcinoma (HCC) is a common malignant tumor with relatively high malignancy and rapid disease progression. Metabolism-related genes (MRGs) are involved in the pathogenesis of HCC. This study explored potential key MRGs and their effect on T-cell immune function in the tumor immune microenvironment to provide new insight for the treatment of HCC. Of 456 differentially expressed MRGs identified from TCGA database, 21 were screened by MCODE and cytoHubba algorithms. From the key module, GAD1, SPP1, WFS1, GOT2, EHHADH, and APOA1 were selected for validation. The six MRGs were closely correlated with survival outcomes and clinicopathological characteristics in HCC. Receiver operating characteristics analysis and Kaplan-Meier plots showed that these genes had good prognostic value for HCC. Gene set enrichment analysis of the six MRGs indicated that they were associated with HCC development. TIMER and GEPIA databases revealed that WFS1 was significantly positively correlated and EHHADH was negatively correlated with tumor immune cell infiltration and immune checkpoint expression. Finally, quantificational real-time polymerase chain reaction (qRT-PCR) confirmed the expression of WFS1 and EHHADH mRNA in our own patients’ cohort samples and four HCC cell lines. Collectively, the present study identified six potential MRG biomarkers associated with the prognosis and tumor immune infiltration of HCC, thus providing new insight into the pathogenesis and treatment of HCC.


INTRODUCTION
Hepatocellular carcinoma (HCC) is one of the most common aggressive malignancies worldwide and the third leading cause of cancer-related death (1). In 2021, estimated approximately 42,230 new cases were diagnosed and 30,230 deaths occurred from HCC in the United States (2). In China, HCC is the most commonly diagnosed cancer and the leading cause of cancer death in men <60 years of age, with 466,100 new cases and 422,100 deaths in 2015 (3). Currently, the treatment of HCC consists mainly of hepatectomy, transplantation, or ablation. Despite advances in the clinical treatment and diagnosis of HCC, recurrence rates remain high and survival rates remain poor (1). Therefore, innovative approaches to HCC treatment are urgently needed to reduce disease recurrence and death.
Advances in genomic and proteomic technologies have improved our understanding of the pathogenesis of malignant tumors. Metabolism refers to the orderly chemical reactions that occur in an organism to sustain life. Dysregulation of cellular metabolism is considered as an important factor associated with carcinogenesis (4). During tumorigenesis, cancer cells produce lactate via glycolysis (Warburg effect), which has a profound effect on the tumor microenvironment (TME) and the proliferation of cancer cells (5). Recent studies showed that alterations in metabolic pathways, such as lipid metabolism, were strongly associated with the development and progression of HCC, and metabolic processes played an important role in the occurrence and progression of HCC (6). The identification of metabolismrelated genes (MRGs) and the development of bioinformatics technologies may lead to the identification of reliable biomarkers for the early diagnosis of HCC, as well as for predicting the progression and prognosis of this disease. In this study, we constructed a six-MRG prognostic signature based on The Cancer Genome Atlas (TCGA) database. We experimentally validated the mRNA expression levels of WFS1 and EHHADH in our own patients' cohort and HCC cell lines, and further explored their expression patterns, potential biological functions, and immune infiltration levels in HCC. Recent studies showed that immune checkpoint inhibitor (ICI) therapies, such as those targeting programmed cell death 1 (PDCD1), cytotoxic T lymphocyte-associated protein 4 (CTLA4), and lymphocyteactivation gene 3 (LAG3) showed significant benefits in terms of survival compared with conventional therapies (7,8). Identification of potential prognostic markers associated with therapeutic benefit may allow the individualization of immunotherapy for patients with HCC. Our findings suggested that WFS1 and EHHADH played essential roles in tumor growth by triggering immune checkpoints in HCC, and are thus promising prognostic biomarkers for patients receiving immunotherapy.

Data Collection and Pre-Processing
HCC expression profiles included 373 HCC tissue samples and 50 peri-tumor samples and corresponding clinical information were obtained from TCGA (https://cancergenome.nih.gov/) and used as an analysis dataset. HCC data from International Cancer Genome Consortium (ICGC) included 212 HCC tissue samples and 177 peri-tumor samples) and four raw HCC-associated gene expression datasets, including GSE14250 (247 HCC tissue samples and 239 peri-tumor samples), GSE36376 (240 HCC tissue samples and 193 peri-tumor samples), GSE54236 (81 HCC tissue samples and 80 peri-tumor samples), and GSE107170 (146 HCC tissue samples and 104 peri-tumor samples), which were downloaded from NCBI-GEO database (https://www.ncbi.nlm. nih.gov/gds/), were used as an external validation dataset. The MRGs were acquired from GeneCards (https://www.genecards. org/), which provides comprehensive information on all annotated and predicted human genes. The term "metabolism" was used as a keyword for the search, and genes with relevance scores >5 were considered as MRGs.

Differentially Expressed MRG Screening and Functional Analysis
Differentially expressed genes (DEGs) in TCGA cohort were screened using the edge R package. Threshold values of log 2 FC>1 and adjusted P <0.05 were considered statistically significant. Intersecting DEGs and MRGs were selected as key genes for HCC metabolism for further analysis. To investigate the biofunctional characteristics and potential pathways of the identified MRGs, Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed by an online annotation, visualization, and integrated discovery database. (DAVID; https://david.ncifcrf. gov/) (9). P <0.05 was deemed statistically significant.

Protein-Protein Interaction (PPI) Network Construction and Hub MRG Identification
The PPI network was constructed based on selected MRGs by the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database and those with co-expression coefficients >0.4 were extracted (9). Thereafter, Cytoscape software was used to visualize and analyze molecular interaction networks (10). The Molecular Complex Detection plugin (MCODE) was then used to screen significant gene clusters, and a score >5 was used as the selection criterion. Hub genes were selected from the intersection of the top 120 genes and calculated with 10 topological analysis methods using the Cytoscape plugin cytoHubba. Finally, the results of the two algorithms were intersected to obtain the final hub MRGs.

UALCAN Database and Kaplan-Meier Analysis
UALCAN is a comprehensive, user-friendly, and interactive web resource based on the TCGA dataset for analyzing cancer OMICS data (11). This database was used to analyze the differential expression of 21 MRGs in tumor/peri-tumor samples, and key genes were selected from the 21 MRGs that were correlated with multiple clinical disease characteristics. The Kaplan-Meier plotter could assess the effect of 54 k genes (mRNA, miRNA, protein) on survival in 21 cancer types including breast (n = 7,830), ovarian (n = 2,190), lung (n = 3,452), and gastric (n = 1,440) cancer (12). Sources for the database include TCGA etc. The prognostic value of the mRNA expression of 19 MRGs in HCC was assessed using the Kaplan-Meier plotter. A log P-value <0.05 was considered statistically significant.

The cBioPortal for Cancer Genomics Analysis
The cBioPortal for Cancer Genomics (http://www.cbioportal. org/) is a freely accessible platform for interactive exploration of multi-cancer genomic datasets such as the TCGA database (13). The frequency of genetic alterations of six hub MRGs in patients with HCC and their association with survival outcomes were explored using cBioPortal.

Construction of a Prognostic MRG Signature and Gene Set Enrichment Analysis (GSEA)
To evaluate the prognostic value of the six selected MRGs, we used the HCC clinical dataset from the TCGA database and implemented the six-gene signature via multivariate Cox regression analysis using the Sangerbox (http://sangerbox.com/ Tool) tool. The integrated prognostic ROC, KM curve plotting tool for the relationship between risk score and gene expression developed by ggplot2 package based on R language on Sangerbox website performed multivariate Cox regression analysis on the six metabolism-related gene prognostic risk models. Risk classification thresholds were generally defaulted to the median. GSEA is a computational method for determining whether a set of a priori defined genes shows statistically significant consistent differences between two biological states. The HCC gene expression matrix from the TCGA database was downloaded and subjected to GSEA to predict potential hallmarks. After executing the permutation test with 1,000 permutations, gene sets with a P-value < 0.05 were considered statistically significant.

Immune Infiltration Analysis
Tumor Immune Estimation Resource (TIMER; https://cistrome. shinyapps.io/timer/) is a comprehensive resource for the systematic analysis of immune infiltrates across diverse cancer types (14). In this study, the expression of six MRGs involved in tumor purity and immune infiltration abundance (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and DCs) in HCC was explored. The connection between immune infiltration and gene copy number variation (CNV) was also analyzed. The immune scores, stromal scores and tumor purity were calculated based on the ESTIMATE algorithm (15). The immunophenscore (IPS) was calculated on a 0-10 scale based on the expression of the representative genes or gene sets of the immunophenogram. Higher IPS was associated with stronger immunogenicity, indicating a better response to ICI. The IPSs of HCC patients were obtained from The Cancer Immunome Atlas (https://tcia. at/home).

Relationship Between Hub MRGs and
Immune Checkpoints of HCC GEPIA (http://gepia.cancer-pku.cn/) is a web tool that provides customizable functions such as tumor/normal differential expression analysis, profiling according to cancer type or pathological stage, patient survival analysis, similar gene detection, correlation analysis, and dimensionality reduction analysis (16). The correlation of screened MRGs with immune checkpoints in HCC was assessed using the GEPIA database. Statistical significance was determined using R >0.1 and P <0.05 as filter criteria.

Clinical Patient Samples Collection and Gene Expression Validation
A total of 10 patients with HCC diagnosed in the General Surgery Department of Tianjin Medical University General Hospital from October 2015 to July 2017 were collected. Tumor tissues and paired non-tumor control tissues have a distance > 5 cm. This study was approved by the ethics committee of Tianjin Medical University General Hospital. All subjects gave informed consent prior to participation in the investigation. Furthermore, four HCC cell lines were selected to further validate the expression levels of WFS1 and EHHADH. Total RNA was extracted using the TRIzol reagent (Invitrogen, CA, USA), and the concentration was measured using a spectrophotometer (NANODROP; Thermo Fisher Scientific, USA). Then, total RNA was reverse-transcribed to synthesize cDNA, which was used to perform qRT-PCR for two MRGs. bactin was used as internal control, and fold change (2− △△CT ) was used to express the relative gene expression levels. The primer sequences are listed in Table 1.

Statistical Analysis
Statistical analysis was performed using GraphPad Prism (version 8.0; San Diego, CA, USA). The Student's t-test was used to compare differences between the tumor and peri-tumor groups. The data were considered statistically significant at P <0.05.

Screening of Differentially Expressed MRGs and Enrichment Analysis
The mRNA expression matrix of TCGA HCC patients was used to identify DEGs. Comparison of HCC and adjacent tissues in TCGA identified 4,786 DEGs, as shown in the volcano plot in Figure 1A. Analysis of the GeneCards website identified 1,857 MRGs with  Figure 1C). Among them, 198 genes were downregulated and 258 genes were upregulated ( Figure 1B); the heat map shows the distribution of the 456 MRGs ( Figure 1D). The function of these DE-MRGs was examined by GO and KEGG analyses. The top 10 GO terms and KEGG pathways are listed in Figure 1E. For biological processes (BPs), DE-MRGs were mainly enriched in oxidation-reduction processes, lipid metabolic processes, metabolic processes, and xenobiotic metabolic processes ( Figure 1E). Cellular component (CC) analysis showed that DEGs were mostly enriched in the extracellular exosome, cytosol, and high-density lipoprotein particle ( Figure 1F). Regarding molecular function (MF), genes were primarily enriched in heme binding, monooxygenase activity, and arachidonic acid epoxygenase activity ( Figure 1G). KEGG pathway enrichment analysis showed that MRGs were significantly associated with material synthesis and material metabolism signaling pathways, such as biosynthesis of amino acids, drug metabolism -cytochrome P450, and carbon metabolism ( Figure 1H).

Hub MRGs Screened Through the PPI Network
The STRING online database and Cytoscape software were used to construct the PPI network, which contained 423 nodes and 2,193 edges ( Figure 2A). Module analysis was performed using the MCODE plugin in Cytoscape software, and three modules were identified, comprising 77 hub genes ( Figures 2B-D). Next, the cytoHubba plugin was used to identify the central gene, and  26 hub genes were selected by intersecting the results from the ten algorithms of cytoHubba. The results of the two algorithms were intersected to obtain the final 21 MRGs ( Figure 2E), which are shown in Figure 2F.

Gene Expression and the Prognostic Value of Identifying Hub MRGs in HCC
The transcriptional levels of the 21 screened MRGs were compared between HCC and peri-tumor tissues using UALCAN. The results showed statistically significant differences in the expression of 19 MRGs except GCG and IL6 ( Figure 3A). Kaplan-Meier analysis showed that the 19 MRGs were related to clinical prognostic outcomes ( Figures 3B, C). The 19 MRGs were subjected to repeat GO functional analysis and KEGG pathway analysis. BP term enrichment mainly involved triglyceride metabolic processes, lipoprotein metabolic processes, and cellular oxidant detoxification ( Figure 3D). The CC terms included extracellular space, extracellular region, and extracellular exosome ( Figure 3E). MF included antioxidant activity, pyridoxal phosphate binding, and receptor binding ( Figure 3F). KEGG pathway analysis showed that the MRGs were associated with carbon metabolism, pathways in cancer, PPAR signaling pathway, and PI3K-Akt signaling pathway ( Figure 3G). Most enriched pathways were involved in tumorigenesis and metabolic pathways.

Validation of the Expression Levels of Six Hub MRGs
Among the 19 MRGs screened, we selected six (GAD1, SPP1, WFS1, GOT2, EHHADH, and APOA1) that have not been reported extensively in HCC, and their transcription levels and association with clinicopathological features were analyzed. GAD1, SPP1, and WFS1 were upregulated, whereas GOT2, EHHADH, and APOA1 were downregulated in HCC samples compared with peritumor controls, and this was consistent with the results in the ICGC, GSE14520, GSE54236, GSE107170, and GSE36376, validation datasets ( Figures 4A-E). To verify the six MRGs at the translational level, protein expression data were obtained from the Human Protein Atlas (HPA) database. The protein expression levels of SPP1, WFS1, GOT2, EHHADH, and APOA1 showed a similar pattern to that of mRNA transcript levels ( Figure 5). Furthermore, the six MRGs were significantly differentially expressed in HCC samples from patients with different tumor stage, grade, age, and sex according to the UALCAN analysis ( Figure 6).

Genetic Alterations of Six Screened MRGs in HCC
The correlation between the six hub MRGs was calculated using the TCGA-LIHC dataset, and Pearson's correlation analysis indicated significant correlations among them ( Figure 7A).

Evaluation and Validation of the Prognostic Signatures of Six MRGs
To evaluate the classification performance of the six-MRG signature, we used the TCGA cohort as a training dataset. The risk scores for each sample were computed using the Sangerbox website, and the samples were divided into high and low risk groups according to the median risk score. Kaplan-Meier curves showed that patients in the high-risk group had markedly worse OS than those in the low-risk group ( Figure 8A). The area under the curve (AUC) of the six gene signatures showed 1-, 3-, and 5year OS rates of 0.72, 0.68, and 0.7, respectively, suggesting that the prognostic model had favorable sensitivity and specificity ( Figure 8C). The risk score, survival time, survival status, and expression heatmap of the six-MRG signature are presented in Figure 8B. Moreover, the GSE14520 dataset were used for the validation of six metabolism-related gene signatures. Consistent with the results of the TCGA database, the Kaplan-Meier curve showed that patients in the high-risk group exhibited markedly worse OS than the low-risk group ( Figure 8D). The AUC of the six gene signatures for 1-, 3-and 5-year OS rates were 0.7, 0.71 and 0.62 ( Figure 8F). The risk score, survival time, survival status, and gene expression heatmap of six prognostic MRGs signature were shown in Figure 8E.

Exploration of Significant Pathways for Six MRGs by GSEA
To further investigate the potential functions of the six MRGs in HCC, we performed hallmark analysis by GSEA using the TCGA-LIHC dataset. GAD1 and SPP1 were enriched in the IL6-JAK-STAT3 signaling pathway, estrogen-response-late pathway, inflammatory-response pathway, TNF-a-signaling-via-NF-kB pathway, and glycolysis pathway ( Figures 9A, B). WFS1 was enriched in the unfolded protein-response and glycolysis pathways ( Figure 9C). GOT2, EHHADH, and APOA1 were all enriched for bile-acid-metabolism, fatty-acid-metabolism, xenobioticmetabolism, peroxisome, and adipogenesis ( Figures 9D-F). The pathways were markedly enriched in high risk samples and most of them were involved in metabolism and carcinogenesis.

Association of the Six MRGs With Intratumor Immune Cell Infiltration
To detect potential correlation between Estimate score/immune/ stromal and six MRGs, the ESTIMATE algorithm was used to conduct the calculations. The results showed that SPP1, GOT2, and EHHADH were significantly associated with immune score or stromal score ( Figure 9G).Then, the role of the six MRGs in the tumor immune microenvironment in HCC was examined by assessing the correlation between the expression of the six MRGs   (GOT2, EHHADH, and APOA1) was significantly negatively correlated with the infiltration of six immune cell subtypes ( Figure 10A). Also, the correlation and estimated statistical significance between six MRGs expression and other immune cell signature infiltration were displayed in Supplementary Table 1.
The results showed that six MRGs were significantly correlated with the majority of immune cells such as NK cells, Tregs, and T cell follicular cell. The relationship between the CNV of six hub MRGs and tumor immune infiltrating cells is shown in Figure 10B.

The Expression of WFS1 and EHHADH Is Correlated With Immune Checkpoints
Immune checkpoints are responsible for tumor immune escape.
Considering the potential function of MRGs in HCC, we examined the relation of the six MRGs with the major immune checkpoints associated with HCC (PDCD1, CD276, CTLA-4, LAG3, TIM3, OX40, and TIGIT) (17,18). The results demonstrated that WFS1 and EHHADH were significantly associated with immune checkpoint biomarkers. As shown in Table 2, relative WFS1 expression was significantly positively correlated with PDCD1, LAG3, TIM3, OX40, TIGIT, and CD276 levels. The tumor samples were then ranked according to the median level of WFS1 expression. The upper median was labeled WFS1-high and the lower median was labeled WFS1-low. The expression of PDCD1, LAG3, TIM3, OX40, TIGIT, and CD276 was significantly higher in the WFS1-high sample group than in the WFS1-low expression group ( Table 3). EHHADH expression was significantly negatively correlated with T cell exhaustion markers in HCC, which showed a significant decrease in the EHHADH-high expression group (Tables 2, 3). We then sought to determine whether the characteristics of the WFS1 and EHHADH genes could provide therapeutic value. We used IPS to predict the effectiveness of ICI in patients which were divided into high-and low-risk groups based on WFS1 and EHHADH gene expression levels. The low-risk group exhibited higher IPS (Supplementary Figures 1A, B). The WFS1 low-risk group exhibited higher IPS-PD1/PDL1/PDL2 scores than the highrisk group, representing higher immunogenicity and thus predicting a better response to ICI (Supplementary Figure  1A). Taken together, these results suggested that WFS1 and EHHADH played a critical role in tumor immune escape, which was involved in the carcinogenesis of HCC.

Experimental Validation of WFS1 and EHHADH mRNA Expression
The results of bioinformatics analysis were confirmed by qRT-PCR in 10 paired HCC with peri-tumor samples and 4 HCC cell lines, the LO2 cell line was used as the control. The results showed that WFS1 was significantly upregulated (Figures 11A, C), whereas EHHADH ( Figures 11B, D) was significantly downregulated in 10 HCC samples and 4 HCC cell lines, which was consistent with the bioinformatics results.

DISCUSSION
The age of incidence of HCC has decreased gradually in association with higher malignancy, rapid disease progression, and lack of obvious clinical symptoms in the early stages of the disease. Therefore, elucidating the molecular mechanisms of hepatocarcinogenesis may provide important clues for finding effective therapeutic targets or for identifying promising prognostic biomarkers. In recent years, there has been an increased interest in the study of metabolic alterations and dysregulation of the TME (5). Increasing evidence suggests that MRGs play a pivotal role in the development and progression of HCC (18). Considering the importance of the metabolic environment in tumor development, it is vital to identify metabolism-related HCC biomarkers. In this study, we detected 456 DE-MRGs in TCGA dataset, including 198 downregulated and 258 upregulated MRGs. The 21 most significant hub MRGs were identified using the cytoHubba and MCODE algorithms with Cytoscape software. Of these, 19 MRGs were most significantly associated with transcriptional levels and clinical prognostic outcomes in HCC. Among these 19 hub MRGs, we selected six that have not been reported extensively in HCC, namely GAD1, SPP1, WFS1, GOT2, EHHADH, and APOA1, to investigate their diagnostic and prognostic value. GAD1 catalyzes the synthesis of g-aminobutyric acid inhibitory neurotransmitters (19), and it is overexpressed in various neoplastic tissues, although its expression in HCC has not been reported. SPP1 is a secreted glycophosphoprotein that is associated with multiple biological functions, including cell adhesion, migration, and invasion (20,21); however, its role in HCC remains unclear. The WFS1 gene, encoding wolframin (WFS1), causes endoplasmic reticulum stress and is linked to a rare autosomal-recessive disorder known as Wolfram syndrome (22). However, its role in the development of HCC remains unclear. Similarly, the prognostic role of GOT2, EHHADH, and APOA1 in HCC remain so far obscure. Compared with peri-tumor samples, GAD1, SPP1, and WFS1 were significantly upregulated in tumor tissues, and GOT2, EHHADH, and APOA1 were significantly downregulated in tumor tissues. The reliability of the expression of the six hub MRGs was further assessed in five datasets, including ICGC, GSE14520, GSE36376, GSE54236, and GSE107170, and the results were consistent with the cohort expression levels. UALCAN analysis indicated that the transcription levels of GAD1, SPP1, and WFS1 were markedly higher in HCC tissues than in peri-tumor control tissues in the subgroup analyses. Increased gene expression was associated with worse survival, suggesting that these three MRGs function as oncogenes in HCC.
In contrast, GOT2, EHHADH, and APOA1 were associated with better survival outcomes and were significantly downregulated in HCC tissues, as well as positively correlated with tumor stage, grade, age, and sex, suggesting that these three MRGs were tumor suppressor genes. Furthermore the gene alteration frequencies of the above six MRGs in HCC were explored. Most of the six MRGs showed genetic alterations including gene amplifications, deep deletions, and missense mutations. Moreover, cBioPortal analysis showed that alterations in the six MRGs correlated with worse OS (p = 1.914E -3 ) and DFS (p = 4.053E -4 ). These results indicated that genetic alterations of the six MRGs may serve as valuable biomarkers for HCC diagnosis and treatment. In this study, we assessed the classification performance of the six-MRG signature model using TCGA database. The ROC curves showed that the six-MRG signature performed well in predicting the prognosis of HCC patients, especially for OS within 1, 3, and 5 years, implying that all six MRGs can be used as biomarkers to predict prognosis with good sensitivity and specificity. The prognostic value of this six MRG signatures was further successfully validated in the GSE14520 dataset. Moreover, the GSEA results suggest that many HALLMARK pathways related to tumorigenesis and metabolism, such as IL6-JAK-STAT3-signaling, TNF-asignaling-via-NF-kB, glycolysis pathway, bile-acid-metabolism, and fatty-acid-metabolism were enriched in the expression groups of the six MRGs, suggesting their contribution to HCC progression through metabolic regulation. The IL6-JAK-STAT3 pathway consists of extracellular IL-6 ligands that activate the IL-6 receptor, which phosphorylates JAK, which in turn phosphorylates STAT3. Phosphorylated STAT3 dimerizes and translocates into the nucleus to induce the expression of genes with various tumor-promoting properties (23). Therefore, exploring biomarkers related to the IL6-JAK-STAT3 pathway may be potentially valuable for the detection and treatment of cancer. In addition, NF-kB is a family of inducible transcription factors that play multiple evolutionarily conserved roles in the immune system (24).TNF-a cytokines induce rapid transcription of genes regulating cell survival, proliferation and differentiation mainly through activation of the NF-kB pathway (25). Previous studies have reported that TNF-a activation of NF-kB is one of the most characteristic tumorigenic cytokines in hepatocellular carcinogenesis (26). Deep understanding of TNFa-signaling-via-NF-kB pathway would throw lights on the interaction of MRGs and tumor progression. Similarly, glycolysis pathway, bile-acid-metabolism, and fatty-acidmetabolism also involved in tumor progression through affecting cellular metabolic processes (27)(28)(29).
Tumor-infiltrating immune cells affect the efficacy of radiotherapy, chemotherapy, or immunotherapy, and thus the prognosis of cancer patients (30)(31)(32). We showed that GAD1, SPP1, and WFS1 were markedly positively correlated with diverse immune cells, including B cells, CD8+T cells, CD4+T cells, macrophages, neutrophils, and dendritic cells in HCC. In contrast, GOT2, EHHADH, and APOA1 were significantly negatively correlated with tumor infiltrating immune cells. These findings indicated that tumor immune infiltration may partially account for the oncogenic role of metabolism-related genes in HCC. Additionally, the efficacy of immunotherapy requires not only sufficient immune cell infiltration in the TME, but also depends on the adequate expression of immune checkpoints (33). Compared with conventional therapies, ICI therapies such as those targeting programmed death ligand 1 (PD-L1), PDCD1, and CTLA4, have shown significant survival benefits (17). The activity of the immune system is primarily regulated by immune cells called T cells. In the TME, the response of T cells is regulated by a group of cell surface molecules called immune checkpoints. Cancer cells can evade the immune system, mainly by overexpressing suppressive immune checkpoints that cause Tcell exhaustion, thereby evading T-cell attack (34). In the HCC tumor microenvironment, PDCD1 is expressed in CD8+ T cells, and when activated, PDCD1 inhibits T cell migration, proliferation and secretion of cytotoxic mediators, thereby blocking the "cancer immune cycle" (35). In human liver cancer tissues, PD-L1 expression was mainly expressed in Kupffer cells, and PD-L1+Kupffer cells interacted with PD-1+CD8+ T cells, resulting in hepatocellular carcinoma effector T cell dysfunction, suggesting that PD-L1/PD-1 immune checkpoints could be used as targets for the treatment of HCC (35,36). LAG-3 and T-cell immunoglobulin and mucin structural domain molecule-3 (Tim-3), are also upregulated on specific CD8+ T cells in various cancer types and are also involved in the progression of hepatocellular carcinoma. Recently, LAG3 expression was previously found to be significantly higher in CD8+ tumor-infiltrating helper T cells and CD8+ cytotoxic T cells than in tumor-free liver tissue and blood from patients with hepatocellular carcinoma (37). T-cell immunoglobulin-3 (Tim-3) is a type I transmembrane protein that is expressed by most tumor-infiltrating Tregs in the tumor setting and appears to represent a specific tissue-resident Treg subpopulation with enhanced suppressive activity (38). It has also been reported that Tim-3 is strongly expressed on CD4+ and CD8+ T cells obtained from hepatocellular carcinoma lesions, and these data suggest that the functional relationship between Tim-3 and different T cells may affect the prognosis of HCC (37). TIGIT was identified as a co-receptor expressed mainly by activated and regulatory T cells (Treg) and NK cells (39). It has been found that TIGIT can contribute to the progression of hepatocellular carcinoma by modulating the immunosuppressive effect of cytokines produced by CD155 on CD8 T cells (40). Similarly, high expression of OX40, CD276, CTLA4 was significantly associated with poor prognosis of hepatocellular carcinoma (41)(42)(43).
Cancer cells are highly dependent on stress adaptation and altered metabolic states for tumor cell proliferation (44). Metabolic changes in tumor cells can also affect the composition and function of tumor-infiltrating lymphocytes  Bold values indicate that P < 0.05 was considered statistically significant.
residing in the same microenvironment, such as T-cell exhaustion, which in turn causes tumor progression (44). Many tumors are capable of evading the immune system, primarily by overexpressing suppressive ligands that inhibit Tcell attack. Consequently, fewer and impaired T cells are found in patients with HCC, which contributes to the progression of this cancer. We therefore assessed the relationship between the six MRGs and immune checkpoints. The results showed that high expression of WFS1 was closely associated with PDCD1, LAG3, TIM3, OX40, TIGIT, and CD276 in HCC, and low expression of EHHADH was negatively correlated with PDCD1, LAG3, TIM3, OX40, TIGIT, CD276, and CTLA-4. This indicated that WFS1 or EHHADH could promote inhibitory receptor induction and exhaustion and that targeting both may increase the efficacy of immunotherapy for HCC.

CONCLUSION
This study identified six MRGs that were closely related to the tumorigenesis and prognosis in HCC by multi-cohort dataset and comprehensive bioinformatics analyses. The current findings suggested that WFS1 and EHHADH exerted their oncogenic effects by activating immune checkpoints, indicating that their expression levels may serve as biomarkers for predicting the efficacy of ICI in HCC patients, which should be further explored in the future. The elucidation of the underlying mechanism and the validation of our observations in clinical trials will help develop more effective treatments for HCC patients.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. The analyzed data in this study can be found here: GEO (http://www.ncbi.nlm.nih.gov/ geo/) and TCGA (https://portal.gdc.cancer.gov/). Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
HR and NZ performed bioinformatics analysis and paper draft. HG, WW, and WL revised the images. XL and SL performed the