Glucose metabolism 3-gene signature associates with poor prognosis in hepatocellular carcinoma as a hallmark of tumor microenvironment evolution and progression


 Background

At present, immunotherapy has become an established treatment for cancer. Recently, we have realized that tumor metabolism has a huge impact on the shaping of immune reactivity in the tumor microenvironment (TME). Exploration of the tumor metabolic and immunocellular response reveals metabolic vulnerability as a therapeutic window for intervention to enhance immunotherapy. The clinical significance of glycolysis and its role in tumor-immune evolution have not been fully explored. Understanding the relationships between glycolysis, TME evolution and disease progression is of great significance for optimizing immunotherapy of liver cancer.
Methods

A glycolysis-related biomarker signature was initially identified in transcriptomic data within The Cancer Genome Atlas (TCGA, n=424). Predicted overall survival (OS) for the gene signature in the GSE54236 dataset (n=161) and our HCC cohort (n=132) was subsequently independently verified. The relative number of immune cells in the microenvironment was calculated using CIBERSORT. CCK-8, IHC, Transwell and Seahorse were used to assess the effects of a 3-gene signature on the malignant phenotype of liver cancer.
Results

A 3-gene signature of glycolysis significantly associated with poor OS was identified. Based on a multi-omics research strategy, we found the 3-gene signature was closely related to the evolution of liver cancer TME and disease progression, and was associated with liver cancer mutation load and immune evolution. High-risk patients have inhibition and depletion of the immune TME and immune escape through antigen presentation inhibition. Remarkably, the 3-gene signature can predict the clinical outcomes of HCC patients with different clinical subtypes.
Conclusions

Our investigation of the 3-gene signature provided evidence for tumor metabolism-immune co-evolution along HCC progression. The 3-gene signature for stratification of liver cancer risk has robust predictive power, both at the RNA and protein levels.


Background
It is well known that hepatocellular carcinoma (HCC) is a heterogeneous tumor, mainly arising from chronic liver in ammation. Liver cancers have different pathogenesis, including chronic viral hepatitis infection, alcoholism and fatty liver [1,2]. HCC remains the third leading cause of cancer death worldwide due to its heterogeneity and limited options for targeted therapy [3]. Cancer immunotherapy is a new treatment to identify and attack cancer cells by activating or reconstructing the immune system of patients. At present, tumor immunotherapy, such as immune checkpoint inhibitors or adoptive cell therapy, has achieved clinical therapeutic effect in some patients with HCC, but most patients do not respond to tumor immunotherapy [4]. For example, only 20% of patients respond to PD-1 immune checkpoint blockade therapy [5,6]. CAR-T cell therapy is currently only used for the treatment of hematological malignancies, and the therapeutic effect for solid tumors is very limited [7]. The main reason is that tumor in ltrating lymphocytes (TILs) are continuously stimulated in the tumor microenvironment, resulting in T cell depletion. Even after PD-1 antibody treatment, depleted tumor in ltrating T cells are di cult to convert into an effective activated state [8]. The depletion of T cells, which is di cult to prevent, has become a bottleneck problem in the eld of tumor immunotherapy [9]. However, recent studies have reported that OXPHOS metabolic reprogramming related to the pyruvate-MPC metabolic pathway can activate depleted T cells in the microenvironment, revealing a new metabolic regulation pathway to combat T cell depletion [10]. This study also showed that depleted T cells did not respond to conventional immunotherapy drugs such as PD-1 antibody but could still respond to metabolic regulation drugs. Therefore, tumor immunotherapy based on metabolic regulation may be a new breakthrough, following immune checkpoint inhibitors. However, how tumor metabolism affects the function of immune cells is still unclear, which necessitates understanding of the relationship between tumor metabolism reprogramming and the immune microenvironment [11].
A hallmark of cancer is to avoid cell death while undergoing sustained growth and cell proliferation, which requires considerable energy costs [12]. Despite the existence of su cient oxygen, tumor cells still prefer to use glycolysis to obtain energy, with faster production of ATP to meet high energy demands, lower levels of reactive oxygen species that induce cellular damage and production of synthetic metabolic precursors required for tumors, known as the Warburg effect [13]. Glucose metabolism by glycolysis during rapid cell proliferation results in loss of nutrients for immune cells such as in ltrating T cells from the TME, acidi cation of the TME due to secretion of lactate, and enrichment of immunosuppressive metabolites [14,15]. [16,17]. Although it has long been known that the intrinsic metabolism of tumor cells can shape TME, the effect of this metabolic change on non-transformed cells (such as invasive immune cells) has not been fully understood [18]. It is not clear whether the changes of glucose metabolism and the interaction between immune cells play a role in immune editing and tumor escape, or whether glucose metabolism is related to tumor genome instability, thus driving tumor evolution. In addition, the clinical effect of glucose metabolism reprogramming in tumor evolution has not been fully explored. Given the multi-step nature of HCC carcinogenesis and disease progression, it is very important to study and understand the immune microenvironment and the evolution of tumor genome.
Our study aims to explore the relationship between glucose metabolism and liver cancer immune microenvironment, and to study its signi cance in tumor evolution and immune evolution, so as to provide a theoretical basis for tumor immunotherapy strategies based on metabolic regulation. Through machine learning and multi-factor cox regression analysis, we constructed a glucose metabolism model and developed a 3-gene signature based on SLC2A1 and its related genes, and veri ed its prognostic and diagnostic value in a variety of clinical subtypes. Using the multi-group research strategy, we found that the 3-gene signature was associated with instability of the tumor genome. With the increasing risk score, tumors remodelled from innate immunity to adaptive immunity, and the direction of immune microenvironment evolved from activation to inhibition / depletion. With this different immune evolution, tumors develop different immune escape strategies, such as antigen presentation inhibition. We also found 3-gene signature-related transcriptome immune networks and different molecular pathways involved in disease progression and immune status. Current studies have shown that with increasing risk scores, immune remodeling becomes another dimension of tumor-immune co-evolution, which can be used as a predictor of tumor progression.

Methods
The Cancer Genome Atlas dataset (TCGA, https://cancergenome.nih.gov/) The Cancer Genome Atlas (TCGA) is a joint project initiated in 2006 currently featuring a total of 36 cancer types, incorporating clinical data, mRNA expression, methylation, and other data for a variety of human cancers, including subtypes of tumors [19]. We downloaded and analyzed the LIHC RNA-Seq count data from the TCGA database (for 374 liver cancer and 50 para-cancer samples), from which the clinical data were extracted (Fig. 1a). The data were used for subsequent analysis after normalization.
As above, all data were normalized. Clinical Patient Tissue Between 2010 and 2014, 132 tumor samples were obtained from liver cancer patients at the Department of Urology, Second A liated Hospital of Wenzhou Medical University, each tumor tissue sample matched with two biopsies of normal tissue from the same patient. Patients ranged in age from 25 to 89 years, with a mean of 57.9 years. Each tissue sample was available both formalin-xed and para n-embedded (FFPE) for routine histopathological examination, and a separate piece of the same tissue was immediately frozen in liquid nitrogen and stored at -80°C for further analysis. Samples from each FFPE tissue block (an average of 20 tissue slices were assembled into one chip, including a total of 7 chips for 132 cases) were assembled into a tissue micro array (TMA; Fig. 1a). All patients were clinically and pathologically con rmed to have HCC. Tumor staging was based on the classi cations of the WHO and American Joint Committee on Cancer (AJCC)/International Union for Cancer Control (UICC). The study was approved by the Board of Directors and Ethics Committee of the Second A liated Hospital of Wenzhou Medical University. All patients participating in the study provided written informed consent.
Patient demographics and the clinicopathological characteristics of all patient cohorts in the study are detailed in Table 1.

Gene set enrichment analysis (GSEA)
To identify potentially enriched pathways, single-gene set enrichment analysis (sGSEA) was performed on interested genes [20]. The analysis used v4.10 of the GSEA tool. Genesets with P-values < 0.05, a normalized enrichment score (NES) > 1 and false discovery rate (FDR) < 0.25 were considered signi cant.
Identi cation of cell types by estimation of relative subsets of RNA transcripts (CIBERSORT) The CIBERSORT tool (https://cibersort.stanford.edu/index.php), generally used for describing cell composition based on the gene expression pro les of complex tissues, was used to enumerate hematopoietic subsets from mixtures of mRNA in the TCGA database [21]. CIBERSORT is superior to other methods in terms of noise, mixtures of unknown content, and closely related cell types.

RNA isolation and RT-qPCR
The total RNA of para-cancerous and liver cancerous tissues of liver cancer patients was extracted using Trizol reagent (Invitrogen). First-strand cDNA was synthesized using reverse transcription, serving as a template for qPCR. A 2 μL cDNA/20 μL reaction volume was performed using a CFX96 optics module connected to a C1000 thermocycler system (Bio-Rad) using 2× SYBR Green Supermix (Vazyme), in accordance with the manufacturer's protocol. The primer sequences are available in the Supplementary section (give Table number?). Thermal cycling parameters were as follows: 30 s at 95℃ for heat activation of DNA polymerase, 10 s at 95℃, and 10 s at 60℃ for 40 cycles of ampli cation. The number of cycles (Ct) required for threshold uorescence was recorded for each reaction. The Ct values for SLC2A1, NCAPH, and KIF11 were normalized to the housekeeping gene β-actin. RT was conducted on each sample in triplicate while each cDNA was evaluated in three separate pieces of tissue.

Western blot analysis
Expression of proteins was detected by Western blot assay as described previously [22].For cancerous tissues and matching normal adjacent tissue for each patient, 20μg total protein were extracted, separated on a 10% SDS-PAGE gel, then electro-transferred onto polyvinylidene uoride membranes (Merck Millipore, Darmstadt, Germany) in Tris-glycine buffer. A blocking buffer (5% skim milk in TBST) was used to block non-speci c reactions on the membrane at room temperature for 1.5 h, shaken well on a slow shaker. Each membrane was incubated overnight with primary antibodies against KIF11, NCAPH, SLC2A1, and β-actin at 4℃, respectively. After washing with TBST for 3 × 10 min, anti-rabbit or antimouse immunoglobulin G conjugated with horseradish peroxidase (HRP), were incubated at room temperature for 1h on a slowly shaking platform. In accordance with the manufacturer's instructions (BioSharp, China), immunoreactive proteins were observed using a Bio-Rad Chemidoc imaging system and ECL reagents. The optical density of protein bands was quanti ed using ImageJ software (National Institutes of Health).

Assessment of IHC staining
Two pathologists were blinded to the experimental details and clinicopathological characteristics and independently scored the TMA IHC. Three representative 200× microscopic images on elds chosen stochastically were acquired for each section using an inverted Leica optical microscope to quantify IHC staining, using Image Image-Pro Plus 6.0 (IPP) software (Media Cybernetics, Rockville, MD). Integrated optical density (IOD) values were calculated using the digital image analysis method described above [42], representing the integral of staining intensity for each pixel in the region multiplied by area, representing the total stained tissue in the region. Whole tissue microarray scanning was used to measure the levels of KIF11, NCAPH, and SLC2A1 proteins in each TMA.

Cell culture and transfection
The hepatocarcinoma cell line MHCC97H was purchased from the Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences. Cells were cultured with Dulbecco's modi ed Eagle's medium (DMEM) (Life Technologies) with 10% fetal bovine serum (FBS) (Life Technologies) and antibiotics (100 U/mL penicillin and 100 μg/mL streptomycin), growing at 37℃ in a humidi ed incubator with 5% CO 2 .
The cells (1 × 10 6 ) were inoculated on the 6 -well plate on the rst day. When cell density reached ~70%, 4 μg DNA or shRNA and 5 μl Lipofectamine 2000 (Invitrogen) were diluted with 125 μl OPTI-MEM medium respectively, and added to the 6-well plates in 1.1 ml complete medium and replaced with complete medium after 24 hours. DNA and shRNA were purchased from Public Protein/Plasmid Library. The experiment was conducted 24-72 hours after transfection.

CCK-8 assay
The transfected MHCC97H cells were seeded in 96-well plates at a density of 1.2 × 10 4 cells per well with six replicate wells. The cells were cultured overnight at 37℃ and the absorbance was detected at 450 nm using the Cell Counting Kit-8 (CCK-8, Dojindo Laboratories, Kumamoto, Japan) to re ect the cell proliferation ability.

Transwell migration and invasion assays
Transwell migration experiments were performed using 24-well Transwell chambers (Corning, USA ) and 8-μm pore size polycarbonate membrane. In the migration experiment, 100μL cell suspension containing 30,000 transfected MHCC97H cells were added to the upper chamber in cell culture medium that did not contain FBS, and 600μL of complete medium containing FBS was added to the lower chamber. After incubation at 37℃ for 18 hours, the chamber was xed and stained, imaged and counted under a 100fold magni cation microscope.

Results
Development and validation of a three-gene prognostic signature based on glucose metabolism for risk strati cation in liver cancer We identi ed a 3-gene risk strati cation signature using the following strategy (Fig. 1b). First, we used R packages ("limma" and "edgeR") to normalize the liver cancer RNA expression data from the TCGA, GSE14520, GSE25097, GSE57957, and GSE64041 datasets, with thresholds for p < 0.05 or 0.001 (based on data quality) [23]. From all cohorts combined, 512 differentially expressed genes (DEGs) were detected (Fig. S1a). The GeneCard database was used to search for the keyword 'glucose metabolism', obtaining 183 genes with scores greater than 30, after which the geneset for REACTOME_GLUCOSE_METABOLISM was then extracted from GSEA. By cross-comparison of DEGs, a total of 7 overlapping DEGs were identi ed. (Fig. S1b,c).
Then using the UALCAN and cBioPortal databases, the data was restricted to that with correlation scores > 30 (or correlation coe cients >0.4) to identify co-expressed genes for SLC2A1 [24,25]. By crosscomparison of DEGs in TGCA, a total of 317 DEGs positively correlated with SLC2A1 were identi ed (co-DEGs, Fig. S1d). Through the cytohubba plug-in for cytoscape, 50 genes were identi ed as hub nodes in the PPI network (Fig. S1e, f) [26]. We used the random forest algorithm to incorporate SLC2A1 and its 50 hub genes to further extract prognostic signature genes. By univariate and multivariate regression analysis, we identi ed that SLC2A1, KIF11, and NCAPH provided a distinct pro le (p < 0.05, Fig. 1c). We next showed the protein levels of these 3 risk factors within hepatocellular carcinoma and paraneoplastic tissue using IHC (Fig. 1d).
We then veri ed the effect of the three-gene prognostic signature in different datasets and different clinical subtypes (Fig. 1e). Speci c data reference supplementary materials are shown in Fig. S2. The results show that the model performs well in internal and external validation sets (LIHC and GSE54236). We found that except for clinical subgroups of women and G4 and T2 cohorts, all other subgroups of patients in the high-risk group displayed signi cantly shorter survival times than those in the other group (P < 0.05). In addition, we found that the prognostic accuracy for Asian Man G2 clinical subgroup, was quite remarkable, with 1 and 3 -year AUCs of 0.86 and 0.83, respectively.
To explore the expression levels of the 3 risk factors, we veri ed SLC2A1, NCAPH and KIF11 at genome, transcriptome, and proteome levels (Fig. 1f). Based on the GEO and TCGA databases, we validated their mRNA and DNA methylation levels with 617 liver samples, nding p-values that were statistically signi cant (P<0.05, Fig. S3). We also validated the diagnostic and prognostic value of three genes (Fig.  S4a,b). In addition, we collected 10 paired HCC and adjacent tissues to verify their protein expression levels by Western blot-based assays (Fig. S4c). We collected 132 cases of liver cancer tissues and veri ed them by IHC and qRT-PCR and the results were consistent. This was achieved by using simple, inexpensive quantitative PCR or IHC -based detection methods, which can be easily converted to clinical use.
Tumor DNA copy number variation events are associated with 3-gene signature.
Next, our goal was to explore 3-gene signature as a hallmark of tumor evolution along the path of tumor DNA copy number variation (CNVs). First, we scored patients for the 3-gene signature and differentiated them between high/low risk cohort by median value. The results demonstrated that CNVs events were highest in high risk group and lowest in the low risk group, indicating higher genome-wide instability and mutational burden in the high risk group (Fig. 2a,b). Deletion of commonly mutated tumor-associated genes, such as PTEN, NOTCH1 and CDKN2A, are associated with 3-gene signature (Fig. 2a,c). In addition, ampli cation events of VEGFA and MUC1 increased in high-risk groups, suggesting that tumors in highrisk group were predisposed to angiogenesis or metastasis (Fig. 2c). Interestingly, there were more CNV mutations in the CTNNB1 gene in the high risk group and no difference in APC mutations, but we found more APC deletion events in the high-risk group, consistent with the deletion f APC hindering degradation of CTNNB1 and causing accumulation of β-catenin in nucleus, potentially leading to abnormal differentiation of high-risk tumor cells (Fig. 2c) [27].
We further explored the relationship between CNVs or mutations and malignant behavior of tumors, such as cell cycle, tumor invasion and other signaling pathways. At the same time, the heat map veri cation of transcription group and malignant phenotype of tumor was performed (Fig. S5). We downloaded cell cycle, EMT, glycolysis and in ammation hallmarks genesets from GSEA database and counted their DNA mutation data. We observed more CNVs and mutation events in high risk group for each of these, implying that genomic aberrations are closely related to many behavioral hallmarks of tumor cells (Fig. 2d) [28].
In line with this, our 3-gene signature illuminated the dynamics of tumor progression process, such as ampli cation of proto-oncogenes and deletion of oncogenes associated with malignant behavior and progression of tumors.
3-gene signature knockdown effects proliferation, energy metabolism, migration and invasion in HCC cells Subsequently, we veri ed the correlation between the 3-gene signature and tumor malignant behavior at the cellular level. Firstly, MHCC97H cells were transfected with the sh-RNA or expression plasmids for these three genes. According to CCK-8 detection, we found that the silencing of SLC2A1, NCAPH and KIF11 inhibited the proliferation of MHCC97H cells, while the overexpression of SLC2A1, NCAPH and KIF11 stimulated it (Figure 3a), con rming that the 3-gene signature genes are involved in regulating the proliferation of tumor cells. In addition, Transwell analysis showed that overexpression of SLC2A1, NCAPH and KIF11 enhanced the migration and invasion of MHCC97H cells (Figure 3b), while the silencing of SLC2A1, NCAPH and KIF11 inhibited it (Figure 3b), indicating that the 3-gene signature promoted migration and invasion of HCC cells.
Based on the above studies, we know that these 3-gene signature genes are associated with DNA copy number variations of proto-oncogenes and tumor suppressor genes, and the ampli cation or deletion of these genes can lead to changes in tumor metabolic phenotype. For example, c-MYC directly regulates glucose and glutamine metabolism. Therefore, we designed the hippocampal energy metabolism experiment to explore the correlation between these three risk characteristics and the energy metabolism level of MHCC97H cells. The results showed that the silencing of SLC2A1, NCAPH or KIF11 inhibited the glycolysis and oxidative phosphorylation of MHCC97H cells, and inhibited ATP synthesis (Fig. 3c), thereby affecting the microenvironment of tumor cells.

3-gene signature is associated with an exhaustive and inhibitory tumor microenvironment
Firstly, we used CIBERSORT algorithm to evaluate the immune microenvironment composition in highversus low-risk groups [21]. We then mapped a comprehensive in ltration landscape of 22 immune cell types in HCC (Fig. 4Sd). To further explore the interactions of tumor immune cell in ltration in high and low risk groups, we mapped the correlation heat map of 22 immune cells in HCC, suggesting a different immunological microenvironment between these groups (Fig. 4a).
To further elucidate the immune cell in ltration differences in the high-versus low-risk groups, immune cells for which the differences were statistically signi cant were divided into two sets based on immune features for innate and adaptive immune subsets [29]. We found that the in ltration levels of dendritic cells, mast cells and NK cells were higher in the low-risk group, suggesting that the immune response of low-risk tumors was mainly innate immunity. On the contrary, we found that the in ltration levels of Tregs cells and T helper cells were higher in the high-risk group, suggesting that the immune response of highrisk tumors was mainly adaptive immunity. With the increase of 3-gene signature score, the immune microenvironment of liver cancer oversteps from innate immunity to adaptive immunity, suggesting the evolution of immune microenvironment. During tumorigenesis, mast cells are one of the rst immune cells to be recruited to the tumor. Similarly, mast cell in ltration is higher in the low risk group, suggesting that the low risk group is at an early stage of tumor evolution (P=0.00047, Fig. 4b) [30].
In addition, we examined the antigen presentation potential of different groups. Although there was no difference in dendritic cell in ltration, lower levels of monocyte cell in ltration was seen in the high-risk group (P=0.00036, Fig. 4b,). Researchers have found that monocytes can activate T cells e ciently and transfer tumor antigens directly to DCs through cell junctions, inducing intense adaptive immunity to achieve anti-tumor effects [31]. Although adaptive immune function is activated in high-risk groups, low levels of monocytes suggest that antigen presentation is inhibited and eventually immune escape occurs. Interestingly, there was a negative correlation between the in ltration of monocytes and the in ltration of M2 macrophages and Tregs in the immune microenvironment of the high risk group, suggesting that the high risk group had an inhibitory immune microenvironment [32]. The higher in ltration of Thelper cells and Treg cells also con rms the depletion and suppression of immune microenvironment in high risk groups [33]. Our IHC results also suggested a higher level of macrophage in ltration in the high risk group, consistent with our data above (Fig. 4c).
Furthermore, we explored the gene expression levels of CD274 (PD-L1), CD80 (B7-1), CTLA4, HAVCR2, LAG3, PDCD1, TIGIT (VSIG9), and VSI (VISTA) in the low-and high-risk groups (Fig. 4d). The expression of these tumor immune checkpoints was higher in high-risk LIHC groups than in the low-risk group (P < 0.0001), suggesting that suppression of the immune microenvironment contributes to decreased antigen presentation and immunosuppression [34]. The use of immune checkpoint inhibitors in the high-risk group is a clinical tool that may be rewarding for clinical management of HCC.

Immune function status networks and transcriptomic pathways of the 3-gene signature
To further emphasize the effect of the 3-gene signature on overall immune activation in tumors, compared with the low risk group, we found 2405 differentially expressed genes, (log(fold change) ≥ 1, p < 0.05). Then we used the ClueGO plug-in in Cytoscape to enrich the immune group. The results showed that up-regulated genes and down-regulated genes were enriched in different immune function networks, suggesting that different risk groups had different immune response patterns. The high-risk group mainly activated the negative regulation of immune system process, adaptive immune response and lymphocyte activation pathways (Fig. 5a). Conversely, the low-risk group mainly activated the positive regulation of immune system process and innate immune response pathways (Fig. 5b). This result further illustrated that high risk group tumours had inhibitory and depleted immune microenvironments.
Next, we performed biological process (GO term) enrichment analysis of these DEGs. We found that organic acid catabolism, carboxylic acid catabolism, fatty acid metabolism process and other metabolic processes were enriched in the low-risk group tumors, suggesting that multiple acid metabolic pathways may be involved in maintaining TME immune status in HCC (Fig. 5c). Conversely, high-risk group tumors were mainly enriched in cell cycle or division processes, which are known to be associated with tumor cell proliferation and diseases development [12]. Volcano maps showed some genes involved in fatty acid metabolism and nuclear division and it can be seen that most genes involved in fatty acid metabolism are down-regulated in high-risk group tumors, suggesting that acid metabolism is related to evolution of immune microenvironment (Fig. 5d).
To further clarify the association of the 3-gene signature with tumor evolution and development, we divided the tumors into 4 grades according to the degree of histologic differentiation. The lower the degree of tumor differentiation, the higher the risk score, which suggests that the signature score is related to tumor heterogeneity and evolution (Fig. 5e). In addition, we collected low, medium, and high-risk tumor tissues and observed them under the microscope. The results showed that the higher the risk score is, the worse the level of organizational differentiation is, consistent with our data above (Fig. 5f).

3-gene signature can robustly identify poor disease prognosis and survival in patients with liver cancer
Although genome and transcriptome 3-gene signature have been related to bad prognosis of many kinds of cancers, the clinical correlation of metabolic regulators is still unknown [35,36]. By exploring the relationship between 3-gene signature score and multiple clinical features, we found it was closely related to clinical staging and pathological staging, which indicates the occurrence and development of tumors and poor prognosis (Fig. 6a). Notably, patients with nonalcoholic fatty liver (MASH/MAFLD) have higher risk scores (Fig. 6a). Researchers found that NASH increased abnormal T cells in the liver, and the activation of immunotherapy promoted the progression of liver cancer [37]. The immune microenvironment of tumors with NASH underwent more complex metabolic reprogramming, which made immune checkpoint inhibitors ineffective for such patients. The 3-gene signature was also associated with other parameters such as liver ablation, person neoplasm cancer status, platelet function, chronic HBV infection and BMI index clinical parameters (Fig. 6a).
Besides, patients in the high-risk group displayed signi cantly shorter survival times than those in the low-risk group (P < 0.0001, Fig. 6b). Next, using univariate Cox analysis, we revealed that there was a signi cant association between risk score and OS, and following multivariate Cox regression analysis after adjusting for age, grade, T, M, or N stage, the 3-gene signature score was demonstrated to be an independent prognostic factor (P = 0.01062, Fig. 6c). These results strongly support the ability of the 3gene signature to accurately stratify poor prognosis in patients with liver cancer using transcriptomic data.
More importantly, in order to be more convenient, inexpensive, rapid and amenable to large-scale use in clinical diagnosis and treatment monitoring, we collected pathological tissue sections of patients with liver cancer, and veri ed the 3-gene signature status at the proteome level by IHC. The results showed that patients in the high-risk group displayed signi cantly lower OS than those in the low-risk group, showing that the diagnostic accuracy for poor prognosis in patients by IHC was also quite remarkable (P = 0.001, HR = 3.57, Fig. 6d). In addition, ROC curves also showed strong predictive power of these 3 genetic features at the proteomic level for the prognosis of HCC patients (Fig. 6e). We found that the 3-gene signature displayed strong predictive capability at both transcriptomic and protein levels, so the constructed 3-gene signature prognosis model was robust and reliable, and easily applicable to clinical practice.
In addition, we tested the predictive ability of the 3-gene signature in both internal and external validation sets. Indeed, we con rmed that patients with high risk scores had poorer OS than those with low risk scores, which further proves the reliability and robustness of the 3-gene signature (Fig. 6f).
The predictive capability of the three-gene diagnostic model in clinical tissue was also validated with areas under the ROC curve, respectively, of 0.957 (TMA dataset) and 0.885 (TCGA dataset) ( (Fig. 6g)).
The results suggest that the 3-gene signature has considerable potential as diagnostic marker of liver cancer.
Finally, in order to establish an accurate method to predict the survival probability of liver cancer patients, we developed an easy-to-use and clinically adaptable, risk nomogram for predicting the survival probability in LIHC patients (Fig. 6h). For example, a patient with a higher risk score, G3 and 60 years old produced a total score of 160 points (100 points for high-risk score, G3 phase of 30 points, age of 30 points), and the predicted 1-year, 3-year and 5-year OS rates are 78.0%, 54.0% and 38.0%, respectively.

Discussion
The relationship between tumor glycolysis processes and microenvironment immune evolution has never been explored before. The study of glycolysis genes on the prognosis of patients with liver cancer has also focused on transcriptomics [38,39]. Based on the multi-group studies at the genome, transcriptome and proteome levels, we found a prognostic model consisting of 3-gene signatures and glycolysis-related genes, the relationship between the model and the evolution of immune microenvironment, and the impact on the prognosis and clinical relevance of patients with different subtypes of liver cancer. We observed that with the increase of 3-gene signature risk score, the tumor immune microenvironment showed an evolutionary trend from innate immunity to adaptive immunity, and from activated immune microenvironment to inhibitory immune microenvironment. In this case, tumors used different mechanisms to achieve immune evasion. The 3-gene signature also identi ed pathways related to tumor progression and immune evolution, such as multiple acid metabolism plays an important role in this process, which was ignored by most researchers in the past [40]. In addition, our data suggest that the 3gene signature scores are associated with poor prognosis and multiple clinical parameters, and were validated in cohorts of different clinical subtypes, which will be more easily translated into easy-to-use, more scienti c and personalized risk assessment and clinical analysis. This also emphasizes that the 3gene signature is a marker of tumor evolution and a key indicator of disease progression. In summary, our study found extensive crosstalk between tumor glucose metabolism and immune microenvironment.
With the reprogramming of glucose metabolism, TME evolution and tumor mutation co-evolution supported the tumor metabolism-immune co-evolution model [41,42].
Tumor glycolysis can affect lactic acid, glucose, and pH value in the immune microenvironment, and these metabolic pressures can affect the formation of tumor immune microenvironment [43][44][45]. Previous studies have focused on the effect of immune cell glycometabolism reprogramming on its own function, but the effect of tumor glycolysis on the immune microenvironment is largely unknown. Watson et al.
found that lactic acid plays an important role in regulating the immune response of tumor microenvironment by supporting the metabolism of tumor in ltrating regulatory T cells [46]. Our study focuses on the effect of tumor glycolysis on the function of the immune microenvironment, and the metabolic pressure caused by tumor glycolysis represents an evolutionary event related to clinical results. In conclusion, we identi ed a 3-gene signatures that could serve as a biomarkers of evolutionary events in tumor-immune progression. It is more important to recognize that the interaction between metabolism and immunity is heterogeneous, dynamic and bidirectional. For example, TME glucose restriction caused by tumor glycolysis can signi cantly affect T cell function. Blocking immune checkpoints can inhibit tumor glycolysis and restore glucose in TME [47]. From our research data, as the level of 3-gene signature glycolysis increases, the TME of HCC changes from innate immunity to adaptive immunity, from activated immune microenvironment to inhibitory immune microenvironment, which is the manifestation of immune landscape remodeling and evolution. On the other hand, with the increase of 3-gene signature glycolysis level, more mutations are accumulated in tumors, and the degree of differentiation is getting worse and worse, which supports the concept of tumor-immune coevolution proposed in previous studies [42].
Our research results also emphasize that tumors with high 3-gene signature glycolysis level have higher expression levels of all common immune checkpoints and higher tumor mutation load, suggesting that immunotherapy strategies for high-risk patients can bring greater returns [48]. In addition, the resistance of tumor high metabolism to metabolic pathway-dependent immunotherapy should also be considered.
For example, CTLA4 and PD1 signals can interfere with glycolysis by reducing the activation of AKT [49].
Melanoma patients with with high expression of glycolysis-related genes show worse progression-free survival after inhibiting PD1 treatment, but the treatment effect is greatly improved after limiting tumor glycolysis [50]. In addition, recent studies on the function of Tregs cells have shown that high glucose conditions can damage the function and stability of Tregs cells. Tregs can adapt to TME by using the glycolysis product lactic acid [46]. Consistently in our ndings, the tumor microenvironment function of Tregs cells at high 3-gene signature glycolysis level is not affected and the in ltration level is higher, which indicates that the metabolic pressure produced by tumor glycolysis has different effects on different immune cells. Therefore, studying how metabolism affects the response to immunotherapy and further understanding of tumor immune-metabolism dynamics not only reveals many insights into immune regulation from the perspective of metabolism, but also reveals its important potential therapeutic applications.
In addition, the 3-gene signature provides a model for predicting tumor progression and prognosis of patients with liver cancer. Although there are many studies on prognosis models of liver cancer at present, and there are many prognosis models with excellent predictive ability, they are rarely used in clinical practice [51,52]. They are mainly affected by three main obstacles. Most studies are based on transcriptome gene expression pro les. Due to its high cost, long turnover time and the need for bioinformatics expertise, it is di cult to apply to real clinical scenarios. Secondly, some researchers are prone to dozens of prognostic models constructed by genes. Although the effect is good, the use is too cumbersome, which limits its clinical application [53]. Thirdly, because RNA is easily degraded and unstable, the prognostic model based on transcriptome analysis does not verify its robustness at the protein level, which limits its clinical application [54,55]. Our prognosis model is composed of fewer genes and has excellent effect, and its stability is veri ed at the protein level, which provides a promising method for predicting the prognosis of liver cancer. By only simple immunohistochemical methods, levels of the three protein can be directly detected in patients and, are simple and stable. The 3-gene signature can be directly applied to real clinical scenarios, and its prediction effect is still robust and reliable.
In summary, our study interpreted the complexity of metabolic-immune interactions in tumors and provides evidence that the 3-gene signature is a marker of tumor-immune coevolution in HCC progression, which provides a theoretical basis for tumor immunotherapy based on metabolic regulation. In addition, the reliability of the 3-gene signature was established and validated, which provides an easy-to-use tool for identifying high-risk HCC patients and predicting their prognosis (risk score determined by immunohistochemistry). Availability of data and materials: All data generated and described in this article are available from the corresponding web servers and are freely available to any scientist wishing to use them for noncommercial purposes, without breaching participant con dentiality. Further information is available from the corresponding author on reasonable request.

Abbreviations
Ethics approval and consent to participate Not applicable.     LIHCSupplementaryMaterials01222022.pdf