Representing a globally major reason for death related to cancer, hepatocellular carcinoma (HCC) has a five-year overall survival (OS) rate of only 18% and is rising in incidence (1, 2). Poor survival is attributed to recurrence, relapse and treatment resistance. Specifically, 50-70% of patients with HCC experience a resurgence within five years after loco-regional therapy and up to 70% of relapsed patients recur within two years (3, 4). Treatment resistance is another common phenomenon in HCC. Immunotherapy (5) and sorafenib (6) have an objective response rate of approximately 15% and 9%, respectively, which are relatively low compared with other tumors. The efficacy is below expectation and high heterogeneity might be one of the reasons. Well-documented evidence has linked the heterogeneity of HCC with different clinical outcomes and diverse levels of sensitivity to treatment (7–9). Nevertheless, no classification accurately sub-classifying HCC patients, predicting prognosis and guiding treatment has been widely recognized in clinical practice up to now. Therefore, it is crucial to find a way to clearly distinguish HCC and elucidate biological and clinical characteristics simultaneously.

Recently, the theory of cancer stem cells (CSCs) has been generally accepted, providing new insights into cancer development and treatments. CSCs are a subgroup of tumor cells which have the potential to renew themselves and differentiate (10). A number of studies have already demonstrated the vital role of CSCs in HCC. Despite being small in number, CSCs are easily spread to distant organs, leading to HCC progression, recurrence (11) and metastasis (12). They are also implicated in therapy resistance, including sorafenib resistance (13–16), TACE refractoriness (17) and immunotherapy resistance (18). Moreover, CSCs could enhance immune evasion and induce an immunosuppressive microenvironment via the up-regulation of inhibitory molecules and the low expression of stimulatory molecules (19). They also attenuate the function of tumor-infiltrating lymphocytes by decreasing the expression of programmed cell death-ligand 1 (PD-L1), which is associated with immunotherapy response as well (20) (21). Several HCC CSC markers have been identified, including cluster of differentiation 90 (CD90), CD24, CD47, CD13, CD133, CD44, intercellular cell adhesion molecule-1 (ICAM1), epithelial cell adhesion molecule (EpCAM), leucine-rich repeat-containing receptor (LGR5), etc. (21). All the identified markers are separately reported to be tightly correlated with more aggressive HCC subtypes and worse outcomes (21). A few previous studies attempted to use CSC markers and gene signatures linked to CSC markers to identify HCC subtypes and predict HCC survival (22, 23). However, the results were unsatisfactory due to the heterogeneity of CSCs, and large-scale analysis is still needed to conclude how CSCs contribute to the prognosis of individual patients. It is necessary to conduct more research since the underlying mechanisms of CSCs remain unclear. A new stemness index, mRNAsi, was invented for the quantification of stemness on the basis of gene expression profiles (GEPs) (10). Ranging from zero (low) to one (high), the mRNAsi value was positively correlated with the similarity between tumor and stem cells. More than that, the mRNAsi value provided a new way of revealing the mechanism of CSCs in HCC, which was used as well during analysis in this article.

It has been universally accepted that energy metabolism reprogramming is a new emerging hallmark for cancer (24). The liver, which is an essential organ of metabolism, is vital in glucose and lipid homeostasis and responsible for more than 80% of protein synthesis. Diagnosed in the majority of HCC patients (68.4%) (25), metabolic associated fatty liver disease (MAFLD) has risen as one of the leading etiologies for HCC and received heated attention from researchers recently. Mounting research has shown that fatty acids, glucose, glutamine and amino acid metabolism pathways experience significant alternations in HCC owing to the energy demand for rapid cell multiplication. Glycolysis, glutamine catabolism as well as fatty acid synthesis and oxidation are enhanced, and the key enzymes included are highly expressed and related to poor clinical outcomes (26–29). Metabolic reprogramming is extremely deeply involved in the maintenance of CSCs compared with that of normal HCC cells. CSCs could preferentially survive a more restricted glucose supply by the increased expression of glucose transporter isoforms 1 (GLUT1) and GLUT3 (30). Stearoyl-coenzyme A desaturase 1 (SCD1) is an enzyme that catalyzes the desaturation of lipids, experiences a particular elevation in EpCAM⁺ alpha-fetoprotein (AFP)⁺ HCC and contributes to sorafenib resistance (31–33). CD13⁺CD133⁺HCC CSCs are deficient in Xanthine dehydrogenase/oxidase, an enzyme catalyzing purine catabolism (34). It is well-known that no studies have combined metabolism alterations with stemness features for the prediction of prognosis and immunotherapy response in HCC to date.

In this study, therefore, the transcriptomic data of patients with HCC were included to verify the following hypothesis: Stemness features were closely correlated with metabolism alterations in HCC patients, and the classification of patients based on a novel stemness-metabolism-related model could predict the clinical outcomes and immunotherapy response of HCC patients. Stemness features and mRNAsi were calculated by the OCLR algorithm. Metabolism-related gene sets were downloaded from MSigDB. Differential and weighted gene co-expression network analyses (WGCNA), univariate cox and least absolute shrinkage and selection operator (LASSO) regressions were performed and 3 most efficient prognostic DEMRGs: PFKP, PDE2A and UGT1A5 were identified. A risk score model was established and CIBERSORT, functional enrichment and copy number variation (CNV) analyses, estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE), and Tumor immune dysfunction and exclusion (TIDE) were all explored between subgroups. A novel stemness-metabolism-related model was further constructed based on the risk score model and the areas under the receiver operating characteristic (ROC) curve (AUCs) of the novel model corresponding to the survival of one, three and five years were 0.744, 0.720 and 0.680 in training dataset, and the AUCs corresponding to the survival of one, three and five years in validation dataset were 0.633, 0.769 and 0.841, respectively. Collectively, a novel stemness-metabolism-related model was proposed and the vital role of metabolism reprogramming and CSCs in predicting the clinical outcomes of patients with HCC and potential immunotherapy response was highlighted. This study presented the first comprehensive analysis of genes genetically connected with the metabolic reprogramming and stemness features of HCC, providing several promising therapeutic targets.

Heterogeneity is an essential and distinct feature of HCC and one of the main causes of poor prognosis. Staying in distinct differentiation states, CSCs maintain the ability to differentiate into diverse tumor cells, which contributes to heterogeneity. Accordingly, the subgroup classification based on CSCs might become a new viable way and shed light on future treatments. It is difficult to describe and quantify CSCs and stemness. Malta et al. adopted an innovative OCLR algorithm and defined mRNAsi as a new signature quantitatively reflecting the degree of oncogenic dedifferentiation (10). Since then, mRNAsi has provided new ideas and possibilities for linking stemness characteristics to clinical prognosis, gene mutations, treatment resistance, tumor immune characteristics, etc. Several studies have further investigated mRNAsi-related genes (55) and probed into the role of mRNAsi in HCC. Zhang et al. demonstrated a survival model using five mRNAsi-related genes (56), and Cai et al. constructed a six-gene prognosis signature (57). Mai et al. developed an HCC stemness risk model as a potential indicator of TACE treatment response (17). Zhang et al. and Xu et al. revealed the role that mRNAsi played in predicting immunotherapy response (58, 59). All studies obtained the same result that higher mRNAsi were correlated with worse outcomes and more advanced clinical stages. This is in accordance with the theory that CSCs are involved in the progression, recurrence, metastasis and treatment resistance of HCC. Because of technical limitations, not all mRNAsi-related genes are suitable and practical drug targets at present. Numerous cellular metabolic target drugs such as gemcitabine, 5-fluorouracil, l-asparaginase and methotrexate (60) have already addressed encouraging anti-cancer therapeutic effects in clinical practice or preclinical experiments. As a result, it was hypothesized that genes related to both stemness and metabolism reprogramming might be appropriate drug target candidates.

Metabolic signatures play a crucial role in tumor subclassification and immunotherapy response prediction. Some previous research had confirmed that metabolic signatures showed good power in tumor sub-phenotype and treatment response prediction in diffuse large B-cell lymphoma (61), ovarian (62) and colorectal cancers (63) as well as clear cell renal cell carcinoma (64). Chen et al. (65) made use of 90 metabolic genes to classify HCC and affirmed that metabolic signatures also showed great power in HCC subclassification. Three subtypes were developed: C1, C2 and C3 with high, low and intermediate metabolic activities, respectively. Further analysis was demonstrated, and the results showed that different subtypes also possessed different distinctions in prognostic value, immune infiltration and clinical characteristics. To be specific, C1 displayed low AFP expression and good clinical outcomes; C2 exhibited high-expression immune checkpoint inhibitors (ICIs), predicting high sensitivity to immunotherapy; C3 demonstrated high AFP expression and the worst prognosis (65). However, 90 genes for a classifier were too costly in clinical practice. In the present study, significant metabolic differences between low- and high-mRNAsi groups were confirmed. Besides, the association of metabolism reprogramming with clinical outcomes was revealed, with the enrichment of lipid-related metabolism pathways showing a relationship with poorer prognosis and carbon- and nitrogen-related metabolism pathways being closely related to better outcomes. The scope was narrowed, and a stemness-metabolic-gene signature with fewer genes was constructed. Additionally, clinical feasibility was improved while maintaining accuracy and sensitivity.

In this study, higher mRNAsi values were correlated with advanced clinical stages and worse clinical outcomes, which also aligned with previous research results (54, 59, 66). Subsequently, 74 DEMRGs were identified by performing differential analysis and WGCNA between low- and high-mRNAsi groups. LASSO and univariate cox regressions were explored, and the three most efficient prognostic genes, including PFKP, PDE2A and UGT1A5, were identified. PFKP is a protein-coding gene that translates into PFKP, an isoform of the rate-limiting enzyme of glycolysis, phosphofructokinase 1 (PFK1), which catalyzes the irreversible conversion of fructose-6-phosphate to fructose-1,6-biphosphate (67). As the important enzyme of glycolysis, PFKP was observed to allow cancer cells to survive under metabolic stress (68). Studies showed that PFKP went through an increase in HCC tissues (69) and was proven to be highly participated in glycolysis remodeling and associated with overall survival in HCC (70, 71). In addition to metabolic reprogramming, PFKP has a close correlation with stemness as well. PFPK served an essential role via LIF-Stat3 signaling to maintain embryonic stem cell (ESC) differentiation (72). The silencing of PFKP decreased the levels of stemness markers and proliferation capabilities in HCC (69). PFKP also played an important role in immune regulation. PFKP influenced stimulation of monocytes with oxidized low-density lipoprotein (73) and the expression level was correlated with interferon-gamma (IFN-γ) expression level (74). Sirtuin 2-PFKP interaction led to decreased light chain-3B activation and repressed phagocytosis (75). PFKP induced PD-L1 expression through EGFR activation and promoted immune evasion in human glioblastoma cells (76). In addition, a study explored the difference between progression HCC patients and partial response/stable HCC patients in response to the first-line combined immunotherapy, and PFKP showed a great difference in the level of mRNA, suggesting its potential in immunotherapy response stratification as well (77). PDE2A is a protein-coding gene that encodes PDE2A, an enzyme which belongs to the phosphodiesterase (PDE) family and hydrolyzes both 3’,5’-cyclic guanosine monophosphate (cGMP) and 3’,5’-cyclic adenine monophosphate (cAMP), mediating crosstalk between cGMP and cAMP signaling cascades (78), regulating mitochondrial clearance (79) and mitochondrial morphology (80, 81). The expression of PDE2A is tissue-specific and PDE2A is widely expressed in the brain and liver (78). PDE2A was demonstrated to correlate with tumorigenesis in osteosarcoma (82) and colorectal cancer (83) and to correlate with cancer stem cell stemness in glioma (84) and HCC (85). A recent study revealed that overexpressed PDE2A was associated with serum AFP level, vascular invasion, histologic grade, and pathologic stage, closely participating in inhibiting HCC cell proliferation, migration, and immune function, which had the potential to be used as a biomarker for HCC prognosis (81), the results of which consisted with our results. In glioma, PDE2A/miR-139 suppressed Wnt/β-catenin signaling by inhibiting cAMP accumulation and Glycogen Synthase Kinase-3β phosphorylation, thereby modulating the self-renewal of glioma stem cells (84). UGT1A5 is a protein-coding gene translated to UGT1A5, a member of the UGT1 family which is mainly implicated in the glucuronidation of bilirubin and phenol and acts as an essential player in the detoxification and metabolism (86, 87). Previous studies had reported UGT1A5 with a relatively low expression level in liver tissues and low enzymatic activity (88–91), while hepatic UGT1A5 expression was also proven highly inducible by multiple activators. UGT1A5 showed a significant age-dependent transcription in children (92). Treatment with rifampicin or 3-methylcholanthrene increased UGT1A5 expression level in human hepatocytes (86). In female efflux transported knockout FVB mice, the expression level of UGT1A5 was severely decreased, which indicated that UGT1A5 expression might be female-predominant at least for mice (93, 94). In the cholestasis mice model, the UGT1A5 expression level changed significantly as well (95). Those three genes were used as the basis for the establishment of a risk score model according to which patients with HCC were segmented into low- and high-risk groups. Differential analysis was performed. Moreover, it was further confirmed that DEGs were mainly involved in important metabolic pathways such as ether lipid, one-carbon and nitrogen metabolism, and disease associated with surfactant metabolism. Though the roles of PFKP, PDE2A and UGT1A5 in HCC weren’t fully understood yet, previous studies and our study suggested that the mechanisms of PFKP, PDE2A and UGT1A5 deserve further investigation. At last, the validation of the risk score model was completed in the GSE76427 dataset independently, and its AUCs corresponding to one-, three- and five-year survival were 0.740, 0.679 and 0.664, respectively, showing high predictive value.

Atezolizumab plus bevacizumab showed meaningful survival benefits in HCC, with a median OS of 19.4 months, further laying out the significance of immunotherapy in HCC treatments (96). Hence, gene signatures linked with the tumor microenvironment and the patterns of immune cell infiltration were analyzed as well, which could serve as important biomarkers to predict immunotherapy response. It has been reported that the tumor microenvironment could be conducive to maintaining CSCs which could modulate the tumor microenvironment and vice versa (21). In this study, the high-risk group demonstrated elevated immune and stromal scores, which meant that stemness features were positively associated with the abundance of stromal and immune cells. Compared with the low-risk group, the high-risk one exhibited an increased number of Tregs, M0 and other infiltrating suppressive immune cells, and higher expression levels of several immunotherapy target molecules, providing support for a negative association between stemness and immunotherapy efficacy. A recent study reported by Zhen Zhang et al. (58) demonstrated that stemness was robustly associated with ICIs through a different analysis method. This is in agreement with the findings of this research that cancer stemness was positively related to ICI resistance and intratumor heterogenicity.

Cumulative data from previous studies delineated an accurate picture of genetic variations in HCC and proved the correlation between gene alterations and antitumor immunity and metabolism (97). The genetic alterations of glucose metabolism, such as glucose-6-phosphatase, catalytic (G6PC), maturity-onset diabetes of the young type 3 (MODY3) and hepatocyte nuclear factor-1 alpha (HNF1A) genes, lead to glycogen storage diseases. That is, specific MODY3 could facilitate the occurrence of genetic liver adenomatosis and transform it into malignant HCC (98, 99). Low- and high-risk groups were found to have the same top gene mutations: tumor protein P53 (TP53), catenin beta-1 (CTNNB1) and titin (TTN) but different TMB and MSI levels and CNV patterns. The high-risk group showed more genomic alterations, indicating more possibility for immunotherapy resistance.

Finally, a nomogram making a combination of gender, age, TNM staging and gene signatures (PFKP, PDE2A and UGT1A5) was constructed for prognosis prediction and the predicted survival probability was closely fitted to the ideal line, indicating good efficacy.

CSCs are responsible for poor clinical outcomes yet are resistant to the majority of current therapies. Therapies capable of eliminating CSCs have aroused great concern, and many efforts have been dedicated to CSC-targeted therapies. CSC biomarkers are promising therapeutic targets, oncolytic measles viruses targeting CD133⁺ cells and EpCAM/CD3 and CD44 antibodies were invented in HCC (100–102). Nonetheless, no single biomarker is presented in all CSCs. For this reason, targeting a single biomarker resulted in the evasion of some CSCs (19, 103). Compared with biomarkers, CSCs share more biological features, with similar metabolic alterations. On this account, targeting altered metabolism-related pathways might be a better option and some drugs with metabolic targets have achieved inspiring results in preclinical and clinical studies. Metformin and Phenformin both restrain electron transport chain complex I and impair mitochondrial energy metabolism, giving rise to cell death in CSCs in pancreatic cancer (104). The down-regulation of the mammalian target of rapamycin (mTOR) signaling pathway which has been demonstrated to be deeply involved in energy homeostasis could reduce CSCs in breast cancer (105). The connection between metabolism reprogramming and stemness was uncovered in this study. Worthy of more in-depth research, DEMRGs have the potential to become targets for novel therapeutics. Moreover, a prognostic model containing both stemness and metabolic features was established before.

Nevertheless, limitations exist in this study. Statistical power was probably low on account of the relatively low sample size in both training and validation datasets. Therefore, the increase of statistical power makes it necessary to verify the predictive value of the novel stemness- and metabolism-related model by more HCC patients in the future. In clinical practice, BCLC staging is more widely used than TNM staging in HCC, but due to the lack of relevant clinical information on BCLC staging in public databases, our prediction model used TNM staging instead. More importantly, further experimental validations of PFKP, PDE2A, and UGT1A5 remain to be conducted at organismal, cellular and molecular levels despite powerful microarray-based bioinformatic analysis.

To sum up, the connection between metabolism reprogramming and cancer stem cells was comprehensively elucidated in this work. In addition, a survival model was established to predict prognosis and immunotherapy response with high accuracy, sensitivity and specificity. It was postulated that three genes could also be the potential therapeutic targets of HCC.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material .

YW conceived the project and wrote the manuscript. XW participated in data analysis. SD reviewed the manuscript and participated in language editing. All authors contributed to the article and approved the submitted version.