ABSTRACT
-
Background/Aims
Previously, we advocated the importance of classifying hepatocellular carcinoma (HCC) based on physiological functions. This study aims to classify HCC by focusing on liver-intrinsic metabolism and glycolytic pathway in cancer cells.
-
Methods
Comprehensive RNA/DNA sequencing, immunohistochemistry, and radiological evaluations were performed on HCC tissues from the training cohort (n=136) and validated in 916 public samples. HCC was classified using hierarchical clustering and compared with previous molecular, histopathological, and hemodynamic classifications.
-
Results
Liver-specific metabolism and glycolysis are mutually exclusive and were divided into two major subclasses: The “rich metabolism” subclass (60.3%) is characterized by enhanced bile acid and fatty acid metabolism, wellto-moderate differentiation, microtrabecular or pseudoglandular pattern, and homogeneous arterial-phase hyperenhancement (APHE), corresponding to Hoshida S3 with favorable prognosis. In IL6-JAK-STAT3-high (25.0%) conditions, upregulated ALB expression, enhanced gluconeogenesis and urea cycle activity, and an inflammatorymicroenvironment are observed. Conversely, the Wnt/β-catenin-high environment (19.9%) features elevated GLUL, APOB and CYP3A4 expression, frequent CTNNB1 (D32–S37) mutations, and an immune-desert/excluded phenotype. The “glycolysis” subclass (39.7%), characterized by histopathological dedifferentiation and downregulated liver-specific metabolism, encompasses subclasses with PI3K/mTOR (20.6%) and NOTCH/TGF-β (19.1%) signaling. These often exhibit TP53 mutations, macrotrabecular massive or compact patterns, inhomogeneous/rim-APHE, and high expression of hypoxia-inducible factors and glucose transporters, corresponding to Hoshida S1/2 with poor prognosis.
-
Conclusions
The loss of liver-specific metabolism correlates with morphological dedifferentiation, indicating cellular dedifferentiation may exhibit both physiological and pathological duality. Key signaling pathways involved in the maturation process from fetal to adult liver and zonation program may play a critical role in defining HCC diversity.
-
Keywords: Carcinoma, hepatocellular; Cell differentiation; Metabolism; Glycolysis; Gene expression profiling
Study Highlights
• This study focused on the dedifferentiation of HCC and demonstrated that it involves not only morphological changes but also physiological dedifferentiation, such as the loss of liver-specific metabolic functions. Tumor size, tumor blood flow, and the expression of hypoxia-inducible factors and glucose transporters were strongly correlated with histological and physiological dedifferentiation, which were also reflected in the tumor immune microenvironment. Furthermore, we revealed that the key signaling pathways governing the maturation of fetal to adult liver and liver zonation in normal hepatocytes may play a critical role in determining the phenotypic diversity of HCC.
Graphical Abstract
INTRODUCTION
Genomic instability is a hallmark of epithelial malignancy, caused by smoking, alcohol consumption, chronic inflammation, and aging [
1]. Through the accumulation of genetic mutations, cancer cells acquire novel functions while simultaneously losing the functional characteristics of their nonmalignant counterparts. These newly acquired capabilities include increased cell proliferation, aggressiveness, immune evasion, angiogenesis, and adaptation to competition for nutrients. As driver mutations accumulate, cancer cells gain the capacity for autonomous cell division, which markedly enhances proliferative potential [
2]. Immune evasion enables cancer cells to persist within the host, leading to the formation of immunosuppressive tumor microenvironments [
3]. As tumors grow, they develop new vascular networks and shift their metabolism, such as enhanced glucose metabolism, to overcome nutrient competition. In 1956, Warburg demonstrated that tumors generate adenosine triphosphate (ATP) through oxidative phosphorylation and by enhancing the glycolytic pathway [
4] This glycolytic shift is pronounced in gastric and lung cancers [
5], and a subset of hepatocellular carcinomas (HCC) has also been shown to exhibit enhanced glycolysis [
6]. Even in lactate-rich environments, both tumors and immunosuppressive cells continue to proliferate due to enhanced glycolysis [
5].
Conversely, cancer cells may also lose normal functions, including their well-differentiated morphology and metabolic processes. During the multistep process of hepatocarcinogenesis, which progresses from regenerative nodules to dysplastic nodules and early HCC [
7], blood flow from the portal vein and hepatic artery often decreases [
8,
9]. As early HCC dedifferentiates into moderately differentiated HCC, new blood vessels form, increasing tumor blood flow.8 This vascular increase is categorized as nonrim-arterial phase hyperenhancement (APHE) [
10]. Tumor blood flow may become heterogeneous during the progression to poorly differentiated HCC, manifesting as inhomogeneous APHE/rim-APHE, particularly when intrahepatic cholangiocarcinoma (iCCA) components are present [
11,
12]. These changes are accompanied by tumor dedifferentiation, characterized by morphological abnormalities such as disorganized cell arrangements and nuclear pleomorphism [
13]. The tumor architecture transitions from a microtrabecular pattern of well-differentiated HCC into a compact or macrotrabecular massive (MTM) pattern associated with poorly differentiated HCC. At this stage, HCC may lose hepatic-specific metabolic functions; however, research on the metabolic capabilities of tumor cells remains limited [
14-
17].
Building upon the concept of HCC classification proposed by Désert et al. [
18], which was influenced by the liver zonation program, we previously emphasized the importance of reclassifying HCC based on liver-specific metabolic functions [
19]. Traditional HCC classifications have primarily been based on cellular origin, pathological morphology, and genetic alterations, with molecular classifications, such as the Hoshida et al. [
20,
21], Calderaro et al. [
22], and Boyault et al. [
23] classifications. However, few attempts have been made to reclassify hepatocellular tumors based on their unique hepatocyte functions, including albumin production, bile acid synthesis, and lipogenesis. This reclassification provides a broader perspective on the biological characteristics of the disease.
In this study, the newly acquired glucose metabolism program (the Warburg effect) was compared with the preservation of amino acid and lipid metabolic functions typically observed in normal hepatocytes. We also evaluated the relationship between this metabolic classification and the existing classifications based on histopathological morphology, tumor hemodynamics, and molecular features.
MATERIALS AND METHODS
Patients and samples
This retrospective cohort study included 1,052 HCC patients. Of these, 136 HCC samples and their non-cancerous liver tissues were obtained at Kindai University Hospital be-tween December 2003 and July 2014 as fresh-frozen, clinicopathologically, and radiologically annotated samples from patients who underwent surgical resection of primary liver cancers. The median follow-up period was 107.3 months (95% confidence interval [CI] 93.5–121.2). HCC was diagnosed in accordance with the American Association for the Study of Liver Diseases [
24]. The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA)-LIHC (n=373) [
25], the International Cancer Genome Consortium (ICGC)-LIRI-JP (n=203), GSE14520 (n=225), and GSE76427 (n=115) cohorts were used to validate our findings.
In this study, we first performed clustering analysis using enrichment scores derived from gene sets registered in the Hallmark collection of the Molecular Signatures Database (MSigDB) (
https://www.gsea-msigdb.org/gsea/msigdb). Gene sets related to metabolism were initially analyzed, followed by those associated with signaling pathways (details are provided in the Supplementary File). In the training cohort, we identified the combination of gene sets that yielded the highest silhouette coefficient. In the validation cohort, we evaluated the validity and reproducibility of this selected combination of gene sets.
Figure 1 shows a flowchart of the study. Details of histopathological assessment, immunohistochemical staining, RNA sequencing, DNA cancer panel, and statistical analysis are provided in the online Supplemental File. This study adhered to the Declaration of Helsinki and was approved by the Kindai University Hospital IRB (Approval R04-179) in 2022. Informed consent for archived samples was waived with an opt-out approach under Japan’s Personal Information Protection Act.
RESULTS
Evaluation of tumor metabolism using MSigDB HALLMARK enrichment scores
Gene set variation analysis (GSVA) was performed on a training dataset and four external validation datasets using 50 gene sets from the Hallmark collection registered in MSigDB. Among the 50 gene sets, six were associated with metabolism (
Supplementary Table 3A). Of these, the combination of the following four gene sets yielded the highest silhouette coefficients: “HALLMARK_BILE_ACID_ METABOLISM,” “HALLMARK_FATTY_ACID_METABOLISM,” “HALLMARK_GLYCOLYSIS,” and “HALLMARK_ XENOBIOTIC_METABOLISM” (
Supplementary Table 3B).
Figure 2A presents the hierarchical clustering (Ward’s method) using the enrichment scores of four gene sets, classifying the training 136 samples into two distinct groups: a glycolysis-enhanced subclass with suppressed liver-specific metabolic functions (glycolysis subclass, n=54, 39.7%) and a subclass with enhanced liver-specific metabolic functions, including bile acid, fatty acid, and xenobiotic metabolism with suppressed glycolytic pathways (rich metabolism subclass, n=82, 60.3%). The average silhouette coefficient was 0.61 (
Supplementary Fig. 2A), and the first two principal components showed clear separation between the groups (
Supplementary Fig. 2B).
The baseline patient characteristics are summarized in
Table 1 and
Supplementary Table 4. Patients in the glycolysis subclass had a significantly larger tumor size (5.00 cm [interquartile range (IQR) 3.00–7.00] vs. 3.50 cm [IQR 2.55–5.00],
P=0.035) and a higher incidence of vascular invasion (46.3% vs. 19.5%,
P=0.002). Despite no significant differences in etiology or liver fibrosis markers, including platelet count and FIB-4 index, the glycolysis subclass exhibited higher neutrophil-lymphocyte ratio (NLR) and Creactive protein (CRP) levels. In terms of tumor markers, alpha fetoprotein (AFP) (median 235.0 vs. 7.0 ng/dL,
P<0.001) and des-γ-carboxy prothrombin (DCP) (median 227.5 vs. 115.5 mAU/mL,
P=0.018) were significantly elevated in the glycolysis subclass.
Histopathologically, the glycolysis subclass showed a significantly higher prevalence of poorly differentiated HCC (42.6% vs. 4.9%,
P<0.001), with a higher incidence of compact (48.1% vs. 21.7%,
P=0.006) and MTM patterns (17.6% vs. 0%,
P=0.001). In contrast, microtrabecular (36.0% vs. 81.1%,
P<0.001) and pseudoglandular structures were less common (1.9% vs. 32.5%,
P<0.001). Immunohistochemistry (IHC) showed a high positivity rate for biliary stem cell markers, such as epithelial cell adhesion molecule (Ep-CAM), cytokeratin 19 (CK19), and Spalt-like4 (SALL4) (
Supplementary Table 4). The Wnt/β-catenin mutation was predominantly wild-type in the glycolysis subclass (87.0%,
P=0.013), whereas
TP53/cell cycle control mutations were observed in 46.3% of cases (
P=0.066) (
Fig. 2A,
Table 1). No other genetic mutations showed significant differences between groups. Most cases in the glycolysis subclass were classified as Hoshida subclass S1/2 (98.1% vs. 17.1%,
P<0.001) (
Fig. 2A,
Table 1).
Radiological findings (
Fig. 2B) revealed a significantly higher incidence of inhomogeneous APHE/rim-APHE in the glycolysis subclass (70.8% vs. 32.4%,
P<0.001). Among 61 patients who underwent gadolinium-ethoxybenzyl-diethylenetriamine (Gd-EOB-DTPA)-enhanced MRI, the glycolysis subclass had significantly fewer nodules with a relative enhancement ratio ≥0.9 (8.0% vs. 38.9%,
P<0.05). Additionally,
18F-fluorodeoxyglucose (FDG) positron emission tomography (PET)-CT scans revealed high tumor-toliver uptake ratios (TLR) ≥2.0 in all five glycolysis subclass cases (100% vs. 0%,
P=0.02).
Survival analysis (Kaplan–Meier, log-rank test) demonstrated a median overall survival (OS) of 58.2 months (95% CI 36.9–NA) in the glycolysis subclass, compared to 92.4 months (95% CI 75.9–NA) in the rich metabolism subclass (
P=0.20) (
Fig. 2C). The median recurrence-free survival (RFS) was 10.2 months (95% CI 4.63–24.0) for the glycolysis subclass and 21.3 months (95% CI 14.6–29.3) for the rich metabolism subclass (
Fig. 2C), with a significantly higher rate of recurrence observed in the glycolysis subclass (
P=0.01).
Differential expression, gene ontology (GO)-BP, Kyoto Encyclopedia of Genes and Genomes (KEGG) gene set enrichment analysis (GSEA), and pathway mapping (
Supplementary Fig. 2C–
2G) revealed high enrichment scores for cell cycle-related pathways in the glycolysis subclass, whereas fatty acid, xenobiotics, and drug metabolism pathways showed lower scores. Pathway mapping indicated a metabolic environment favoring the conversion of glucose to pyruvate and lactate (
Supplementary Fig. 3A–
3D).
Figure 3A shows hierarchical clustering (Ward’s method) based on the signaling pathway enrichment scores for HCC samples in the rich metabolism subclass (n=82, 60.3%). The pathways of MSigDB used are listed in
Supplementary Table 5, and the top gene sets based on silhouette coefficients were HALLMARK_IL6_JAK_ STAT3_SIGNALING, HALLMARK_KRAS_SIGNALING_UP, and HALLMARK_WNT_BETA_CATENIN_SIGNALING (
Supplementary Table 3B). The silhouette coefficients for the three identified clusters are shown in
Supplementary Figure 4A. These three clusters were designated as Wnt/β-catenin-High (n=27, 19.9%), K-RAS-High (n=21, 15.4%), and IL6-JAK-STAT3-High (n=34, 25%) based on the most expressed pathway in each group.
Supplementary Table 6 compares the characteristics of the three groups. The Wnt/β-catenin-High group had a significantly higher proportion of males (92.6%,
P=0.022) and a higher frequency of
CTNNB1 (D32–S37) mutations (48.1% vs. 14.3% vs. 0%,
P=0.001). In the IL6-JAK-STAT3-High group had a smaller maximum tumor diameter (4.00 cm [IQR 3.00, 5.15] vs. 4.50 cm [IQR 3.00, 6.00] vs. 3.00 cm [IQR 2.15, 3.90],
P=0.022), and a lower frequency of pseudoglandular pattern (48.1% vs. 41.2% vs. 13.4%,
P=0.041). No significant differences were observed in etiology, NLR, CRP, tumor markers, tumor differentiation, or microtrabecular structure frequency. The K-RAS-High group showed intermediate traits between those of the other two groups. RNA analysis revealed complementary
ALB and
GLUL expression (
Fig. 3B).
GLUL expression was elevated in the Wnt/β-catenin-High group, while
ALB was higher in the IL6-JAK-STAT3-High group.
Supplementary Figure 4B illustrates immune cell proportions and tumor immune microenvironment (TIME). The IL6-JAK-STAT3-High group had significantly higher immune/stromal score calculated by ESTIMATE [
26], total immune cell counts, CD8
+ T cell counts by quanTIseq algorithm [
27,
28], and CD8
+ T cell counts by microscopy. The OS and RFS rates are shown in
Figure 3C. The median OS was 51.4 months (45.7–123.2) for the Wnt/β-catenin-High, 80.3 months (66.9–NA) for the K-RASHigh, and 156.8 months (88.1–NA) for the IL6-JAK-STAT3-High group (
P=0.06). The median PFS was 29.3 months (12.4–61.6) for the Wnt/β-catenin-High, 15.9 months (12.5– 29.3) for the K-RAS-High, and 26.1 months (14.6–NA) for the IL6-JAK-STAT3-High group (
P=0.05).
Reclassification of the glycolysis subclass
Cluster analysis reclassified HCCs into the glycolysis subclass (n=54, 39.7%) using Hallmark signaling pathways from MSigDB. The optimal classification was based on the highest silhouette coefficient when applying four pathways: HALLMARK_MTORC1_SIGNALING, HALLMARK_ NOTCH_SIGNALING, HALLMARK_PI3K_AKT_MTOR_ SIGNALING, and HALLMARK_TGF_BETA_SIGNALING (
Supplementary Table 3B). The hierarchical clustering results are presented in
Figure 4, with silhouette analysis in
Supplementary Figure 5A. The two subgroups identified were PI3K/mTOR-High (n=28, 20.6%) and NOTCH/TGF-β-High (n=26, 19.1%).
Supplementary Table 7 outlines the patient characteristics of these subgroups. The PI3K/mTOR-High group had higher proportion of males (85.7% vs. 57.7%,
P=0.046), larger tumor sizes, more frequent vascular invasion, and elevated tumor markers. No signifi-cant differences in histopathological features or genetic mutations were observed between the two groups. Immune cell infiltration, as estimated by the quanTIseq algorithm (
Supplementary Fig. 5B), showed no significant differences in CD8
+ T and regulatory T cell (Treg) counts, and ESTIMATE immune/stromal scores. The OS and RFS comparisons are shown in
Supplementary Fig. 5C. The median OS was 51.8 months for the PI3K/mTOR-High (95% CI 21.6–NA) and 97.3 months (95% CI 50.6–NA) for the NOTCH/TGF-β-High group. The median RFS was 10.2 months (95% CI 3.6–29.1) for the PI3K/mTOR-High, and 10.3 months (95% CI 6.5–36.7) for the NOTCH/TGF-β-High group, with no significant differences.
We present the patient characteristics, histopathological features, and genetic mutations for the five HCC subclasses (
Fig. 5A,
Table 2, and
Supplementary Table 8). Representative samples from each group were evaluated using hematoxylin and eosin (H&E) staining, and IHC for Ep-CAM, CK19, and carbonic anhydrase 9 (CA9) are shown in
Figure 5B,
5C. The approximate curves of enrichment scores and expression levels of metabolism-related genes (
GLUL,
CYP3A4,
APOB,
ALB) are shown in
Supplementary Figure 6A.
Supplementary Figure 6B illustrates
MYC family genes, which are involved in cell proliferation and apoptosis, hypoxia-inducible factors (
HIF1A,
CA9,
ESM1), and glucose transporters (
SLC2A1-5). A comparison of the TIME is depicted in
Supplementary Figure 6C, showing CD8
+ T cells, Tregs, and M1/M2 macrophages.
Supplementary Figure 6D provides ridge plots of the Hallmark gene set enrichment scores, and OS and RFS data are in
Supplementary Figure 6E.
The Wnt/β-catenin-High group (n=27, 19.9%) is a part of the rich metabolism subclass and has a significantly higher frequency of CTNNB1 (D32–S37) mutations. This group exhibited a higher prevalence of normal hepatocyte-like cell arrangements (microtrabecular patterns, P<0.001), and pseudoglandular structures (48.1%, P<0.001). It had the highest proportion of males (92.6%, P=0.011), and the lowest values for AFP (4.00 ng/dL [IQR 3.00, 61.50]) and DCP (90.50 mAU/mL [IQR 28.0, 213.5], both P<0.05). Tumors in this group rarely showed poorly differentiated features (14.8%, P<0.001) or positive for biliary stem cell markers EpCAM, CK19, and SALL4 (P<0.001), with 85.1% corresponding to the Hoshida S3 subclass. The ESTIMATE immune/stromal scores were the lowest, and IHC revealed low expression of immune checkpoint molecules: 0% antiprogrammed cell death protein 1 (PD-1) positivity, 11.1% programmed death-ligand 1 (PD-L1) positivity, 0% lymphocyte-activation gene (LAG)-3 positivity, and 0% T-cell immunoglobulin mucin (TIM)-3 positivity (all P<0.05). The RNA expression of GLUL and CYP3A4 was high, whereas MYC family genes, hypoxia-inducible factors, and glucose transporters were low. Hallmark enrichment scores indicat-ed higher activity in DNA repair-related pathways and low levels in inflammatory response, interferon gamma response, epithelial-mesenchymal transition (EMT), and angiogenesis.
The IL6-JAK-STAT3-High group (n=34, 25.0%), also a part of the rich metabolism subclass, had the smallest median tumor diameter at 3.00 cm (IQR 2.15, 3.90) (P=0.005), the lowest incidence of vascular invasion (14.7%, P=0.001), and relatively low tumor marker levels, with AFP at 7.00 ng/dL (IQR 3.25, 72.00) and DCP at 80.00 mAU/mL (IQR 25.5, 304.0) (P<0.05). Poorly differentiated HCC were rare (17.6%) and pseudoglandular structures were uncommon (13.4%). This group showed a higher proportion of steatohepatitis HCC (16.7%, P=0.098) and Hoshida S3 subclass (76.5%, P<0.001). Although 29.4% had Wnt/β-catenin activating mutations, none involved the D32–S37 region of CTNNB1 (P<0.001). TIME analysis indicated high ESTIMATE immune/stromal scores and M1 macrophage infiltration. RNA showed relatively high ALB, low GLUL, and CYP3A4 expression, and low MYC family genes, hypoxiainducible factors, and glucose transporters. The Hallmark scores were high for inflammatory and interferon-gamma responses, and this group had the best RFS among the five groups (P=0.02). The K-RAS-High group (n=21, 15.4%) exhibited a phenotype intermediate between the Wnt/β-catenin-High and IL6-JAK-STAT3-High groups.
The PI3K/mTOR-High group (n=28, 20.6%), a part of the glycolysis subclass, had the largest median tumor diameter (5.50 cm [IQR 3.00, 10.35], P=0.005) and a high incidence of vascular invasion (60.7%, P=0.001). Poorly differentiated HCC was common (64.3%, P<0.001), with a predominance of compact (51.8%, P=0.031) and MTM (25.9%, P=0.001) patterns, and loss of microtrabecular structures. This group had a higher frequency of TP53/cell cycle control mutations (57.1%, P=0.091) and Hoshida S1/2 subclass (96.4%, P<0.001). This group elevated NLR of 2.57 (IQR 1.73, 4.25), and the highest median values for CRP (0.53 mg/L [IQR 0.14, 4.21]), AFP (333.0 ng/dL [IQR 18.00, 4068.3]), and DCP (729.0 mAU/mL [IQR 47.5, 10714.5]) (P<0.05). RNA analysis revealed downregulation of metabolic genes (GLUL, CYP3A4, APOB, ALB) and upregulation of hypoxia-inducible factors and glucose transporters. IHC revealed CA9-positive staining. The Hallmark signatures related to DNA repair and the G2M checkpoint were significantly en-riched. PD-1 and PD-L1 positivity rates of 32.1% and 42.9%, respectively. The ESTIMATE immune/stromal scores were mixed, with both high and low scores observed within this group.
The NOTCH/TGF-β-High group (n=26, 19.1%), a part of the glycolysis subclass, had few males (57.7%) and a higher median NLR of 2.17 (IQR 1.69, 3.95). Poorly differentiated HCC was observed in 38.5%, with higher frequencies of compact (44.0%), scirrhous (12.5%), and MTM (8.3%) patterns. All patients were classified as Hoshida S1/2. Gene expression showed downregulation of metabolic genes (GLUL, CYP3A4, APOB, ALB), and upregulation of MYC family genes, hypoxia-inducible factors, and glucose transporters (SLC2A1, SLC2A3, and SLC2A5). IHC revealed CA9-positive staining. The Hallmark enrichment scores were high for DNA repair and angiogenesis. In TIME, there was reduced CD8+ T-cell infiltration, slightly increased Tregs, and upregulation of immune checkpoint molecules (PD-1: 26.9%, PD-L1: 38.5%, LAG-3: 34.6%, TIM-3: 34.6%, all P<0.05).
External validation
External validation analysis was conducted on TCGA-LIHC (n=373), ICGC-LIRI-JP (n=203), GSE14520 (n=225), and GSE76427 (n=115) cohorts. These external cohorts were classified into five HCC subclasses based on metabolic functions and signaling pathway activation, following the same methodology applied to the training cohort. The results are shown in
Figure 5D. Across all cohorts, subclasses were identified based on either enhanced glycolytic pathway activity or elevated bile acid, fatty acid, and xenobiotic metabolism. Notably, glycolysis and liver-specific metabolic functions were mutually exclusive. The rich metabolism subclass was predominantly associated with the Hoshida S3 subclass, while the glycolysis subclass corresponded mainly to the S1/S2 subclass. Furthermore, comparison of OS between the two groups revealed a trend toward poorer prognosis in the glycolysis subclass (shown in
Fig. 5E). Within the rich metabolism subclass, three distinct groups were identified based on the activation of Wnt/β-catenin, K-RAS, and IL6-JAK-STAT3 signaling pathways. The Wnt/β-catenin-High group exhibited low ESTIMATE immune and stromal scores, whereas the IL6-JAK-STAT3-High group demonstrated elevated scores. Conversely, in the glycolysis subclass, two subgroups were identified based on activation of PI3K/mTOR and NOTCH/TGF-β signaling pathways, with ESTIMATE immune and stromal scores displaying heterogeneity.
Furthermore, we performed GSEA using each dataset to examine the association between metabolic pathways. The metabolic pathways analyzed included amino acid metabolism/catabolism, gluconeogenesis, and urea cycle, which are critical for metabolic reprogramming in cancer. Our findings revealed distinct subclass-dependent alterations in these pathways. The dot plot visualization demonstrated significant variations in normalized enrichment scores (NES) and adjusted
P-values (–log10 transformed) among five subclasses, suggesting differential metabolic regulation in each subgroup (
Fig. 5F). Meta-analysis across six datasets revealed that NOTCH/TGF-β-High exhibited a suppressive trend in all pathways, whereas IL6-JAKSTAT3-High showed significant activation in gluconeogenesis, urea cycle, amino acid metabolism, and amino acid catabolism (FDR <0.05) (see
Supplementary Table 10,
Fig. 5F).
As a representative example, the analysis workflow for TCGA-LIHC is presented in
Supplementary Figure 7A–
7D. Detailed patient characteristics are shown in
Supplementary Table 9. Notably, the IL6-JAK-STAT3-High group exhibited a higher frequency of mutations in the JAK/STAT signaling genes (11.0%,
P=0.039) and a greater prevalence of HCV infection (
P=0.043).
DISCUSSION
This study elucidated that the newly acquired glucose metabolism program (glycolysis) in tumors and the maintenance of metabolic functions characteristic of normal hepatocytes were mutually exclusive. This metabolic shift closely aligns with histopathological features, such as disrupted cellular arrangement and nuclear pleomorphism. When the loss of liver-specific metabolism in cancer cells is defined as physiological dedifferentiation, the process of cellular dedifferentiation may exhibit both pathological and physiological duality (see graphic abstract). This concept, which has received limited attention, may offer a critical perspective in discussing tumor plasticity.
Pathological and physiological insights into differentiated or dedifferentiated HCC subclasses
In this study, HCC was classified into two major subtypes based on physiological-pathological differentiation and dedifferentiation. The first subtype is “rich metabolism” subclass, characterized by small tumor size, low serum tumor marker levels, a microtrabecular pattern, and well-to-moderately differentiated tumors (
Fig. 6A). These tumors correspond to Hoshida S3, with homogeneous APHE, suggesting perfusion by uniform tumor blood flow [
9]. They also exhibit higher enhancement nodules on hepatobiliary phase (HBP) images of Gd-EOB-DTPA-enhanced MRI and low
18F-FDG uptake on PET-CT, with downregulated glucose transporter gene expression. GO and KEGG analyses revealed enrichment in drug, fatty acid, and bile acid metabolism.
The second subtype is “glycolysis” subclass with suppressed liver-specific metabolic function, consisting of larger tumors, high serum AFP levels, compact or MTM types, and poorly differentiated features (
Fig. 6A). These tumors correspond to Hoshida S1/2 and are frequently associated with
TP53 mutations, upregulated cell cycle pathways, positive stem cell markers, high proliferative activity, and increased vascular invasion, reflecting their aggressive behavior. Imaging studies reveal inhomogeneous or rim-APHE, suggesting partly hypoxic blood flow regulation [
9]. These tumors exhibit hypointense in HBP on Gd-EOB-DTPA-enhanced MRI and increased
18F-FDG uptake on PETCT, with high expression of glucose transporter genes, particularly
SLC2A1, SLC2A3, and
SLC2A5. They also show elevated levels of hypoxia-inducible factors. KEGG pathway analysis demonstrated upregulated glycolysis with abundant lactate production.
Among the highly differentiated HCC subgroup, the Wnt/β-catenin-High group often has
CTNNB1 mutation in exon 3 (D32–S37 region), pseudoglandular structures, and a higher male prevalence, aligning with reports of androgen receptor influence on these structures under Wnt/β-catenin signaling [
29]. Its TIME often exhibits an immune desert/excluded phenotype, consistent with previous studies [
30-
32]. This subclass is marked by preserved metabolic functions, including ammonia (
GLUL), drug (
CYP3A4), fatty acid (
APOB), and bile acid metabolism (
Fig. 3B,
Supplementary Fig. 6A). In perivenous normal hepatocytes, Wnt/β-catenin activation drives differentiation into cells for bile acid, fatty acid, ammonia, and drug metabolism [
33].
In contrast, HCCs with activated K-RAS or IL6-JAKSTAT3 signaling show upregulated
ALB expression, enhanced gluconeogenesis and urea cycle activity, while
GLUL,
CYP3A4, and
APOB are downregulated, resembling periportal hepatocytes (
Fig. 3B,
Supplementary Fig. 6A). The IL6-JAK-STAT3-High group often includes smaller steatohepatitic-type HCCs, with immune cell-rich TIME, enhanced IFN-γ signaling, and favorable prognosis and features frequent
JAK/STAT pathway mutation. Periportal normal hepatocytes influenced by RAS signaling activation and Wnt/β-catenin signaling suppression acquire the phenotype characterized by active albumin production [
33].
Among the poorly differentiated HCC subgroup, subclasses with activated NOTCH/TGF-β signaling pathways exhibited increased CK19 expression and a higher prevalence of scirrhous-type HCC phenotypes, which resemble iCCA. Conversely, activation of the PI3K/mTOR signaling pathway was associated with features reminiscent of fetal liver and an aggressive phenotype with vascular invasion. Activation of TGF-β and NOTCH pathways is linked to increased
HNF1B and decreased
HNF4A, associated with iCCA development [
19]. HCC subclasses resembling iCCA, including combined HCC/iCCA and scirrhous-type HCC, activate these pathways [
34,
35] and show glandular structures, abundant extracellular matrix, and the expression of CK19, EPCAM, and EMT-related genes. Moreover, Boyault G3, known for MTM features, also activates the MAP Kinase and PI3K/AKT pathways, similar to the Boyault G1/G3 [
22].
We previously emphasized the importance of reclassifying HCC by tumor metabolic functions, highlighting that similar signaling pathways regulate both HCC and normal hepatobiliary cell phenotypes via the liver zonation program [
19]. These observations align with prior reports on HCC classification based on liver zonation principles [
18]. However, our study extends this framework by incorporating an evaluation of tumor metabolic function. Interestingly, our study also demonstrated that HCC subclasses closely resemble hepatocellular adenoma (HCA). HCCs with activated Wnt/β-catenin signaling shared features with
CTNNB1-mutated HCAs, including bile acid and ammonia metabolism. Moreover, HCCs with activated IL6-JAK-STAT3 signaling pathways resembled inflammatory HCAs in their specific signaling pathway and TIME (
Supplementary Fig. 8).
We revealed that the signaling pathways involved in maturation of fetal to adult liver and liver zonation in normal hepatocytes may play a key role in determining the diversity of HCC phenotypes. This framework provides a compelling perspective for understanding the spatial and temporal heterogeneity of HCC. Results from all five cohorts collectively summarize the defining characteristics of these five HCC subclasses (
Fig. 6B).
Cancer cells generate ATP not only through mitochondrial oxidative phosphorylation but also via glycolysis, even under aerobic conditions [
4]. This Warburg effect is characterized by high glucose transporter expression and enhanced glycolysis. Glucose is metabolized into fructose-6-phosphate and fructose-1,6-bisphosphate, producing abundant pyruvate and lactate. The glycolysis subclass exhibited high expression of hypoxia-inducible factors (e.g.,
HIF1A,
CA9), and glucose transporters (e.g., GLUT1 [
SLC2A1], GLUT3 [
SLC2A3], GLUT5 [
SLC2A5]) (
Supplementary Fig. 6B). KEGG pathway mapping revealed activation of key glycolytic enzymes, including hexokinase, phosphofructokinase, glyceraldehyde-3-phosphate dehydrogenase, phosphoglycerate kinase 1/2, phosphoglycerate mutase 1/2, enolase 1/2/3, and pyruvate kinase M1/2 (
Supplementary Fig. 3A,
3B), indicating rapid glucose conversion to pyruvate in tumor cells. This metabolic environment also affects TIME, creating a survival environment favorable to immunosuppressive cells such as Tregs [
5].
Metastatic liver tumors exhibit increased expression of glycolysis-related molecules via upregulation of
HIF1A expression, leading to elevated lactate levels [
35]. Hypoxia triggers metabolic reprogramming, shifting from oxidative phosphorylation to glycolysis.
HIF1A activation, often due to cancer-related mutations, promotes therapeutic resistance [
36] and upregulates the aldolase A, a key enzyme that supports tumor growth [
37]. The correlation between glucose transporters, particularly GLUT1 and GLUT3, and
18F-FDG uptake in HCC is well established, with high accumulation detectable by
18F-FDG/PET-CT [
38-
40]. Gene expression related to glycolysis and the HIF1A signaling also correlates with
18F-FDG uptake, consistent with our findings [
6]. Enhanced glycolysis and lipid metabolism are mutually exclusive, as observed in our study [
6]. Additionally, increased expression of genes related to the mTOR pathway correlates with
18F-FDG uptake [
41], consistent with the PI3K/mTORHigh subclass in our cohort. These findings suggest that HCC with upregulated glucose transporter expression and activated glycolysis may correspond to PET-positive HCC. However, future research on HCC metabolic functions with high
18F-FDG uptake is required.
Despite these important considerations, several limitations exist. First, although we reclassified HCC based on metabolic functions and signaling pathways and demonstrated their correlations, we could not establish a causal relationship. Specifically, we did not examine whether activation of these signaling pathways directly induces phenotypic changes in HCC cells. Second, although we showed a correlation between the preservation of metabolic functions and suppressed glycolysis and conversely, the loss of metabolic functions with enhanced glycolysis and increased hypoxia-inducible factor expression, the mechanisms linking these processes remain unclear, and the precise pathways driving metabolic changes in HCC were not explored. Third, the number of 18F-FDG/PET-CT cases was limited, and the association between glycolysis and the imaging findings was not thoroughly investigated. Further studies are required to elucidate these relationships.
Conclusion
In HCC, pathological dedifferentiation is associated with the loss of liver-specific metabolism and a shift towards a glycolysis-dominant metabolic environment. This shift is driven by the absence of homogeneous tumor blood flow and increased expression of hypoxia-inducible factors and glucose transporters. This pathological and physiological dedifferentiation stratifies HCC into two major subtypes: well-to-moderately differentiated HCC, in which key signaling pathways involved in liver zonation define the phenotype, and dedifferentiated HCC, which exhibits tumor characteristics resembling fetal liver or biliary epithelial cells. These findings highlight the role of hepatocyte-specific maturation processes in shaping the phenotypic diversity of HCC.
FOOTNOTES
-
Authors’ contribution
Conceptualization, T.A. and N.N.; Methodology, T.A. and N.N.; Software, T.A.; Validation, T.A., N.F. and Y.K.; Formal Analysis, T.A., K.S. and K.N.; Investigation, T.A., N.N. and M. Takita; Data Curation, M. Tsurusaki, K.H., M.M., H.C., M. Takita, S.H., H.I, K.U., Y.M, A.T., K.K., T.M., T.N., and I.M.; Writing – Original Draft Preparation, T.A. and N.N.; Writing – Review & Editing, T.A., M.T., K.H., and N.N.; Visualization, T.A. and M. Tsurusaki; Supervision, M.K.; Project Administration, T.A. and M.K.
-
Acknowledgements
The authors thank Yuji Nishikawa (Asahikawa Medical University, Japan), Atsushi Miyajima, Tohru Itoh (Laboratory of Cell Growth and Differentiation, Institute for Quantitative Biosciences, The University of Tokyo, Japan), Michiie Sakamoto (School of Medicine, International University of Health and Welfare, Japan), Koji Kadota (Department of Chemistry and Biotechnology, Graduate School of Engineering, The University of Tokyo, Japan), Takashi Ishikura (Field Bioinformatics Specialist, Technical Support, Thermo Fisher Scientific, Japan), and Toshio Kitamura (Division of Cellular Therapy, IMS, The University of Tokyo, Japan) for their valuable contributions, expertise, and assistance throughout the study. We acknowledge the use of AI in writing R scripts for data analysis and English language editing services provided. M. Tsurusaki. We express our gratitude to the editors and reviewers for their insightful comments and suggestions, which have significantly improved the quality of this manuscript.
This work was supported in part by a Grant-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (KAKENHI: 21K07184, N. Nishida), (24K10393, N. Nishida), and a grant from SRF (T. Aoki and N. Nishida).
-
Conflicts of Interest
Y.K. received research grant from Eisai. K.S. received lecture fee from Life Technologies Japan Ltd., Yodosha Co., Ltd., Qiagen, Inc., Takeda Pharmaceutical Co., Ltd., and Nippon Kayaku Co., Ltd. K.U. received honoraria from Eisai, Takeda and Chugai Pharmaceuticals. K.N. received grants from Nichirei Biosciences Inc., West Japan Oncology Group, Eli Lilly Japan K.K., Hitachi, Ltd., Nippon Boehringer Ingelheim Co., Ltd., osakaminami hospital, Sysmex Corporation, Otsuka Pharmaceutical, Thoracic Oncology Research Group, University Public Corporation Osaka, and Okayama University. He has also received lecture fee from AstraZeneca K.K., Chugai Pharmaceutical Co., Ltd., DAIICHI SANKYO, ELI LILLY JAPAN K.K., Guardant Health Inc., Invitae Japan, Janssen Pharmaceutical K.K., Maruho Co., Ltd., MSD K.K., Nichirei Biosciences Inc., Novartis Pharma K.K., Ono Pharmaceutical Co., Ltd., Otsuka Pharmaceutical Co., Ltd., and SymBio Pharmaceuticals Limited. M.K. received grants from Taiho Pharmaceuticals, Chugai Pharmaceuticals, Otsuka, Takeda, GE Healthcare, AbbVie, Astellas Pharma, Eisai, Eli Lilly, and AstraZeneca. He has also received grants and personal fees from Eisai, Roche, Chugai, and AstraZeneca. The other authors disclose no conflicts.
SUPPLEMENTAL MATERIAL
Supplementary material is available at Clinical and Molecular Hepatology website (
http://www.e-cmh.org).
Supplementary Table 3A. Metabolism gene sets listed in the Hallmark collection of MSigDB
Supplementary Table 3B. Gene sets of Hallmark metabolism and signaling
cmh-2024-1088-Supplementary-Table-3.pdf
Supplementary Table 10A. Meta-analysis results of GSEA based on the hepatocyte functions across six independent datasets. Adjusted P-value of amino acid metabolism
Supplementary Table 10B. Meta-analysis results of GSEA based on the hepatocyte functions across six independent datasets. Adjusted P-value of amino acid catabolism
Supplementary Table 10C. Meta-analysis results of GSEA based on the hepatocyte functions across six independent datasets. Adjusted P-value of gluconeogenesis
Supplementary Table 10D. Meta-analysis results of GSEA based on the hepatocyte functions across six independent datasets. Adjusted P-value of urea cycle
cmh-2024-1088-Supplementary-Table-10.pdf
Supplementary Figure 1.
Nearest template prediction of the Hoshida classification. Generation of GCT files from RNA-seq data was performed, including (A) 136 training samples, (B) TCGA-LIHC (n=373), (C) ICGC-LIRI-JP (n=203), (D) GSE14520 (n=225), and (E) GSE76427 (n=115). NTP analysis was conducted, and the results are displayed as sample bar plots, heatmaps, and BH-FDR-adjusted values.
Supplementary Figure 2.
(A) Result of silhouette analysis for glycolysis subclass vs. rich metabolism subclass Plot of silhouette scores. Cluster 1 (n=82) represented the metabolism-rich subclass, whereas cluster 2 (n=54) represented the glycolysis subclass. (B) The scatter plot of the first and second principal components. Principal component analysis (PCA) plots with the first and second principal components as axes. Samples plotted in green represent the rich metabolism subclass (n=82) and red represent the glycolysis subclass (n=54). (C) The result of expression variability analysis. Differentially expressed genes (DEGs) were identified in the glycolysis and rich metabolism subclasses using the TCC package with a false discovery rate (FDR) cutoff of 0.05. An MA plot of expression variability analysis is shown. DEGs are highlighted in magenta (9,959 genes), whereas non-differentially expressed genes (non-DEGs) are shown in black. (D) GO enrichment analysis of differentially expressed genes identified by TCC. A total of 9,959 differentially expressed genes were identified using TCC, and a gene vector was created for GO analysis by selecting genes with an absolute Log2FC greater than 2. The P-value and q-value cutoffs were set at 0.05, and the top five biological process (BP) categories were listed. The x-axis represents the number of genes involved in each biological process. The color of the bars corresponds to the P-values adjusted for multiple testing using the Benjamini-Hochberg method (p.adjust). (E) Ridge plot of GSEA for GO-BP analysis. Gene set enrichment analysis. A GSEA was performed using a ranked gene list to identify the enriched biological processes (BP in the GO database). The analysis included gene sets with sizes ranging from 20 to 500 genes. The P-value cutoff was set at 0.05. The top five enriched GO-BP categories are shown in the ridge plot. Each ridge represents the distribution of ranked genes for a particular GO-BP category, with the height indicating gene density. The x-axis denotes the enrichment score (ES), which reflects the degree to which genes in the GO category were overrepresented at the top or bottom of the ranked gene list. The bars are color-coded based on the level of statistical significance. (F) The result of KEGG gene set enrichment analysis. The bar plot illustrates the top 5 upregulated and top 5 downregulated KEGG pathways based on normalized enrichment scores (NES) from Gene Set Enrichment Analysis (GSEA). The y-axis represents the KEGG pathway descriptions, and the xaxis shows the NES values. Bars are colored according to the q-value, with a gradient from red (high q-value) to blue (low q-value). (G) Dot plot of KEGG pathway enrichment analysis. This plot shows the top 15 KEGG pathways enriched in each gene set. The x-axis represents the gene ratio (the proportion of genes in the pathway relative to the total number of genes in that pathway) and the y-axis lists the KEGG pathway descriptions. Each dot represents a KEGG pathway and is colored according to the -log10(P-value), indicating the statistical significance of enrichment. Larger dots denote pathways with a higher gene ratio, whereas color intensity reflects the significance level.
Supplementary Figure 3.
(A) KEGG pathway mapping of glycolysis/glycogenesis (all 20,802 genes). This figure illustrates the glycolysis/ gluconeogenesis pathway (KEGG ID: “hsa00010”), focusing on changes in gene expression. The gene expression log fold change (logFC) was mapped to the pathway. All 20,802 genes were used, with a color representation of the gene expression levels observed in the glycolysis subclass. The color and size of the gene nodes within the pathway correspond to logFC. The color gradient indicates the extent of the expression increase or decrease, with warmer colors representing higher expression levels and cooler colors representing lower expression levels. Compound nodes (metabolites) in the pathway were displayed using EC numbers. Please refer to the list below for comparison.
HK3; hexokinase 3 [KO:K00844] [EC:2.7.1.1]
HK1; hexokinase 1 [KO:K00844] [EC:2.7.1.1]
HK2; hexokinase 2 [KO:K00844] [EC:2.7.1.1]
HKDC1; hexokinase domain containing 1 [KO:K00844] [EC:2.7.1.1]
GCK, glucokinase [KO:K12407] [EC:2.7.1.2]; PFKM, phosphofructokinase; muscle [KO:K00850] [EC:2.7.1.11]
PFKP; phosphofructokinase, platelet [KO:K00850] [EC:2.7.1.11]
PFKL, phosphofructokinase; liver type [KO:K00850] [EC, 2.7.1.11]; GAPDH, glyceraldehyde-3-phosphate dehydrogenase [KO:K00134] [EC:1.2.1.12]
PGK2; phosphoglycerate kinase 2 [KO:K00927] [EC:2.7.2.3].
PGK1; phosphoglycerate kinase 1 [KO:K00927] [EC:2.7.2.3]
PGM1; phosphoglucomutase 1 [KO:K01835] [EC:5.4.2.2]
PGAM1; phosphoglycerate mutase 1 [KO:K01834] [EC:5.4.2.11]
PGAM2; phosphoglycerate mutase 2 [KO:K01834] [EC:5.4.2.11]
FBP1; fructose-bisphosphatase 1 [KO:K03841] [EC:3.1.3.11]
FBP2; fructose-bisphosphatase 2 [KO:K03841] [EC:3.1.3.11]
ALDOC, aldolase, fructose bisphosphate C [KO:K01623] [EC:4.1.2.13]
ALDOA: aldolase, fructose bisphosphate A [KO:K01623] [EC:4.1.2.13].
ALDOB, aldolase, fructose bisphosphate B [KO:K01623] [EC:4.1.2.13],
ENO3, enolase 3 [KO:K01689] [EC:4.2.1.11]
ENO2; enolase 2 [KO:K01689] [EC:4.2.1.11]
ENO1; enolase 1 [KO:K01689] [EC:4.2.1.11]
ENO4, enolase 4 [KO:K27394] [EC:4.2.1.11]; PKM, pyruvate kinase M1/2 [KO:K00873] [EC:2.7.1.40]
PKLR; pyruvate kinase L/R [KO:K12406] [EC:2.7.1.40]
LDHA; lactate dehydrogenase A [KO:K00016] [EC:1.1.1.27]
LDHB; lactate dehydrogenase B [KO:K00016] [EC:1.1.1.27]
LDHC; lactate dehydrogenase C [KO:K00016] [EC:1.1.1.27]
(B) KEGG pathway mapping of glycolysis/glycogenesis (9,959 DEGs). The figure illustrates the glycolysis/gluconeogenesis pathway (KEGG pathway ID: “hsa00010”), with gene expression log fold changes (logFC) mapped onto the pathway. Gene expression levels in the glycolysis subclass were color-coded based on the 9,959 differentially expressed genes identified by TCC. (C) KEGG pathway mapping of pyruvate metabolism (all 20,802 genes). The figure depicts the pyruvate metabolism pathway (KEGG pathway ID: “hsa00620”), with gene expression log fold changes (logFC) mapped onto the pathway. Gene expression levels in the glycolysis subclass were colorcoded based on all 20,802 genes, including both differentially expressed genes and non-differentially expressed genes. (D) KEGG pathway mapping of HIF signaling pathway (all 20,802 genes). The figure illustrates the HIF signaling pathway (KEGG pathway ID: “hsa04066”), with gene expression log fold changes (logFC) mapped onto the pathway. Gene expression levels in the glycolysis subclass were color-coded based on all 20,802 genes, including both differentially expressed and non-differentially expressed genes.
cmh-2024-1088-Supplementary-Figure-3.pdf
Supplementary Figure 4.
(A) Result of silhouette analysis for Wnt/b-catenin-High vs. K-RAS-High vs. IL6-JAK-STAT3-High group including rich metabolism subclass. Cluster 1 represents Wnt/β-catenin-High group (n=27), Cluster 2 represents IL6-JAK-STAT3-High group (n=34), and Cluster 3 represents K-RAS-High group (n=21). (B) Tumor immune microenvironment for Wnt/b-catenin-High vs. KRAS-High vs. IL6-JAK-STAT3-High group including rich metabolism subclass. Violin plots of immune cell infiltration estimates were obtained using the MCP-counter and CIBERSORT as well as the immune and stromal scores from ESTIMATE.
Supplementary Figure 5.
(A) Result of silhouette analysis for PI3K/mTOR-High vs. NOTCH/TGF-b-High group including glycolysis subclass. Cluster 1 represents PI3K/mTOR-High group (n=28), and Cluster 2 represents NOTCH/TGF-β-High group (n=26). (B) Tumor immune microenvironment for PI3K/mTOR-High vs. NOTCH/TGF-b-High group including glycolysis subclass. Violin plots of immune cell infiltration estimates were obtained using the MCP-counter and CIBERSORT as well as the immune and stromal scores from ESTIMATE. (C) Overall survival and recurrence free survival for PI3K/mTOR-High vs. NOTCH/TGF-b-High group including glycolysis subclass. The Kaplan-Meier survival curves are presented for three clusters based on distinct signaling pathways: NOTCH/TGF-b-High (n=26, red) and PI3K/mTOR-High (n=28, blue). Censoring is indicated on the curves and 95% confidence intervals are displayed. A log-rank test was performed to compare the survival distribution among clusters. The median survival and recurrence-free times were calculated, and a risk table is included in the figure to show the number of patients at risk over time.
Supplementary Figure 6.
(A) Approximate curve for Hallmark enrichment scores and RNA expression levels of metabolism related genes. Curves approximating the enrichment scores and RNA expression levels (Z-scores) for each sample were plotted. (B) Approximate curve for RNA expression levels of hypoxia induced factors and glucose transporters. Curves approximating RNA expression levels (Z-scores) for each sample were plotted. In the glycolysis subclass, the expression levels of the MYC family members (MYCN with 95%C. I.), hypoxia-inducible factors (HIF1A with 95% C.I. and CA9), and glucose transporters (SLC2A1 with 95%C.I., SLC2A3, and SLC2A5) were relatively high. The heatmap was scaled by column, with color gradients indicating enrichment scores for each metabolism, ranging from low values (dark green) to high values (yellow). (C) Violin plots for immune cell infiltration. Immune cell infiltration ratios estimated by quanTIseq and ESTIMATE immune/stromal scores. (D) Ridge plots for Hallmark enrichment scores. Ridge plots of enrichment scores for hallmark gene sets of MSigDB, showing those with significant differences among the five groups, as determined by the Mann-Whitney U test (P<0.05). (E) Overall survival and recurrence-free survival for five groups. Kaplan-Meier survival curves are presented for the five HCC subclasses. Censoring is indicated on the curves, and survival distributions among clusters were compared using the log-rank test. A risk table is included in the figure showing the number of patients at risk over time. Among the five subgroups, the IL6-JAK-STAT3-High group demonstrated the best recurrence-free survival (RFS) with statistical significance.
Supplementary Figure 7.
(A) External validation using TCGA-LIHC: classification based on tumor metabolic functions. In TCGA-LIHC (n=373), we first performed hierarchical clustering analysis (Ward’s method) based on the enrichment scores of metabolism-related gene sets from the hallmark collection of MSigDB. The results are shown in Supplementary Figure 7A, revealing distinct groups: one characterized by enhanced glycolysis and suppression of other metabolic functions (glycolysis subclass, n=163, 43.7%), and another characterized by suppressed glycolysis and preserved hepatocyte metabolic functions (rich metabolism subclass, n=210, 56.3%). The average silhouette coefficient was 0.56, indicating moderate clustering quality. The heatmap was scaled by column, with color gradients indicating enrichment scores for metabolism, ranging from low values (dark green) to high values (yellow). (B) External validation using TCGA-LIHC: reclassification of the rich metabolism subclass by activated signaling pathways. Next, the rich metabolic subclasses were stratified based on the activation of specific signaling pathways. The 210 cases classified into the rich metabolism subclass were divided into groups with activated Wnt/β-catenin signaling (n=64, 17.2%), activated K-RAS signaling (n=24, 6.4%), and activated IL6-JAK-STAT3 signaling (n=122, 32.7%). The heatmap is scaled by columns, with color gradients indicating enrichment scores for each signaling pathway, ranging from low values (dark green) to high values (yellow). (C) External validation using TCGA-LIHC: reclassification of the glycolysis subclass by activated signaling pathways. Among the 163 cases classified into the glycolysis subclass from TCGA-LIHC, the cases were further divided into two groups based on the enrichment of the four signaling pathways. The subclasses were classified according to the activation of the PI3K/mTOR (n=89, 23.9%) and NOTCH/TGF-β signaling pathways (n=74, 19.8%). The heatmap is scaled by columns, with color gradients indicating enrichment scores for each signaling pathway, ranging from low values (dark green) to high values (yellow). (D) External validation using TCGA-LIHC: reclassification of the glycolysis subclass by activated signaling pathways. External Validation Using TCGA-LIHC cohort: MSigDB enrichment scores, genetic mutations, and ESTIMATE immune/stromal scores are shown in a heatmap. The heatmap is scaled by column, with color gradients indicating enrichment scores for each signaling pathway, ranging from low values (dark green) to high values (yellow). Gene mutations are indicated in dark blue for positive cases and white for negative cases. The ESTIMATE stromal/immune scores are displayed using a color scale, with warmer colors indicating higher scores and cooler colors indicating lower scores. HCC, hepatocellular carcinoma; TGFB, Transforming Growth Factor-β; PI3K/mTOR, Phosphoinositide 3-kinase/ Mechanistic target of rapamycin; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumors; H&E, Hematoxylin and Eosin stain; CA9, Carbonic Anhydrase IX; CK19, cytokeratin 19; EpCAM, anti-epithelial cell adhesion molecule; TCGA-LIHC, The Cancer Genome Atlas Liver Hepatocellular Carcinoma.
Supplementary Figure 8.
Homology between HCC, HCA, and normal hepatocytes. Within the rich metabolism subclass, Wnt/β-catenin-High and IL6-JAK-STAT3/K-RAS-High groups exhibit tumor phenotypes similar to CTNNB1-mutated HCA and inflammatory HCA, respectively. In normal hepatocytes, bile acid production and GLUL expression are upregulated following Wnt/β-catenin signaling activation. Ha-RAS signaling activation enhances albumin production. The combination of signaling pathways and phenotypes observed in normal hepatocytes closely resembles those seen in HCC and HCA. HCC, hepatocellular carcinoma; RNA, ribonucleic acid; HCA, hepatocellular adenoma; IL6-JAK-STAT3, Interleukin-6 - Janus Kinase - Signal Transducer and Activator of Transcription 3; TCGA-LIHC, The Cancer Genome Atlas Liver Hepatocellular Carcinoma; TIME, tumor immune microenvironment; AFP, α-fetoprotein; DCP, des-γ-carboxy prothrombin; CA9, Carbonic Anhydrase IX; HIF, Hypoxia Inducible Factor.
Figure 1.Study design. The study comprised a training cohort of 136 HCC samples and validation cohorts of 916 samples. Metabolic gene sets from 50 Hallmark gene sets in MSigDB were analyzed in the training cohort using hierarchical clustering and silhouette analysis. HCC was further subclassified based on signaling pathway-related gene sets and enrichment scores. DNA mutation analysis, as well as pathological and radiological evaluations, was conducted on the same cohort. APHE, arterial-phase hyperenhancement; CECT, Contrast-Enhanced Computed Tomography; CT, computed tomography; FDG, fluorodeoxyglucose; Gd-EOB-DTPA, gadolinium-ethoxybenzyl-diethylenetriamine; GSVA, gene set variation analysis; HBP, hepatobiliary phase; HCC, hepatocellular carcinoma; ICGC, International Cancer Genome Consortium; IHC, immunohistochemistry; PET, positron emission tomography; RER, relative enhancement ratio; SUV, standardized uptake value; TCGA-LIHC, The Cancer Genome Atlas Liver Hepatocellular Carcinoma; TIME, tumor immune microenvironment; TLR, tumor-to-liver uptake ratios.
Figure 2.Subclassification of HCC based on tumor metabolic function. (A) Heatmap and hierarchical clustering analysis. Hierarchical clustering of 136 samples based on enrichment scores from 4 metabolic gene sets (fatty acid, bile acid, xenobiotic metabolism, glycolysis). Samples were divided into the glycolysis subclass (red band) and rich metabolism subclass (green band). Associated Hoshida subclass, Wnt/β-catenin, and TP53/cell cycle mutations are displayed. (B) CT/MRI APHE image. Representative CT/MRI and 18F-FDG/PETCT images show glycolysis subclass with inhomogeneous/rim-APHE and high 18F-FDG uptake (TLR ≥2.0), while the rich metabolism subclass exhibits homogeneous APHE and no 18F-FDG uptake (TLR <2.0). Gd-EOB-DTPA-enhanced MRI reveals higher HBP enhancement in the rich metabolism subclass. (C) Overall survival and recurrence-free survival (Kaplan–Meier, log-rank test). Kaplan–Meier curves for OS and RFS compare the glycolysis subclass (n=54, red) and the rich metabolism subclass (n=82, green). APHE, arterial phase hyperenhancement; CT, computed tomography; Gd-EOB-DTPA, gadolinium-ethoxybenzyl-diethylenetriamine; HBP, hepatobiliary phase; HCC, hepatocellular carcinoma; MRI, magnetic resonance imaging; MSigDB, molecular signature database; OS, overall survival; RFS, recurrence-free survival; 18F-FDG/PET-CT, 18F-fluorodeoxyglucose positron emission tomography-computed tomography.
Figure 3.Subclassification of rich metabolism subclass HCC. (A) Heatmap and hierarchical clustering analysis. This heatmap shows enrichment scores of signaling pathways in the rich metabolism subclass. Metabolic gene expression is scaled as Z-scores, with a color gradient from dark green (low) to yellow (high). Tumor-infiltrating immune cell proportions (quanTIseq), gene mutation data (white to light slate blue), and ESTIMATE stromal/immune scores (warm to cool colors) are also displayed. (B) Violin plot. The expression levels of GLUL and ALB were compared across the three groups. GLUL expression increased with Wnt/β-catenin activation, and ALB expression with IL6-JAK-STAT3 activation. (C) Overall survival and recurrence-free survival (Kaplan–Meier, log-rank test). Kaplan–Meier curves compare three clusters: Wnt/β-catenin-High (n=27, red), IL6-JAK-STAT3-High (n=34, blue), and K-RAS-High (n=21, green). IL6-JAKSTAT3-High showed better overall (P=0.06) and recurrence-free survival (P=0.05). Risk tables and 95% confidence intervals are included. ESTIMATE, estimation of stromal and immune cells in malignant tumors; HCC, hepatocellular carcinoma; IL6-JAK-STAT3, interleukin-6-janus kinase-signal transducer and activator of transcription 3; OS, overall survival; RFS, recurrence-free survival.
Figure 4.Subclassification of glycolysis subclass HCC. Heatmap and hierarchical clustering analysis. The heatmap classifies the glycolysis subclass with high silhouette coefficients, further divided into NOTCH/TGF-β-High (n=26) and PI3K/mTOR-High (n=28) via Ward.D2 clustering. Enrichment scores are scaled by column, with a gradient from dark green (low) to yellow (high). Tumor-infiltrating immune cell proportions (quanTIseq), mutation data (white to light slate blue), and ESTIMATE stromal/immune scores (cool to warm colors) are also shown. ESTIMATE, estimation of stromal and immune cells in malignant tumors; HCC, hepatocellular carcinoma; PI3K/mTOR, phosphoinositide 3-kinase/mechanistic target of rapamycin; TGFB, transforming growth factor-β.
Figure 5.Comparison of five subclasses. (A) Heatmap: Five subclasses are shown with MSigDB enrichment scores, Hoshida classifications, mutations, TIME (quanTIseq, ESTIMATE), and histopathological data. Color gradients represent enrichment scores (dark green to yellow). Histological classifications (H&E staining) are marked in blue; necrotic samples are gray. (B) Histopathological evaluation: Representative H&E, EpCAM, and CK19 staining. Glycolysis subclasses (NOTCH/TGF-β-High, PI3K/mTOR-High) often show MTM/compact structures with EpCAM/CK19 positivity, while rich metabolism subclasses retain cord-like structures with negative staining. (C) Immunohistochemistry: Glycolysis subclasses show strong CA9 positivity, indicating hypoxia, while rich metabolism subclasses are negative. (D) Comparison with external cohorts: Heatmaps classify training and validation cohorts into glycolysis (subclasses 1, 2) and rich metabolism (subclasses 3–5). Metabolic enrichment scores (low in dark purple, high in orange), signaling pathway (low in dark green, high in yellow), Hoshida subclasses, and immune/stromal scores (low in blue, high in red) are displayed. Subclass 1 corresponds to NOTCH/TGF-β-High, subclass 2 to PI3K/mTOR-High, subclass 3 to IL6-JAK-STAT3-High, subclass 4 to K-RAS-High, and subclass 5 to Wnt/β-catenin-High. (E) Overall survival of external cohorts: Kaplan–Meier curves show better OS for the rich metabolism subclass (green) compared to glycolysis (red). (F) Meta-analysis results of GSEA based on metabolic pathways among the training and validation cohorts: The dot plot visualizes the results of GSEA for amino acid metabolism/catabolism, gluconeogenesis, and urea cycle across five subclasses. The x-axis represents the subclasses: S1 corresponds to NOTCH/TGF-β-High, S2 to PI3K/mTOR-High, S3 to IL6-JAK-STAT3-High, S4 to KRAS-High, and S5 to Wnt/β-catenin-High. Meanwhile, the y-axis shows the adjusted P-value on a –log10 scale, where positive values indicate activation and negative values suggest suppression of the pathway. The dots represent each dataset, and their size corresponds to the NES. Black and gold diamonds indicate the mean NES for each subclass. The analysis was conducted by comparing each subclass against all others (e.g., subclass 1 vs. others, subclass 2 vs. others, etc.). Asterisks (*) indicate the statistical significance of the meta-analysis at FDR <0.05 (see Supplementary Table 10). CA9, carbonic anhydrase IX; CK19, cytokeratin 19; EpCAM, anti-epithelial cell adhesion molecule; ESTIMATE, estimation of stromal and immune cells in malignant tumors; GSEA, gene set enrichment analysis; NES, normalized enrichment score; H&E, hematoxylin and eosin stain; HCC, hepatocellular carcinoma; PI3K/mTOR, phosphoinositide 3-kinase/ mechanistic target of rapamycin; TCGA-LIHC, the cancer genome atlas liver hepatocellular carcinoma; TGFB, transforming growth
factor-β.
Figure 6.Summary of the training/validation cohort. (A) Schematic of tumor blood flow and dedifferentiation of cancer cells. Small tumors with homogeneous blood supply retain liver-specific metabolism and are well-differentiated microtrabecular types and inflamed TIME. As tumors grow and blood flow becomes heterogeneous, tumor markers increase, microtrabecular structures are lost, and biliary stem cell markers turn positive. Poorly differentiated tumors lose normal hepatocyte metabolic functions and exhibit enhanced glycolysis. (B) Summarization of metabolic functions, RNA expression, gene mutations, histopathology, clinical traits, and TIME across five subclasses. AFP, α-fetoprotein; CA9, carbonic anhydrase IX; DCP, des-γ-carboxy prothrombin; HCC, hepatocellular carcinoma; HIF, hypoxia inducible factor; IL6-JAK-STAT3, interleukin-6-janus kinase-signal transducer and activator of transcription 3; TCGA-LIHC, the cancer genome atlas liver hepatocellular carcinoma; TIME, tumor immune microenvironment.
Table 1.The baseline patients’ characteristics
Table 1.
|
Glycolysis subclass |
Rich metabolism subclass |
P-value |
|
Number |
54 |
82 |
|
|
Tumor size (cm) |
5.00 (3.00, 7.00) |
3.50 (2.55, 5.00) |
0.035 |
|
Vascular invasion, yes |
25 (46.3) |
16 (19.5) |
0.002 |
|
NLR |
2.43 (1.69, 3.99) |
1.87 (1.23, 2.94) |
0.007 |
|
CRP |
0.32 (0.11, 3.51) |
0.20 (0.06, 1.14) |
0.059 |
|
AFP |
235.00 (22.25, 3,548.25) |
7.00 (4.00, 38.50) |
<0.001 |
|
DCP |
227.50 (43.50, 4,394.50) |
115.50 (28.50, 546.50) |
0.018 |
|
CT/MRI APHE, yes |
|
|
<0.001 |
|
Inhomogeneous/rim |
34 (70.8) |
24 (32.4) |
|
|
Homogeneous |
14 (29.2) |
50 (67.6) |
|
|
18F-FDG/PET-CT, yes |
|
|
0.02 |
|
High TLR (>2.0) |
5 (100) |
0 (0) |
|
|
Low TLR (<2.0) |
0 (0) |
4 (100) |
|
|
Tumor differentiation |
|
|
<0.001 |
|
Well |
0 (0) |
6 (7.3) |
|
|
Moderate |
31 (57.4) |
72 (87.8) |
|
|
Poorly |
23 (42.6) |
4 (4.9) |
|
|
Microtrabecular type, yes |
18 (36.0) |
60 (81.1) |
<0.001 |
|
Pseudoglandular type, yes |
1 (1.9) |
24 (32.5) |
<0.001 |
|
Compact type, yes |
25 (48.1) |
16 (21.7) |
0.006 |
|
MTM type, yes |
9 (17.6) |
0 (0.0) |
0.001 |
|
W/B mut variant |
|
|
0.013 |
|
APC |
1 (1.9) |
4 (4.9) |
|
|
CTNNB1, N387/T951 |
0 (0.0) |
2 (2.4) |
|
|
CTNNB1, S45/T41 |
1 (1.9) |
11 (13.4) |
|
|
CTNNB1, D32-S37 |
5 (9.3) |
16 (19.5) |
|
|
Hoshida subclass |
|
|
<0.001 |
|
S1/S2 |
53 (98.1) |
14 (17.1) |
|
|
S3 |
1 (1.9) |
68 (82.9) |
|
Table 2.Patients’ characteristics for the five groups
Table 2.
|
NOTCH/TGFβ-High |
PI3K/mTOR-High |
IL6-JAK-STAT3-High |
K-RAS-High |
Wnt/β-catenin-High |
P-value |
|
Number |
26 (19.1) |
28 (20.6) |
34 (25.0) |
21 (15.4) |
27 (19.9) |
|
|
Sex, male |
15 (57.7) |
24 (85.7) |
21 (61.8) |
15 (71.4) |
25 (92.6) |
0.011 |
|
Tumor size |
3.85 (2.50, 6.00) |
5.50 (3.00, 10.35) |
3.00 (2.15, 3.90) |
4.50 (3.00, 6.00) |
4.00 (3.00, 5.15) |
0.005 |
|
Vascular invasion |
8 (30.8) |
17 (60.7) |
5 (14.7) |
3 (14.3) |
8 (29.6) |
0.001 |
|
AFP |
161.50 (30.0, 2,020.5) |
333.0 (18.00, 4,068.3) |
7.00 (3.25, 72.00) |
11.00 (5.00, 19.00) |
4.00 (3.00, 61.50) |
<0.001 |
|
DCP |
141.00 (39.0, 1,395.0) |
729.0 (47.5, 10,714.5) |
80.00 (25.5, 304.0) |
289.5 (34.0, 1,048.5) |
90.50 (28.0, 213.5) |
0.017 |
|
Tumor differentiation |
|
|
|
|
|
<0.001 |
|
Well |
0 (0) |
0 (0) |
3 (8.8) |
1 (4.8) |
0 (0) |
|
|
Moderate |
16 (61.5) |
10 (35.7) |
25 (73.5) |
19 (90.5) |
23 (85.2) |
|
|
Poorly |
10 (38.5) |
18 (64.3) |
6 (17.6) |
1 (4.8) |
4 (14.8) |
|
|
Microtrabecular type, yes |
12 (52.2) |
6 (22.2) |
24 (80.0) |
15 (88.2) |
21 (77.8) |
<0.001 |
|
Pseudoglandular type, yes |
0 (0) |
1 (3.7) |
4 (13.4) |
7 (41.2) |
13 (48.1) |
<0.001 |
|
MTM type, yes |
2 (8.3) |
7 (25.9) |
0 (0) |
0 (0) |
0 (0) |
0.001 |
|
W/B mut, yes |
1 (3.8) |
6 (21.4) |
10 (29.4) |
6 (28.6) |
17 (63.0) |
<0.001 |
|
TP53/cell cycle control mut, yes |
9 (34.6) |
16 (57.1) |
12 (35.3) |
5 (23.8) |
7 (25.9) |
0.091 |
|
Hoshida subclass prediction |
|
|
|
|
|
<0.001 |
|
S1/S2 |
26 (100) |
27 (96.4) |
8 (23.5) |
2 (9.5) |
4 (14.8) |
|
|
S3 |
0 (0) |
1 (3.6) |
26 (76.5) |
19 (90.5) |
23 (85.2) |
|
Abbreviations
adenomatous polyposis coli
aspartate aminotransferase
carbamoyl phosphate synthetase
des-γ-carboxy prothrombin
epithelial-mesenchymal transition
anti-epithelial cell adhesion molecule
gadolinium-ethoxybenzyl-diethylenetriamine
immune checkpoint inhibitor
interleukin-6/janus kinase/signal transducer and activator of transcription
Kyoto Encyclopedia of Genes and Genomes
lymphocyte-activation gene-3
liver imaging reporting and data system
magnetic resonance imaging
mechanistic target of rapamycin
negative for hepatitis B surface antigen and hepatitis C antibody
neutrophil-lymphocyte ratio
phosphoenolpyruvate carboxy kinase
positron emission tomography
programmed cell death protein 1
programmed death-ligand 1
relative enhancement ratio
standardized uptake value
transforming growth factor-β
tumor infiltrating cytotoxic T-lymphocytes
T-cell immunoglobulin mucin-3
tumor immune microenvironment
tumor-to-liver uptake ratios
REFERENCES
- 1. Alexandrov LB, Nik-Zainal S, Wedge DC, Aparicio SA, Behjati S, Biankin AV, et al. Signatures of mutational processes in human cancer. Nature 2013;500:415-421.
- 2. Vogelstein B, Kinzler KW. The path to cancer --Three strikes and you’re out. N Engl J Med 2015;373:1895-1898.
- 3. Dunn GP, Bruce AT, Ikeda H, Old LJ, Schreiber RD. Cancer immunoediting: from immunosurveillance to tumor escape. Nat Immunol 2002;3:991-998.
- 4. Warburg O. On the origin of cancer cells. Science 1956;123:309-314.
- 5. Kumagai S, Koyama S, Itahashi K, Tanegashima T, Lin YT, Togashi Y, et al. Lactic acid promotes PD-1 expression in regulatory T cells in highly glycolytic tumor microenvironments. Cancer Cell 2022;40:201-218.e209.
- 6. Lee H, Choi JY, Joung JG, Joh JW, Kim JM, Hyun SH. Metabolism-associated gene signatures for FDG avidity on PET/CT and prognostic validation in hepatocellular carcinoma. Front Oncol 2022;12:845900.
- 7. Sakamoto M. Pathology of early hepatocellular carcinoma. Hepatol Res 2007;37 Suppl 2:S135-138.
- 8. Matsui O, Gabata T, Kobayashi S, Terayama N, Sanada J, Kouda W, et al. Imaging of multistep human hepatocarcinogenesis. Hepatol Res 2007;37 Suppl 2:S200-205.
- 9. Matsui O, Kobayashi S, Sanada J, Kouda W, Ryu Y, Kozaka K, et al. Hepatocelluar nodules in liver cirrhosis: hemodynamic evaluation (angiography-assisted CT) with special reference to multi-step hepatocarcinogenesis. Abdom Imaging 2011;36:264-272.
- 10. Liver Imaging Reporting and Data System (LI-RADS) v2018 Core algorithm. SITECORE web site, <https://www.acr.org/-/media/ACR/Files/RADS/LI-RADS/LI-RADS-2018-Core.pdf>. Accessed 11 Jun 2024.
- 11. An C, Park S, Chung YE, Kim DY, Kim SS, Kim MJ, et al. Curative resection of single primary hepatic malignancy: liver imaging reporting and data system category LR-M portends a worse prognosis. AJR Am J Roentgenol 2017;209:576-583.
- 12. An C, Kim DW, Park YN, Chung YE, Rhee H, Kim MJ. Single hepatocellular carcinoma: preoperative MR imaging to predict early recurrence after curative resection. Radiology 2015;276:433-443.
- 13. Sakamoto M, Effendi K, Masugi Y. Molecular diagnosis of multistage hepatocarcinogenesis. Jpn J Clin Oncol 2010;40:891-896.
- 14. Wei X, Michelakos T, He Q, Wang X, Chen Y, Kontos F, et al. Association of tumor cell metabolic subtype and immune response with the clinical course of hepatocellular carcinoma. Oncologist 2023;28:e1031-e1042.
- 15. Liu ZZ, Yan LN, Dong CN, Ma N, Yuan MN, Zhou J, et al. Cytochrome P450 family members are associated with fastgrowing hepatocellular carcinoma and patient survival: an integrated analysis of gene expression profiles. Saudi J Gastroenterol 2019;25:167-175.
- 16. Bai J, Tang R, Zhou K, Chang J, Wang H, Zhang Q, et al. An asparagine metabolism-based classification reveals the metabolic and immune heterogeneity of hepatocellular carcinoma. BMC Med Genomics 2022;15:222.
- 17. He Z, Chen Q, He W, Cao J, Yao S, Huang Q, et al. Hepatocellular carcinoma subtypes based on metabolic pathways reveals potential therapeutic targets. Front Oncol 2023;13:1086604.
- 18. Désert R, Rohart F, Canal F, Sicard M, Desille M, Renaud S, et al. Human hepatocellular carcinomas with a periportal phenotype have the lowest potential for early recurrence after curative resection. Hepatology 2017;66:1502-1518.
- 19. Aoki T, Nishida N, Minami Y, Kudo M. The impact of normal hepatobiliary cell zonation programs on the phenotypes and functions of primary liver tumors. Liver Cancer 2025;14:92-103.
- 20. Hoshida Y, Nijman SM, Kobayashi M, Chan JA, Brunet JP, Chiang DY, et al. Integrative transcriptome analysis reveals common molecular subclasses of human hepatocellular carcinoma. Cancer Res 2009;69:7385-7392.
- 21. Hoshida Y, Toffanin S, Lachenmayer A, Villanueva A, Minguez B, Llovet JM. Molecular classification and novel targets in hepatocellular carcinoma: recent advancements. Semin Liver Dis 2010;30:35-51.
- 22. Calderaro J, Couchy G, Imbeaud S, Amaddeo G, Letouzé E, Blanc JF, et al. Histological subtypes of hepatocellular carcinoma are related to gene mutations and molecular tumour classification. J Hepatol 2017;67:727-738.
- 23. Boyault S, Rickman DS, de Reyniès A, Balabaud C, Rebouissou S, Jeannot E, et al. Transcriptome classification of HCC is related to gene alterations and to new therapeutic targets. Hepatology 2007;45:42-52.
- 24. Singal AG, Llovet JM, Yarchoan M, Mehta N, Heimbach JK, Dawson LA, et al. AASLD Practice Guidance on prevention, diagnosis, and treatment of hepatocellular carcinoma. Hepatology 2023;78:1922-1965.
- 25. Cancer Genome Atlas Research Network. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell 2017;169:1327-1341.e1323.
- 26. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612.
- 27. Finotello F, Trajanoski Z. Quantifying tumor-infiltrating immune cells from transcriptomics data. Cancer Immunol Immunother 2018;67:1031-1040.
- 28. Finotello F, Mayer C, Plattner C, Laschober G, Rieder D, Hackl H, et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNAseq data. Genome Med 2019;11:34.
- 29. Torbenson M, McCabe CE, O’Brien DR, Yin J, Bainter T, Tran NH, et al. Morphological heterogeneity in beta-catenin-mutated hepatocellular carcinomas: implications for tumor molecular classification. Hum Pathol 2022;119:15-27.
- 30. Sia D, Jiao Y, Martinez-Quetglas I, Kuchuk O, Villacorta-Martin C, Castro de Moura M, et al. Identification of an immune-specific class of hepatocellular carcinoma, based on molecular features. Gastroenterology 2017;153:812-826.
- 31. Montironi C, Castet F, Haber PK, Pinyol R, Torres-Martin M, Torrens L, et al. Inflamed and non-inflamed classes of HCC: a revised immunogenomic classification. Gut 2023;72:129-140.
- 32. Aoki T, Nishida N, Kurebayashi Y, Sakai K, Morita M, Chishina H, et al. Two distinct characteristics of immune microenvironment in human hepatocellular carcinoma with Wnt/β-catenin mutations. Liver Cancer 2024;13:285-305.
- 33. Halpern KB, Shenhav R, Matcovitch-Natan O, Toth B, Lemze D, Golan M, et al. Single-cell spatial reconstruction reveals global division of labour in the mammalian liver. Nature 2017;542:352-356.
- 34. Moeini A, Sia D, Zhang Z, Camprecios G, Stueck A, Dong H, et al. Mixed hepatocellular cholangiocarcinoma tumors: cholangiolocellular carcinoma is a distinct molecular entity. J Hepatol 2017;66:952-39616.
- 35. Seok JY, Na DC, Woo HG, Roncalli M, Kwon SM, Yoo JE, et al. A fibrous stromal component in hepatocellular carcinoma reveals a cholangiocarcinoma-like gene expression trait and epithelial-mesenchymal transition. Hepatology 2012;55:1776-1786.
- 36. Semenza GL. HIF-1 mediates metabolic responses to intratumoral hypoxia and oncogenic mutations. J Clin Invest 2013;123:3664-3671.
- 37. Niu Y, Lin Z, Wan A, Sun L, Yan S, Liang H, et al. Loss-of-function genetic screening identifies aldolase A as an essential driver for liver cancer cell growth under hypoxia. Hepatology 2021;74:1461-1479.
- 38. Xia H, Chen J, Gao H, Kong SN, Deivasigamani A, Shi M, et al. Hypoxia-induced modulation of glucose transporter expression impacts (18)F-fluorodeoxyglucose PET-CT imaging in hepatocellular carcinoma. Eur J Nucl Med Mol Imaging 2020;47:787-797.
- 39. Hatano E, Ikai I, Higashi T, Teramukai S, Torizuka T, Saga T, et al. Preoperative positron emission tomography with fluorine-18-fluorodeoxyglucose is predictive of prognosis in patients with hepatocellular carcinoma after resection. World J Surg 2006;30:1736-1741.
- 40. Li YC, Yang CS, Zhou WL, Li HS, Han YJ, Wang QS, et al. Low glucose metabolism in hepatocellular carcinoma with GPC3 expression. World J Gastroenterol 2018;24:494-503.
- 41. An J, Oh M, Kim SY, Oh YJ, Oh B, Oh JH, et al. PET-based radiogenomics supports mTOR pathway targeting for hepatocellular carcinoma. Clin Cancer Res 2022;28:1821-1831.
, Naoshi Nishida1
