A Novel Lactate Metabolism-Related Gene Signature for Predicting Clinical Outcome and Tumor Microenvironment in Hepatocellular Carcinoma

Abstract

Hepatocellular carcinoma (HCC) is the main subtype of primary liver cancer with high malignancy and poor prognosis. Metabolic reprogramming is a hallmark of cancer and has great importance on the tumor microenvironment (TME). As an abundant metabolite, lactate plays a crucial role in cancer progression and the immunosuppressive TME. Nonetheless, the potential roles of lactate in HCC remain unclear. In this study, we downloaded transcriptomic data of HCC patients with corresponding clinical information from the TCGA and ICGC portals. The TCGA-HCC dataset used as the training cohort, while the ICGC-LIRI-JP dataset was served as an external validation cohort. Cox regression analysis and the LASSO regression model were combined to construct the lactate metabolism-related gene signature (LMRGS). Then, we assessed the clinical significance of LMRGS in HCC. Besides, enriched molecular functions, tumor mutation burden (TMB), infiltrating immune cells, and immune checkpoint were comprehensively analyzed in different LMRGS subgroups. In total, 66 differentially expressed lactate metabolism-related genes (LMRGs) were screened. The functions of LMRGs were mainly enriched in mitochondrial activity and metabolic processes. The LMRGS comprised of six key LMRGs (FKTN, PDSS1, PET117, PUS1, RARS1, and RNASEH1) had significant clinical value for independently predicting the prognosis of HCC patients. The overall survival and median survival of patients in the LMRGS-high group were significantly shorter than in the LMRGS-low group. In addition, there were differences in TMB between the two LMRGS subgroups. The probability of genetic mutations was higher in the LMRGS-high group. Most importantly, the LMRGS reflected the TME characteristics. In the LMRGS-high group, the immune microenvironment presented a suppressed state, accompanied by more inhibitory immune cell infiltration, including follicular helper T cells and regulatory T cells. Additionally, the expression of inhibitory checkpoint molecules was much higher in the LMRGS-high group. Our study suggested that the LMRGS was a robust biomarker to predict the clinical outcomes and evaluate the TME of patients with HCC.

Introduction

Hepatocellular carcinoma (HCC) is the most common histological type of primary liver cancer, the third leading cause of cancer death worldwide (Sung et al., 2021). As a highly malignant tumor, the 5-year survival rate of HCC is less than 18% (Villanueva, 2019). Treatment options for HCC include hepatic resection, liver transplantation, image-guided ablation, transarterial therapies, chemotherapy, and molecularly targeted therapy (Llovet et al., 2021). Clinically, patients with HCC are often treated by a combination of several modalities. However, the therapeutic outcomes of advanced HCC remain unsatisfactory. Even after successful tumor eradication, the recurrence rate of HCC is remarkably high. Recently, immunotherapy has been shown to improve the clinical efficacy of advanced HCC. Unlike the mechanism of action of conventional therapy, immunotherapy is based on activating the patient’s own immune system to fight against tumors (Ringelhan et al., 2018). Cancer metabolism plays an essential role in affecting the anti-tumor immune response through modulating the interaction between tumor cells and the tumor microenvironment (TME) (Bader et al., 2020). Therefore, it is vital to identify a metabolism-related signature to assess the TME and improve the treatment efficacy of immunotherapy.

Metabolic alterations of tumor cells not only favor cell proliferation but also have profound influences on anti-tumor immunity through the release of metabolites, especially lactate (Xia et al., 2021). Unlike normal cells, tumor cells metabolize glucose to produce lactate even under adequate oxygen conditions. The accumulation of lactate provides an acidic microenvironment that benefits tumor growth and progression. Besides, lactate produced by aerobic glycolysis can be secreted into the extracellular environment as a signaling molecule to regulate intercellular interactions (Liao et al., 2021). In gastric cancer, lactate derived from tumor cells mediates the up-regulation of BDNF expression in cancer-associated fibroblasts by activating the NF-κB pathway, eventually resulting in acquired resistance (Jin et al., 2021). Alterations in lactate metabolism have been shown to be associated with cell invasion, migration, angiogenesis, drug resistance, and immune escape. High levels of lactate in the TME promote differentiation of tumor-associated macrophages to the M2 subtype, while activated macrophages facilitate tumor invasion through the CCL17/CCR4/mTORC1 signaling axis (Zhang et al., 2021). Lactate-induced PD-L1 up-regulation on neutrophils impairs T cell cytotoxicity in HCC (Deng et al., 2021). In addition, tumor cell-derived lactate induces the expression of GPR81 in dendritic cells via paracrine mode to inhibit the antigen presentation function of immune cells (Brown et al., 2020). Moreover, lactate has an important role in epigenetic regulation. Some studies have demonstrated that histone lysine lactylation takes part in modulating gene transcription (Izzo and Wellen, 2019; Yu et al., 2021). Given the vital role of lactate in oncogenesis and the immunosuppressive TME, targeting its metabolism promises to become an effective means for cancer treatment.

In this study, we screened the key lactate metabolism-related genes (LMRGs) and constructed a prognostic signature to predict the survival outcome. Next, we comprehensively analyze the tumor mutation burden (TMB) features in different subgroups. Then, the association between the TME and lactate metabolism-related gene signature (LMRGS) was explored using the R software package. We focused on the infiltrating immune cells in the TME and characterized the differential immune microenvironment in LMRGS subgroups. The results indicated that the LMRGS had a high value for evaluating the prognosis and reflecting the TME in HCC.

Materials and Methods

Data Acquisition

RNA transcriptome sequencing data, somatic mutation profile, and corresponding clinical information of HCC were obtained from the TCGA data portal (https://portal.gdc.cancer.gov/). In this study, the TCGA-HCC cohort was served as the training set. To verify the training set results, we downloaded an independent dataset of HCC from the ICGC website (https://dcc.icgc.org/releases/current/Projects/LIRI-JP). Therefore, the ICGC-LIRI-JP cohort was used as a validation set. The detailed clinical information of HCC patients from two cohorts was summarized in Table 1.

Differentially Expressed LMRGs and Transcription Factors

The 289 LMRGs were retrieved from the Molecular Signatures database (Liberzon et al., 2015). Transcription factors associated with cancer were downloaded from the Cistrome (Zheng et al., 2019). To identify the LMRGs and transcription factors involved in the progression of HCC, we carried out differential expression analysis between 50 normal tissues and 374 tumor tissues in the TCGA-HCC cohort. Genes with |log2 fold change (FC) | > 1 and false discovery rate (FDR) < 0.05 were defined as differentially expressed. For further understanding the biological function and pathway of differentially expressed LMRGs and transcription factors, we used the “clusterprofiler” package in R (version 4.1.0) to carry out the GO and KEGG enrichment analyses.

Construction and Assessment of LMRGS

The differentially expressed LMRGs were subjected to univariate Cox regression analysis to determine the LMRGs with prognostic value. To avoid overfitting, we further performed the LASSO Cox regression (iteration = 1000) using the “glmnet” package (Friedman et al., 2010; Liu et al., 2021a; Liu et al., 2021b). After screening by LASSO regression, the selected LMRGs were applied to establish the LMRGS through the multivariate Cox regression analysis. The LMRGS score was calculated as the following formula: LMRGS score = expression level of gene₁ × coefficient of gene₁ + expression level of gene₂ × coefficient of gene₂ + … + expression level of gene_n × coefficient of gene_n. We classified HCC patients into two subgroups according to the median LMRGS score, including the LMRGS-high and the LMRGS-low groups. Principal component analysis (PCA) was used to evaluate the classification accuracy of the signature. For assessing the prognostic value of the LMRGS, we conducted the Kaplan–Meier (KM) survival analysis to compare the overall survival (OS) and median survival time between the two LMRGS groups. The time-dependent ROC curve was performed by the “timeROC” package in R. We also applied the Cox proportional hazards regression model to identify the LMRGS as an independent predictor for OS. To explore the influence of the LMRGS on HCC progression, we clarify the association between the LMRGS and clinicopathologic factors, such as TNM stage, pathological grade, fibrosis, vascular invasion, and virus infection.

Establishing a Nomogram

To predict the one-, three-, 5-year survival rate of HCC patients, we constructed a nomogram based on the LMRGS and significant clinicopathologic parameters (Iasonos et al., 2008). The calibration curve was used to estimate the consistency between predicted survival and actual survival. The time-dependent ROC curve was applied to evaluate the specificity and sensitivity of the model.

Calculation of TMB

For calculating the TMB of each HCC tumor sample, we selected the somatic mutation data processed by the VarScan platform in the TCGA-HCC cohort. Then, we compared the difference of TMB between the LMRGS-high and the LMRGS-low groups. Visualization of somatic mutations in the two LMRGS groups was performed by the R package “maftools”. Moreover, we explored the impact of the LMRGS score combined with the TMB on the survival of HCC.

Comprehensive Analysis of TME in Different LMRGS Subgroups

The TME is mainly composed of stromal cells and immune cells (Gysler and Drapkin, 2021). Firstly, we used the ESTIMATE algorithm to calculate the stromal score of all samples (Yoshihara et al., 2013). ESTIMATE is a prevalent R package, which is widely utilized in the cancer-related studies (Liu et al., 2021c; 2021d; 2021e). Then, the single sample gene set enrichment analysis (ssGSEA) was performed to derive the immune enrichment score based on the 29 immune gene sets (Bindea et al., 2013). To identify the immune infiltration features of HCC samples, we imported their gene expression profiles to the CIBERSORTx website with 1000 permutations (https://cibersortx.stanford.edu/). According to the obtained results, we compared the relative fractions of 22 tumor-infiltrating immune cells in the two LMRGS subgroups. Moreover, correlation analysis was carried out to clarify the relationship between the immune cell and the LMRGS score. Immune checkpoints expression and immune function have crucial influences on the treatment responses of immunotherapy. For further investigating the effect of the LMRGS score on immunotherapy, comparisons between the two LMRGS subgroups were analyzed to evaluate the differences of immune checkpoints and immune function.

Gene Set Enrichment Analysis

The HCC samples were stratified into high- and low-LMRGS score groups as described above. To determine the primary signaling pathways and hallmark gene sets involved in the signature, we uploaded sample grouping and gene expression files into the GSEA software (version 4.1.0) to conduct enrichment analysis.

Statistical Analysis

All data analysis and visualization were completed by R software. If the data did not follow a normal distribution and the variance was uninformed, the differences between groups were compared by the Wilcoxon rank-sum test or Kruskal–Wallis test. The Cox regression model was used to perform univariate and multivariate analyses. The log-rank test was performed to evaluate the survival difference. Correlation analyses of LMRGS score and immune infiltration cells were conducted by Spearman’s rank correlation test. In this study, p-value < 0.05 was considered statistically significant as indicated.

Results

Identification of LMRGs

Through the differential gene screening analysis, we obtained 66 differentially expressed LMRGs, including three down-regulated and 63 up-regulated genes. The heat map displayed the expression of LMRGs in HCC samples and normal samples (Figure 1A). The differential expression of down-regulated and up-regulated LMRGs was represented in the volcano plot (Figure 1B). The 66 differentially expressed LMRGs were further analyzed by functional enrichment analysis. The primary biological processes (BP) of LMRGs were involved in mitochondrial genome maintenance, mitochondrial respiratory chain complex assembly, electron transport chain, and metabolic process. For cellular components (CC), the LMRGs primarily existed in the mitochondrial inner membrane, respiratory chain complex, and mitochondrial respirasome. The molecular functions (MF) of LMRGs were mainly enriched in electron transfer activity, NADH dehydrogenase activity, and oxidoreductase activity (Figure 1C). Signaling pathway analysis indicated that the differentially expressed LMRGs were related to thermogenesis, diabetic cardiomyopathy, oxidative phosphorylation, non-alcoholic fatty liver disease, and reactive oxygen species (Figure 1D). The above results showed that the LMRGs were mainly associated with metabolic processes and oxidation responses.

FIGURE 1

Development of the LMRGS

To identify the LMRGs correlated with OS, we performed the univariate Cox regression analysis. A total of 29 LMRGs were related to prognosis (Figure 2A). After selection by LASSO regression, only 10 LMRGs were subjected to multivariate Cox regression analysis to construct the LMRGS (Figures 2B,C). Based on the coefficient and the expression of six crucial genes involved in the LMRGS, we calculated the LMRGS score (Figure 2D). The LMRGS score of every HCC patient was obtained as follows: LMRGS score = FKTN expression × 0.2496 + PDSS1 expression × 0.0881 + PET117 expression × 0.0648 + PUS1 expression × 0.0567 + RARS1 expression × 0.0362 + RNASEH1 expression × 0.0928.

FIGURE 2

Prognostic Significance of the LMRGS

Taking the median LMRGS score as cut-off, we divided the HCC patients into two subgroups: LMRGS-high and LMRGS-low groups (Figure 3A). The heat map showed the differential expression of six crucial genes in the two LMRGS subgroups (Figure 3B). The LMRGS score and survival status of every HCC patient were displayed in Figure 3C. KM analysis indicated that patients with the high LMRGS score had shorter OS and median survival than patients with the low LMRGS score (Figure 3D). According to the different clinical characteristics, subgroup survival analysis also confirmed this result (Supplementary Figure S1). As shown in Figure 3E, the area under curve (AUC) value of 1 year, 3 years, and 5 years for ROC analysis was 0.768, 0.691, and 0.666, respectively, in the TCGA cohort. Moreover, the univariate and multivariate regression analyses demonstrated that the LMRGS score was an independent risk factor for OS (HR = 3.576, 95%CI = 2.105–6.074, p = 2.44E-06) (Table 2). The correlation of the LMRGS score and clinicopathological factors was clarified in the TCGA cohort. The results suggested that the LMRGS score was closely associated with pathological grade, clinical stage, vascular invasion, and virus infection (Supplementary Figure S2). The above results indicated that the LMRGS score played a vital role in HCC progression.

FIGURE 3

A Nomogram for Predicting Survival

To accurately predict the probability of OS, we established a nomogram that integrated the LMRGS score and other clinicopathological features, including age, gender, and TNM stage (Figure 4A). We could estimate the survival rate of 1 year, 3 years, and 5 years based on the total points. The calibration curve demonstrated that the prediction value was highly consistent with the actual value (Figure 4B). The time-dependent ROC curve also indicated that this nomogram had high accuracy for predicting survival (Figure 4C).

FIGURE 4

Regulation Network of Transcription Factors

There exist close interactions between the LMRGs and transcription factors. For exploring the relationship, we carried out the co-expression analysis. As displayed in Figure 5A, we identified 52 differential expressed transcription factors co-expressed with six significant LMRGs. The main functions of co-expressed transcription factors were chromatin remodeling and histone modification (Figure 5B). KEGG analysis revelated that these transcription factors mainly participated in the cell cycle, cellular senescence, and Hippo signaling pathway (Figure 5C).

FIGURE 5

Association With TMB

In the TCGA training cohort, we calculated the TMB of each HCC patient. We found that the TMB was higher in the LMRGS-high group (Figures 6A,B). Then, mutant situations of different LMRGS subgroups were visualized by the waterfall plots (Figure 6C). For the entire dataset, the top 10 mutated genes in HCC were TP53, CTNNB1, TTN, MUC16, ALB, PCLO, APOB, RYR2, MUC4, and FLG. Missense mutations were the most common somatic mutational types. The mutation frequency of samples was higher in the LMRGS-high group. Moreover, patients with high LMRGS scores had a higher mutation probability of crucial genes, especially TP53. Subsequently, we performed KM analysis to evaluate the influence of the LMRGS score combined with the TMB on survival. The result showed that the survival time of the high-TMB group was shorter than the low-TMB group (Figure 6D). More importantly, patients with a low LMRGS score and low TMB had a significantly longer OS than patients with a high LMRGS score and high TMB (Figure 6E). In the ICGC validation cohort, we also analyzed the mutation profiles of all samples. There existed no TMB difference among the two LMRGS subgroups (Supplementary Figures S3A,B). However, the mutation frequencies of prevalently mutated genes in HCC were higher in the LMRGS-high group (Supplementary Figure S3C). Survival analysis results of the LMRGS score combined with the TMB were consistent with the training cohort (Supplementary Figures 3D,E).

FIGURE 6

TME Characteristics in Different LMRGS Subgroups

Stromal cells and immune cells in the TME have profound impacts on tumor progression, treatment efficacy, and clinical outcomes. The heatmap shown in Figure 7A and Supplementary Figure S4A displayed the stromal score and immune activity of all samples. We found that the abundance of stromal cells was relatively higher in the LMRGS-low group. In addition, the LMRGS-low group had higher immune scores than the LMRGS-high group (Figure 7B and Supplementary Figure S4B). As shown in Figure 7C and Supplementary Figure S4C, there were differences in immune function between the LMRGS-high and LMRGS-low groups. The activity of cytolysis and IFN response was higher in the LMRGS-low group. In the LMRGS-high group, there was a higher expression of MHC class Ⅰ.

FIGURE 7

To comprehensively analyze the immune microenvironment, we used the CIBERSORTx to calculate the infiltration degree of 22 immune cells. The immune landscape of TCGA-HCC samples was shown in Figure 8A. By comparing the immune cell profiles, we found that follicular helper T (Tfh) cells, regulatory T cells (Tregs), and M0 macrophages were significantly increased in the LMRGS-high group. On the contrary, resting NK cells, monocytes, resting mast cells, and activated mast cells infiltrated more in the LMRGS-low group (Figure 8B). Apart from immune cells, we further explored the correlation of immune molecular and the LMRGS score. In our results, the LMRGS score was positively associated with the expression of immune checkpoints, including PD-1, CTLA4, LAG3, TIM3, and TIGIT (Figure 9).

FIGURE 8

FIGURE 9

GSEA of the LMRGS

To explore the molecular mechanisms involved in the LMRGS, GSEA was used to analyze the TCGA cohort. Enrichment results of hallmark revealed that DNA repair, E2F targets, G2M checkpoint, glycolysis, mitotic spindle, mTOR signaling, MYC targets, and unfolded protein response were activated by the LMRGS-high group (Figure 10A). Besides, the LMRGS also participated in regulating the transcription factors, DNA repair, cell cycle, and metabolism-related signaling pathways (Figure 10B).

FIGURE 10

Discussion

Despite some advances in diagnosis and treatment, HCC is still cancer with high morbidity and mortality (Forner et al., 2018). As inflammation-driven cancer, there is an intricate interplay between the TME and HCC development (Ringelhan et al., 2018). Increasing evidence indicates that metabolic changes of tumors can sculpt their microenvironment, and then the remodeled TME confer a growth advantage to tumor cells (Dimri et al., 2020; Li et al., 2021). Aerobic glycolysis is a vital hallmark of tumor metabolic reprogramming. Glucose is not completely oxidized but metabolized to produce lactate, even in the presence of oxygen (Palsson-McDermott and O'Neill, 2013). Recently, some studies have reported that there is lactate accumulation in tumors (Yu et al., 2021). Lactate is now considered an essential energy substance for tumor metabolism and plays an indispensable role in restructuring the TME (Certo et al., 2021). Hence, we constructed a novel LMRGS based on LMRGs in this study. The results suggested that the LMRGS was an independent prognostic factor for OS. In addition, the LMRGS proved to have substantial value for predicting the TME in HCC.

The LMRGS was composed of six crucial genes, including FKTN, PDSS1, PET117, PUS1, RARS1, and RNASEH1. FKTN participates in protein glycosylation modification (Kanagawa et al., 2016). A study of gastric cancer indicated that higher FKTN expression is associated with tumor progression, which may be due to the protein encoded by FKTN promoting the interaction between tumor cells and the extracellular matrix (Oo et al., 2016). PDSS1 is the critical enzyme in CoQ10 biosynthesis, mediating metabolism and mitochondrial function. The mutation of PDSS1 has an impact on ATP production and oxidative stress (Mollet et al., 2007). As for PET117, mainly distributed in the mitochondrial matrix, it is related to oxidative phosphorylation via influencing the biogenesis of cytochrome c oxidase (Vidoni et al., 2017). PUS1 involves in the structural modification of mRNA and is correlated with mitochondrial disorders (Carlile et al., 2019). RARS1 fusion with MAD1L1 has been reported to stimulate the FUBP1/c-Myc signaling pathway, inducing tumorigenesis in nasopharyngeal carcinoma (Zhong et al., 2018). RNASEH1 plays a vital role in maintaining the stability of mitochondrial DNA under oxidative stress (Renaudin et al., 2021). However, the role of six essential genes remains unclear in HCC. To clarify the regulation mechanism of these crucial genes, we performed co-expression analysis between transcription factors and six genes. A total of 52 co-expressed transcription factors were identified, and their functions were mainly reflected in chromatin remodeling and histone modification.

Genomic alterations are the main intrinsic drivers of tumor heterogeneity (Müller et al., 2020). To further understand the molecular features, we compared the gene mutations in different LMRGS groups. As suggested by the results, missense mutation was the most common type of mutations. In the LMRGS-high group, the TP53 gene had the highest mutation rate, while CTNNB1 and TTN were the most frequently mutated genes in the LMRGS-low group. TP53 is not only playing a central role in response to genotoxic stress but also in regulating metabolic homeostasis (Levine, 2020). Increasing evidence reveals the critical functions of TP53 in cellular metabolism (Wang et al., 2018; Kim et al., 2019). The dysfunction of p53 protein encoded by TP53 affects the tumor initiation and progression by mediating the metabolism of tumor cells (Lonetto et al., 2019). The poor prognosis of the LMRGS-high group could be due to TP53 hypermutation. A study reported that increased lactate better meets the metabolic needs of tumor cells and thus favors cell proliferation in p53 mutated tumor cells (Boidot et al., 2012). CTNNB1 and TTN also have links with the malignant transformation of liver cells (Jhunjhunwala et al., 2014). However, patients in the LMRGS-low group had lower probabilities of genetic mutations than those in the LMRGS-high group. Based on the gene mutations of the whole genome, the TMB of every patient was calculated. We found that patients with high TMB and high LMRGS scores had the worst clinical outcomes, which might be because of the genome instability caused by the high TMB (Ferguson et al., 2015).

Complex TME influences tumor progression and response to treatment. There were great differences in the TME between the two LMRGS subgroups, especially in the tumor immune microenvironment. The two groups showed different immune function statuses, including cytolytic activity, MHC class I expression, and IFN response. Cytolytic activity of immune cells reflects the ability to kill tumor cells. Transcriptome hypomethylation of CD8⁺ T cells activates cytolytic activity and effector function, which in turn enhances anti-tumor responses (Loo Yau et al., 2021). In HCC, patients with a high cytolytic activity score have favorable TME and more robust immunogenicity, resulting in better prognoses (Takahashi et al., 2020). Increased expression of MHC class I with high T cell infiltration benefits the prognosis of patients with liver metastases from colon cancer (Turcotte et al., 2014). In our analysis, MHC class I expression was higher in patients with high LMRGS scores. Consequently, the impact of MHC class I expression on the prognosis of HCC patients needs to be further clarified. Besides, the activation of IFN response is an essential link to anti-tumor immunity (Takahashi et al., 2021). As it could be seen, patients in the LMRGS-low group had better anti-tumor immune activity.

Tumor-infiltrating immune cells are one of the most important components in the TME, which can be affected by the lactate level (Certo et al., 2021). Low glucose and high lactate accumulation in the TME have immunosuppressive effects. Under lactate-rich conditions, reducing NAD⁺ to NADH by lactate dehydrogenase (LDH) leads to blocked production of GAPDH and PDGH, which in turn impairs effector T cell proliferation dependent on post-GAPDH glycolytic intermediates (Quinn et al., 2020). Tregs have inhibitory effects on immune response and antigen activation, facilitating cancer progression. Increased aerobic glycolytic activity creates a lactate-enrich microenvironment that favors Tregs survival and contributes immunosuppressive functions (Wang et al., 2017). Moreover, elevated lactate levels in the TME can supply potential nutrition to Tregs, which is due to lactate reversal to pyruvate and NADH in the presence of LDH (Lochner et al., 2015). A study suggested that inhibiting glycolysis and promoting oxidative phosphorylation recover the differentiation of Tfh cells and reduce inflammatory damage (Dong et al., 2019). Another interesting study found that high lactate accumulation decreases the PH of the microenvironment, then promotes NK cell apoptosis and inhibits its natural killer function (Harmon et al., 2019). B cells are of great significance in humoral immune responses through antibody production. Altered intra- and extracellular metabolic signaling can affect the immune regulatory function of B cells (Rosser and Mauri, 2021). Monocytes and mast cells play a vital role in regulating immune responses, and they can alter the TME toward anti-tumor immunity when fully triggered (Guilliams et al., 2018; Dudeck et al., 2019). In addition, macrophages have two central polarization states, including M1 and M2. Different TME leads M0 macrophages polarization to different states, resulting in very opposed effects. M1 macrophages polarization contributes to the immunity against the tumor, while M2 macrophages promote cancer progression and treatment resistance (Chen et al., 2021). Lactate derived from tumors leads to M2 macrophages polarization via activating the mTORC2 and ERK signaling pathways (Zhang et al., 2021). The results of our study were consistent with these conclusions. The infiltration levels of B cells, NK cells, monocytes, and mast cells were higher in the LMRGS-low group. Conversely, Tfh cells, Tregs, and M0 macrophages were more abundant in the LMRGS-high group. The results indicated that the immune cells of patients in the LMRGS-high group were affected by lactate metabolism, so the TME was more inclined to an immunosuppressive state.

Apart from the accumulation of immune cells that negatively regulate immune activity, the immunosuppressive TME is also associated with the up-regulated expression of inhibitory immune checkpoints (Sangro et al., 2021). We further explored the differences in the expression of inhibitory molecules between the LMRGS subgroups. In the LMRGS-high group, inhibitory immune checkpoint expressions were significantly higher, including PD-1, CTLA4, LAG3, TIM3, TIGIT, and CD96. In addition, the LMRGS score was positively correlated with PD-1, CTLA4, TIM3, and TIGIT. Recently, immunotherapy targeting inhibitory immune checkpoints has shown promising efficacy in treating advanced HCC (Yau et al., 2020). The expression level of the immune checkpoint is the predictive biomarker of immunotherapy response. From our results, we speculated that patients with high LMRGS scores might gain more benefit from immunotherapy. Besides, TMB associated with neoantigen production is an essential factor in driving anti-tumor immunity. High TMB increases the efficiency of stimulating host immune response (Shum et al., 2021). In our study, HCC patients with high LMRGS scores had high expression of inhibitory immune checkpoints and high TMB. Thus, the LMRGS might have a good value for precisely predicting which patients could respond to immunotherapy.

This study developed a novel LMRGS to predict the prognosis and TME in HCC. Notably, there are certain limitations in the present study. Firstly, the specific molecular functions of six genes involved in the LMRGS remain unclear. There need further experiments to elucidate the role of genes in HCC. Secondly, the LMRGS was constructed and validated using the retrospective data. In the future, we need to carry out multicenter prospective studies to validate the clinical value.

In summary, our study constructed a novel LMRGS with a high value for predicting prognosis and reflecting the TME in HCC. The LMRGS was closely associated with clinical outcomes and was an independent prognostic indicator. In addition, patients with different LMRGS scores had different TME statuses, including infiltration degree of stromal cells and immune cells, immune activity, and expression of immune checkpoints. Thus, the LMRGS was a promising biomarker to speculate molecular and immune features in HCC, which might provide new therapeutic strategies for HCC treatment.

References

• 1

  BaderJ. E.VossK.RathmellJ. C. (2020). Targeting Metabolism to Improve the Tumor Microenvironment for Cancer Immunotherapy. Mol. Cel78, 1019–1033. 10.1016/j.molcel.2020.05.034

  • CrossRef
  • Google Scholar

• 2

  BindeaG.MlecnikB.TosoliniM.KirilovskyA.WaldnerM.ObenaufA. C.et al (2013). Spatiotemporal Dynamics of Intratumoral Immune Cells Reveal the Immune Landscape in Human Cancer. Immunity39, 782–795. 10.1016/j.immuni.2013.10.003

  • CrossRef
  • Google Scholar

• 3

  BoidotR.VégranF.MeulleA.Le BretonA.DessyC.SonveauxP.et al (2012). Regulation of Monocarboxylate Transporter MCT1 Expression by P53 Mediates Inward and Outward Lactate Fluxes in Tumors. Cancer Res.72, 939–948. 10.1158/0008-5472.CAN-11-2474

  • CrossRef
  • Google Scholar

• 4

  BrownT. P.BhattacharjeeP.RamachandranS.SivaprakasamS.RisticB.SikderM. O. F.et al (2020). The Lactate Receptor GPR81 Promotes Breast Cancer Growth via a Paracrine Mechanism Involving Antigen-Presenting Cells in the Tumor Microenvironment. Oncogene39, 3292–3304. 10.1038/s41388-020-1216-5

  • CrossRef
  • Google Scholar

• 5

  CarlileT. M.MartinezN. M.SchaeningC.SuA.BellT. A.ZinshteynB.et al (2019). mRNA Structure Determines Modification by Pseudouridine Synthase 1. Nat. Chem. Biol.15, 966–974. 10.1038/s41589-019-0353-z

  • CrossRef
  • Google Scholar

• 6

  CertoM.TsaiC.-H.PucinoV.HoP.-C.MauroC. (2021). Lactate Modulation of Immune Responses in Inflammatory versus Tumour Microenvironments. Nat. Rev. Immunol.21, 151–161. 10.1038/s41577-020-0406-2

  • CrossRef
  • Google Scholar

• 7

  ChenH.JiangS.ZhangP.RenZ.WenJ. (2021). Exosomes Synergized with PIONs@E6 Enhance Their Immunity against Hepatocellular Carcinoma via Promoting M1 Macrophages Polarization. Int. Immunopharmacology99, 107960. 10.1016/j.intimp.2021.107960

  • CrossRef
  • Google Scholar

• 8

  DengH.KanA.LyuN.HeM.HuangX.QiaoS.et al (2021). Tumor-derived Lactate Inhibit the Efficacy of Lenvatinib through Regulating PD-L1 Expression on Neutrophil in Hepatocellular Carcinoma. J. Immunother. Cancer9, e002305. 10.1136/jitc-2020-002305

  • CrossRef
  • Google Scholar

• 9

  DimriM.HumphriesA.LaknaurA.ElattarS.LeeT. J.SharmaA.et al (2020). NAD(P)H Quinone Dehydrogenase 1 Ablation Inhibits Activation of the Phosphoinositide 3‐Kinase/Akt Serine/Threonine Kinase and Mitogen‐Activated Protein Kinase/Extracellular Signal‐Regulated Kinase Pathways and Blocks Metabolic Adaptation in Hepatocellular Carcinoma. Hepatology71, 549–568. 10.1002/hep.30818

  • CrossRef
  • Google Scholar

• 10

  DongL.HeY.ZhouS.CaoY.LiY.BiY.et al (2019). HIF1α-Dependent Metabolic Signals Control the Differentiation of Follicular Helper T Cells. Cells8, 1450. 10.3390/cells8111450

  • CrossRef
  • Google Scholar

• 11

  DudeckA.KöberleM.GoldmannO.MeyerN.DudeckJ.LemmensS.et al (2019). Mast Cells as Protectors of Health. J. Allergy Clin. Immunol.144, S4–S18. 10.1016/j.jaci.2018.10.054

  • CrossRef
  • Google Scholar

• 12

  FergusonL. R.ChenH.CollinsA. R.ConnellM.DamiaG.DasguptaS.et al (2015). Genomic Instability in Human Cancer: Molecular Insights and Opportunities for Therapeutic Attack and Prevention through Diet and Nutrition. Semin. Cancer Biol.35, S5–S24. 10.1016/j.semcancer.2015.03.005

  • CrossRef
  • Google Scholar

• 13

  FornerA.ReigM.BruixJ. (2018). Hepatocellular Carcinoma. The Lancet391, 1301–1314. 10.1016/S0140-6736(18)30010-2

  • CrossRef
  • Google Scholar

• 14

  FriedmanJ.HastieT.TibshiraniR. (2010). Regularization Paths for Generalized Linear Models via Coordinate Descent. J. Stat. Softw.33, 1–22. 10.18637/jss.v033.i01

  • CrossRef
  • Google Scholar

• 15

  GuilliamsM.MildnerA.YonaS. (2018). Developmental and Functional Heterogeneity of Monocytes. Immunity49, 595–613. 10.1016/j.immuni.2018.10.005

  • CrossRef
  • Google Scholar

• 16

  GyslerS. M.DrapkinR. (2021). Tumor Innervation: Peripheral Nerves Take Control of the Tumor Microenvironment. J. Clin. Invest.131, e147276. 10.1172/JCI147276

  • CrossRef
  • Google Scholar

• 17

  HarmonC.RobinsonM. W.HandF.AlmuailiD.MentorK.HoulihanD. D.et al (2019). Lactate-mediated Acidification of Tumor Microenvironment Induces Apoptosis of Liver-Resident NK Cells in Colorectal Liver Metastasis. Cancer Immunol. Res.7, 335–346. 10.1158/2326-6066.CIR-18-0481

  • CrossRef
  • Google Scholar

• 18

  IasonosA.SchragD.RajG. V.PanageasK. S. (2008). How to Build and Interpret a Nomogram for Cancer Prognosis. Jco26, 1364–1370. 10.1200/JCO.2007.12.9791

  • CrossRef
  • Google Scholar

• 19

  IzzoL. T.WellenK. E. (2019). Histone Lactylation Links Metabolism and Gene Regulation. Nature574, 492–493. 10.1038/d41586-019-03122-1

  • CrossRef
  • Google Scholar

• 20

  JhunjhunwalaS.JiangZ.StawiskiE. W.GnadF.LiuJ.MaybaO.et al (2014). Diverse Modes of Genomic Alteration in Hepatocellular Carcinoma. Genome Biol.15, 436. 10.1186/s13059-014-0436-9

  • CrossRef
  • Google Scholar

• 21

  JinZ.LuY.WuX.PanT.YuZ.HouJ.et al (2021). The Cross-Talk between Tumor Cells and Activated Fibroblasts Mediated by lactate/BDNF/TrkB Signaling Promotes Acquired Resistance to Anlotinib in Human Gastric Cancer. Redox Biol.46, 102076. 10.1016/j.redox.2021.102076

  • CrossRef
  • Google Scholar

• 22

  KanagawaM.KobayashiK.TajiriM.ManyaH.KugaA.YamaguchiY.et al (2016). Identification of a post-translational Modification with Ribitol-Phosphate and its Defect in Muscular Dystrophy. Cel Rep.14, 2209–2223. 10.1016/j.celrep.2016.02.017

  • CrossRef
  • Google Scholar

• 23

  KimJ.YuL.ChenW.XuY.WuM.TodorovaD.et al (2019). Wild-type P53 Promotes Cancer Metabolic Switch by Inducing PUMA-dependent Suppression of Oxidative Phosphorylation. Cancer Cell35, 191e8–203. 10.1016/j.ccell.2018.12.012

  • CrossRef
  • Google Scholar

• 24

  LevineA. J. (2020). p53: 800 Million Years of Evolution and 40 Years of Discovery. Nat. Rev. Cancer20, 471–480. 10.1038/s41568-020-0262-1

  • CrossRef
  • Google Scholar

• 25

  LiS.LiuY.BaiY.ChenM.ChengD.WuM.et al (2021). Ras Homolog Family Member F, Filopodia Associated Promotes Hepatocellular Carcinoma Metastasis by Altering the Metabolic Status of Cancer Cells through RAB3D. Hepatology73, 2361–2379. 10.1002/hep.31641

  • CrossRef
  • Google Scholar

• 26

  LiaoZ.-X.KempsonI. M.HsiehC.-C.TsengS.-J.YangP.-C. (2021). Potential Therapeutics Using Tumor-Secreted Lactate in Nonsmall Cell Lung Cancer. Drug Discov. Today26, 2508–2514. 10.1016/j.drudis.2021.07.014

  • CrossRef
  • Google Scholar

• 27

  LiberzonA.BirgerC.ThorvaldsdóttirH.GhandiM.MesirovJ. P.TamayoP. (2015). The Molecular Signatures Database Hallmark Gene Set Collection. Cel Syst.1, 417–425. 10.1016/j.cels.2015.12.004

  • CrossRef
  • Google Scholar

• 28

  LiuZ.LiuL.JiaoD.GuoC.WangL.LiZ.et al (2021c). Association of RYR2 Mutation with Tumor Mutation burden, Prognosis, and Antitumor Immunity in Patients with Esophageal Adenocarcinoma. Front. Genet.12, 669694. 10.3389/fgene.2021.669694

  • CrossRef
  • Google Scholar

• 29

  LiuZ.LuT.LiJ.WangL.XuK.DangQ.et al (2021b). Development and Clinical Validation of a Novel Six-Gene Signature for Accurately Predicting the Recurrence Risk of Patients with Stage II/III Colorectal Cancer. Cancer Cel Int21, 359. 10.1186/s12935-021-02070-z

  • CrossRef
  • Google Scholar

• 30

  LiuZ.LuT.LiJ.WangL.XuK.DangQ.et al (2021a). Clinical Significance and Inflammatory Landscape of aNovel Recurrence-Associated Immune Signature in Stage II/III Colorectal Cancer. Front. Immunol.12, 702594. 10.3389/fimmu.2021.702594

  • CrossRef
  • Google Scholar

• 31

  LiuZ.LuT.WangL.LiuL.LiL.HanX. (2021d). Comprehensive Molecular Analyses of a Novel Mutational Signature Classification System with Regard to Prognosis, Genomic Alterations, and Immune Landscape in Glioma. Front. Mol. Biosci.8, 682084. 10.3389/fmolb.2021.682084

  • CrossRef
  • Google Scholar

• 32

  LiuZ.WangL.LiuL.LuT.JiaoD.SunY.et al (2021e). The Identification and Validation of Two Heterogenous Subtypes and a Risk Signature Based on Ferroptosis in Hepatocellular Carcinoma. Front. Oncol.11, 619242. 10.3389/fonc.2021.619242

  • CrossRef
  • Google Scholar

• 33

  LlovetJ. M.KelleyR. K.VillanuevaA.SingalA. G.PikarskyE.RoayaieS.et al (2021). Hepatocellular Carcinoma. Nat. Rev. Dis. Primers7, 6. 10.1038/s41572-020-00240-3

  • CrossRef
  • Google Scholar

• 34

  LochnerM.BerodL.SparwasserT. (2015). Fatty Acid Metabolism in the Regulation of T Cell Function. Trends Immunol.36, 81–91. 10.1016/j.it.2014.12.005

  • CrossRef
  • Google Scholar

• 35

  LonettoG.KoifmanG.SilbermanA.AtteryA.SolomonH.Levin-ZaidmanS.et al (2019). Mutant P53-dependent Mitochondrial Metabolic Alterations in a Mesenchymal Stem Cell-Based Model of Progressive Malignancy. Cell Death Differ26, 1566–1581. 10.1038/s41418-018-0227-z

  • CrossRef
  • Google Scholar

• 36

  Loo YauH.BellE.EttayebiI.de AlmeidaF. C.BoukhaledG. M.ShenS. Y.et al (2021). DNA Hypomethylating Agents Increase Activation and Cytolytic Activity of CD8+ T Cells. Mol. Cel81, 1469e1468–1483. 10.1016/j.molcel.2021.01.038

  • CrossRef
  • Google Scholar

• 37

  MolletJ.GiurgeaI.SchlemmerD.DallnerG.ChretienD.DelahoddeA.et al (2007). Prenyldiphosphate Synthase, Subunit 1 (PDSS1) and OH-benzoate Polyprenyltransferase (COQ2) Mutations in Ubiquinone Deficiency and Oxidative Phosphorylation Disorders. J. Clin. Invest.117, 765–772. 10.1172/JCI29089

  • CrossRef
  • Google Scholar

• 38

  MüllerM.BirdT. G.NaultJ.-C. (2020). The Landscape of Gene Mutations in Cirrhosis and Hepatocellular Carcinoma. J. Hepatol.72, 990–1002. 10.1016/j.jhep.2020.01.019

  • CrossRef
  • Google Scholar

• 39

  OoH. Z.SentaniK.MukaiS.HattoriT.ShinmeiS.GotoK.et al (2016). Fukutin, Identified by the Escherichia coli Ampicillin Secretion Trap (CAST) Method, Participates in Tumor Progression in Gastric Cancer. Gastric Cancer19, 443–452. 10.1007/s10120-015-0511-2

  • CrossRef
  • Google Scholar

• 40

  Palsson-McDermottE. M.O'NeillL. A. J. (2013). The Warburg Effect Then and Now: from Cancer to Inflammatory Diseases. Bioessays35, 965–973. 10.1002/bies.201300084

  • CrossRef
  • Google Scholar

• 41

  QuinnW. J.JiaoJ.TeSlaaT.StadanlickJ.WangZ.WangL.et al (2020). Lactate Limits T Cell Proliferation via the NAD(H) Redox State. Cel Rep.33, 108500. 10.1016/j.celrep.2020.108500

  • CrossRef
  • Google Scholar

• 42

  RenaudinX.LeeM.ShehataM.SurmannE.-M.VenkitaramanA. R. (2021). BRCA2 Deficiency Reveals that Oxidative Stress Impairs RNaseH1 Function to Cripple Mitochondrial DNA Maintenance. Cel Rep.36, 109478. 10.1016/j.celrep.2021.109478

  • CrossRef
  • Google Scholar

• 43

  RingelhanM.PfisterD.O’ConnorT.PikarskyE.HeikenwalderM. (2018). The Immunology of Hepatocellular Carcinoma. Nat. Immunol.19, 222–232. 10.1038/s41590-018-0044-z

  • CrossRef
  • Google Scholar

• 44

  RosserE. C.MauriC. (2021). The Emerging Field of Regulatory B Cell Immunometabolism. Cel Metab.33, 1088–1097. 10.1016/j.cmet.2021.05.008

  • CrossRef
  • Google Scholar

• 45

  SangroB.SarobeP.Hervás-StubbsS.MeleroI. (2021). Advances in Immunotherapy for Hepatocellular Carcinoma. Nat. Rev. Gastroenterol. Hepatol.18, 525–543. 10.1038/s41575-021-00438-0

  • CrossRef
  • Google Scholar

• 46

  ShumB.LarkinJ.TurajlicS. (2021). Predictive Biomarkers for Response to Immune Checkpoint Inhibition. Semin. Cancer Biol.S1044-579X (21), 00097-3. 10.1016/j.semcancer.2021.03.036

  • CrossRef
  • Google Scholar

• 47

  SungH.FerlayJ.SiegelR. L.LaversanneM.SoerjomataramI.JemalA.et al (2021). Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A. Cancer J. Clin.71, 209–249. 10.3322/caac.21660

  • CrossRef
  • Google Scholar

• 48

  TakahashiH.KawaguchiT.YanL.PengX.QiQ.MorrisL. G. T.et al (2020). Immune Cytolytic Activity for Comprehensive Understanding of Immune Landscape in Hepatocellular Carcinoma. Cancers12, 1221. 10.3390/cancers12051221

  • CrossRef
  • Google Scholar

• 49

  TakahashiM.LioC.-W. J.CampeauA.StegerM.AyF.MannM.et al (2021). The Tumor Suppressor Kinase DAPK3 Drives Tumor-Intrinsic Immunity through the STING-IFN-β Pathway. Nat. Immunol.22, 485–496. 10.1038/s41590-021-00896-3

  • CrossRef
  • Google Scholar

• 50

  TurcotteS.KatzS. C.ShiaJ.JarnaginW. R.KinghamT. P.AllenP. J.et al (2014). Tumor MHC Class I Expression Improves the Prognostic Value of T-Cell Density in Resected Colorectal Liver Metastases. Cancer Immunol. Res.2, 530–537. 10.1158/2326-6066.CIR-13-0180

  • CrossRef
  • Google Scholar

• 51

  VidoniS.HarbourM. E.Guerrero-CastilloS.SignesA.DingS.FearnleyI. M.et al (2017). MR-1S Interacts with PET100 and PET117 in Module-Based Assembly of Human Cytochrome C Oxidase. Cel Rep.18, 1727–1738. 10.1016/j.celrep.2017.01.044

  • CrossRef
  • Google Scholar

• 52

  VillanuevaA. (2019). Hepatocellular Carcinoma. N. Engl. J. Med.380, 1450–1462. 10.1056/NEJMra1713263

  • CrossRef
  • Google Scholar

• 53

  WangH.FrancoF.HoP.-C. (2017). Metabolic Regulation of Tregs in Cancer: Opportunities for Immunotherapy. Trends Cancer3, 583–592. 10.1016/j.trecan.2017.06.005

  • CrossRef
  • Google Scholar

• 54

  WangS.PengZ.WangS.YangL.ChenY.KongX.et al (2018). KRAB-type Zinc-finger Proteins PITA and PISA Specifically Regulate P53-dependent Glycolysis and Mitochondrial Respiration. Cell Res28, 572–592. 10.1038/s41422-018-0008-8

  • CrossRef
  • Google Scholar

• 55

  XiaL.OyangL.LinJ.TanS.HanY.WuN.et al (2021). The Cancer Metabolic Reprogramming and Immune Response. Mol. Cancer20, 28. 10.1186/s12943-021-01316-8

  • CrossRef
  • Google Scholar

• 56

  YauT.KangY.-K.KimT.-Y.El-KhoueiryA. B.SantoroA.SangroB.et al (2020). Efficacy and Safety of Nivolumab Plus Ipilimumab in Patients with Advanced Hepatocellular Carcinoma Previously Treated with Sorafenib. JAMA Oncol.6, e204564. 10.1001/jamaoncol.2020.4564

  • CrossRef
  • Google Scholar

• 57

  YoshiharaK.ShahmoradgoliM.MartínezE.VegesnaR.KimH.Torres-GarciaW.et al (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun.4, 2612. 10.1038/ncomms3612

  • CrossRef
  • Google Scholar

• 58

  YuJ.ChaiP.XieM.GeS.RuanJ.FanX.et al (2021). Histone Lactylation Drives Oncogenesis by Facilitating m6A Reader Protein YTHDF2 Expression in Ocular Melanoma. Genome Biol.22, 85. 10.1186/s13059-021-02308-z

  • CrossRef
  • Google Scholar

• 59

  ZhangA.XuY.XuH.RenJ.MengT.NiY.et al (2021). Lactate-induced M2 Polarization of Tumor-Associated Macrophages Promotes the Invasion of Pituitary Adenoma by Secreting CCL17. Theranostics11, 3839–3852. 10.7150/thno.53749

  • CrossRef
  • Google Scholar

• 60

  ZhengR.WanC.MeiS.QinQ.WuQ.SunH.et al (2019). Cistrome Data Browser: Expanded Datasets and New Tools for Gene Regulatory Analysis. Nucleic Acids Res.47, D729–D735. 10.1093/nar/gky1094

  • CrossRef
  • Google Scholar

• 61

  ZhongQ.LiuZ.-H.LinZ.-R.HuZ.-D.YuanL.LiuY.-m.et al (2018). The RARS-Mad1l1 Fusion Gene Induces Cancer Stem Cell-like Properties and Therapeutic Resistance in Nasopharyngeal Carcinoma. Clin. Cancer Res.24, 659–673. 10.1158/1078-0432.CCR-17-0352

  • CrossRef
  • Google Scholar