The normalized level 3 RNA-sequencing (RNA-seq) data and corresponding clinical information of HCC patients were collected from The Cancer Genome Atlas (TCGA) database portal. RNA-seq data (FPKM values) were transformed into transcripts per kilobase million (TPM) values, which are closer in format to those resulting from microarrays and more comparable between samples. Mutational data for all samples with RNA-seq data available was also downloaded. Furthermore, other HCC datasets were obtained from the GEO public database.

Genes belonging to the KEGG subset of canonical pathways were searched in the Molecular Signature Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb) (27), and genes involved in the metabolism pathways were collected. Genes displaying a median TPM < 1 across patients were considered extremely low expression and were removed from the consensus clustering analysis. Consensus clustering was performed on metabolic genes using ConsensusClusterPlus R package (parameters: 160 reps=50, pItem=0.8, pFeature=1). The optimum cluster number was determined by testing 2 to 10 clusters, and based on CDF and Δ(K). The collected pathways and their respective genes were subjected clustering analysis are listed in Supplementary Table S1 .

Hierarchical clustering is used to classify samples into different categories based on the sample’s gene expression matrix, and label different sample categories (28). Taking the expression value of the gene as an input feature and the classification of the sample as an output value, a random forest algorithm was used to construct a prediction model. We used the out-of-bag error rate as an index to evaluate the importance of the features and rank them, and take the first 100 and the first 50 features as the features for subsequent analysis. Scikit-learn was used as a tool for related calculations (29).

Level-3 RNA-sequencing (RNA-seq) data containing 374 HCC samples and 50 normal controls was downloaded from The Cancer Genome Atlas (TCGA). The differentially expressed genes (DEGs) were identified using the Edge R package for R software. |Fold change| > 2 and adjusted P-values < 0.05 were set as the thresholds.

The Gene set enrichment analysis (GSEA) was performed to identify enrichment degree of gene set among tumors of a certain metabolic subtype and compared with other samles (27). Significance was considered for values of corrected P as recommended by the software (https://www.gsea-msigdb.org/gsea/downloads.jsp). Gene sets were downloaded from the MSigDB. The GSEA results were merged by using GSVA R package.

Gene Set Variation Analysis (GSVA) is a pathway/gene set-based analysis approach that provides an overall pathway or gene-set activity score for each sample (30). We estimated the activities of metabolic pathways through GSVA. GSVA z-scores for KEGG metabolic gene-sets from the MSigDB and other gene-sets were calculated using the GSVA R package. We considered that the activity a given metabolic pathway was different between two groups if the median of GSVA values across one group differed significantly (Mann-Whitney q-value<0.1) and more than 0.2 from that measured across another group. Besides, we estimated T cell exhaustion scores for each HCC sample by using GSVA and a gene set including signature genes of exhausted CD8+ T cells (31).

GSVA and the single sample Gene Set Enrichment Analysis (ssGSEA) are two most widely used methods that carry out sample level enrichment analysis. Both are unsupervised gene set enrichment methods that compute an enrichment score integrating the collective expression of a given gene set relative to the other genes in the sample. A previous study revealed a significant consistency of the two methods in evaluating immune populations across 28 cancer types (32). Besides, they also compared different signatures identifying immune cell populations reported in previous articles and constructed a new set of gene signatures with better discrimination than that shown in previous sets. According to their experience, we chose GSVA method because its advantage in reducing the noise of the data; the immune signatures were also retrieved from the same publication (32).

Tumor Immune Single Cell Hub (TISCH, http://tisch.comp-genomics.org), a single cell RNA sequencing (scRNA-seq) application platform that can systematically and comprehensively study the heterogeneity of tumor microenvironment (TME) (33), was used to explore the relationship between PPT1 and the TME on cell level.

The human HCC cell lines Huh-7 (JCRB0403) and JHH-7 (JCRB1031) were obtained from the Health Science Research Resources Bank (Osaka, Japan). Cells were grown in Dulbecco’s Modified Essential Media (DMEM, Hyclone) supplemented with 10% fetal bovine serum (FBS) and penicillin/streptomycin (Gibco, Grand Island, NY, USA) in a humidified incubator at 37°C in a 5% CO₂ atmosphere.

Human pyruvate kinase M2 (PKM2) cDNA was cloned into the pCDH-CMV-MCS-EF1-Neo vector and palmitoyl-protein thioesterase 1 (PPT1) was cloned into the pCDH-CMV-MCS-EF1-Puro vector by IGE (Guangzhou, China). After confirming the PKM2 and PPT1 sequence by sequencing, the plasmid was co-transfected into HEK293T cells (ATCC, RRID: CVCL_0063) with the lentivirus packaging plasmids psPAX2 and pMD2G to produce lentivirus particles. To knock down the PKM2 and PPT1 expression, the shRNA sequence of PKM2 was cloned into the pLKO.1-Neo lentiviral vector and the shRNA sequence of PPT1 was cloned into the pLKO.1-Puro lentiviral vector by IGE (Guangzhou, China). The oligonucleotide sequences of shRNAs are listed in Supplementary Table S5 . Similarly, the lentiviral vector was transfected into HEK293T cells together with packaging vectors pMD2.G and psPAX2 to produce lentivirus particles. The media containing retroviruses were collected 72h after transfection, centrifuged to remove cell debris, and then filtered. The lentivirus particles were then used to transduce Huh-7 and JHH-7 cells. Control cells were transfected with a control empty vector.

Total protein was extracted by using RIPA lysis buffer containing a cocktail of protease and phosphatase inhibitors. Using the BCA Protein Assay Kit (CWBio, China, Cat#CW0014S), the concentration of the protein was determined. Following the previous process, 20 μg total proteins extracted from Huh-7 cells were separated electrophoretically in 10% SDS polyacrylamide gels. Separated proteins are transferred to polyvinylidene fluoride (PVDF) membranes. The membrane was then blocked with 5% BSA at room temperature for 1h before being incubated with the primary antibodies overnight at 4°C. Next day, the membranes were incubated for 1h at room temperature with secondary antibodies. Eventually, the immunoblotting was tested using the ECL Chemiluminescence Kit (Thermo Fisher Scientific, USA, Cat#32109). The antibodies are shown in Supplementary Table S4 .

MTT assay was used to evaluate cell viability. Cells were seeded in 96-well plates (2 × 10³/well) (in triplicate for each condition) and were cultured for different hours (34). The 20 µl MTT (Sigma, Saint Louis, USA) stock solution (5 mg/ml) was added to each well and incubated for 2h, the absorbance was measured at a wavelength of 570 nm.

ECAR was estimated by using Seahorse assays with a Seahorse XFp Analyzer, 2.0 × 10⁴ Huh-7 cells per well were plated onto 8-well plates and glycolysis stress test was performed following the manufacturer’s specifications (35, 36). The assay DMEM media was free of glucose and pyruvate. The concentration of drugs used for test were as following: glucose (10.0 µM), oligomycin A (1.0 µM), and 2-deoxy-D-glucose (2-DG, 100 µM).

5 x 10⁵ HCC cells were incubated in phenol red-free medium at 37°C, 5% CO₂ for 1h. The supernatant medium was collected. The lactate secreted by the cells was measured using the Lactate Colorimetric/Fluorometric Assay Kit according to the manufacturer’s protocol (BioVision, Inc.) (37).

All quantitative data were conducted at least three independent experiments and expressed as the mean ± SD. For comparisons between two groups, statistical significance for normal distribution data was estimated by unpaired Student t-test, otherwise Mann-Whitney U test was applied. For comparisons of more than two groups, One-way ANOVA multiple comparison was used. Correlation coefficients were computed by Spearman and distance correlation analyses. The influence of a single factor on the survival was evaluated through the Cox proportional hazard model, Kaplan-Meier survival curves, or Log-rank test. All statistical analyses were conducted using R (https://www.r-project.org/) or SPSS software (version 17.0). Two-sided P values of less than 0.05 were considered statistically significant.

We obtained four robust subtypes of HCC in TCGA dataset ( Supplementary Table S2 ). A heatmap based on one hundred genes obtained with machine learning was shown in Figure 1A . The robust classification revealed that HCC tumors displayed highly heterogeneous expression of genes that were directly involved in metabolism. The mHCC1 was defined as the largest group (144/371; 38.8%), followed by mHCC2 (138/371; 37.2%), mHCC4 (51/371; 13.7%), and mHCC3 (38/371; 10.2%). These subtypes are henceforth termed mHCC1 (cluster 1), mHCC2 (cluster 2), mHCC3 (cluster 3), and mHCC4 (cluster 4). We noted distinct patterns of gene mutations across the four different metabolic subtypes ( Figures 1B, C ). For example, mutations in TP53 were enriched in mHCC1, and 76% of MHCC4 tumors showed mutations in CTNNB1. In addition, co-occurrence analysis using the five most frequent mutations in HCC indicated different relationships among the five common mutations across four subtypes. A mutual exclusion between TP53 mutations and CTNNB1 mutations in mHCC1 tumors was observed, while co-occurrence between TTN mutations and mutations in CTNNB1 or ALB were found in mHCC2. These results indicated that there are multiple metabolic phenotypes in HCC driven by different oncogenic mutations. The mHCC1 was associated with the shortest survival while the prognosis of other three subtypes were similar ( Figures 1D, E ).

To identify the underlying biological functional characteristics of each subtype, signature genes in each group were collected. KEGG analysis of the signature genes was conducted and the results were visualized in Figure 2 . Among those signature genes dysregulated in each subtype/cluster, we identified that the mHCC1 had much larger proportion of genes enriched in metabolism-related pathways than that in other clusters; especially including “carbohydrate metabolism”, “lipid metabolism”, and “amino acid metabolism”, which were major reprogrammed metabolic pathways in malignant cells. In addition, we noted that most genes that were enriched in “signal transduction”, “signaling molecules and interaction”, and “cellular community” came from mHCC4, followed by mHCC2 and mHCC3, but barely from mHCC1. Besides, immune-related pathways were more likely to be enriched by genes from mHCC3 and mHCC4, such as “Immune system”, “Infectious disease: bacterial”, “Infectious disease: viral”, “Infectious disease: parasite”, and “Immune disease”.

We next supposed that subtypes produced through gene clustering had specific features in different metabolic pathways and employed a sample-level gene set enrichment method (GSVA) to compute the GSVA enrichment scores of the selected metabolic pathways. The differences in GSVA scores for each pathway between normal liver and four subtypes were calculated. As shown in Figure 3A , compared with other subtypes, mHCC1 displayed significantly dysregulated fold change in GSVA scores with reference to those in normal livers. It was well known that HCC cells are metabolically distinct from normal hepatocytes and express different metabolic enzymes (38, 39). Our next gene set enrichment analysis (GSEA) using a gene set involving metabolic genes of normal liver revealed a significant absence of normal metabolic genes in samples of mHCC1 and mHCC4, but no significant differences in mHCC2 and mHCC3 ( Figure 3B ). These results were consistent with the KEGG analysis visualized in Figure 2 and suggested different types of metabolism impairments among subtypes. Similarly, as shown in Figures 3C, D , there was a remarkable difference between GSVA scores for most metabolic pathways of the mHCC1 and other three clusters as well as normal liver samples; in contrast, GSVA scores of the most pathways in both mHCC2 and mHCC3 showed limited differences compared with normal liver samples, which is consistent with the GSEA results.

We next compared GSVA scores of major metabolism reactions that are critical to carcinogenesis, including carbohydrate metabolism ( Figure 3E ), amino acid metabolism ( Figure 3F ), and lipid metabolism ( Figure 3G ). Overall, compared with normal liver samples, tumors in mHCC1 were most significantly differentially enriched for almost all the selected pathways. The GSVA enrichment pattern in different pathways of carbohydrate metabolism was identical, mHCC1 had most significant differences in GSVA scores than in normal liver, followed by mHCC4, while mHCC2 and mHCC3 had similar scores compared with each other. In the respect of amino acid metabolism and lipid metabolism, mHCC2, mHCC3, and mHCC4 had similar enrichment in many pathways but mHCC1 still showed a significantly differential enrichment for all pathways. Specifically, mHCC4 displayed much more significant enrichment in “Pentose phosphate pathway” than other clusters and showed a similar enrichment in both “Alanine aspartate and glutamate pathway” and “Glycolipid metabolism pathway” compared with mHCC1. In addition to the characteristic metabolic pathways specifically present above, there are some amino acid metabolic pathways that are significantly reduced in mHCC1.

To complete the landscape of the metabolic pathways across the four subtypes, we further evaluated whether each subtype had unique metabolic pathways enriched. If a pathway was especially different (Mann-Whitney q-value<0.05 and difference between medians of GSVA was more than 0.2) in only one cluster, it was regarded as a characteristic metabolic pathway in corresponding subtype. As shown in Figures 4A–C , there were six characteristic pathways in total identified in mHCC1, mHCC3, and mHCC4. “Lysine Degradation” and “On-Carbon-Pool by Folate” were identified as characteristic metabolic pathways of mHCC1, “Oxidative Phosphorylation” and “O-Glycan Biosynthesis” in mHCC3, and two sub-pathways of “Glycosphingolipid Biosynthesis” in mHCC4, while none in mHCC2. In addition to the characteristic metabolic pathways specifically present above, there are some amino acid metabolic pathways that are significantly reduced in mHCC1.

The above observations led us to question how these differences among metabolic subtypes relate to clinical outcome of the patients. Therefore, the prognostic value of metabolic genes in each subtype was investigated using log-rank test and univariate Cox regression model. The results of survival analysis are summarized in Supplementary Table S3 . Remarkably, only a small portion of metabolic genes with prognostic significance were shared by multi-subtypes, and none was shared by all four subtypes ( Figures 5A, B ), suggesting distinct role of genes across subtypes even for the same metabolic pathway. Furthermore, survival analysis based on genes involved in the six characteristic metabolic pathways also revealed specific prognostic value of metabolic genes in their corresponding clusters ( Figures 5C, D ). These data suggested subtype-specific roles for these unique tumor metabolic pathways as mechanisms contributing to tumor progression and identified some genes of unique metabolic pathways as potential targetable metabolic vulnerabilities in subtype-specific manner.

The immune response is associated with dramatic modifications in tissue metabolism, and emerging evidences suggested the tremendous impact of tumor metabolic reprogramming on the immune response (40). Therefore, we assessed the relationship between the metabolic subtypes and previously defined scores which indicating the abundance of various immune cell populations. The immune landscape revealed significant intra-cluster heterogeneity ( Figure 6A ). Consistent with an immune-cold phenotype, tumors in mHCC3 had the lowest rate of leukocytes infiltration. Contrary to the immune-cold mHCC3, tumors in mHCC1 had the highest cell infiltration fraction and it was mainly characterized by higher infiltration of lymphocytes. Interestingly, although tumors in mHCC1 displayed the highest infiltration of cytotoxic cells including CD8+ T cells, Tgd, and NK cells, mHCC1 was the subtype with both poorest prognosis and highest T cell exhaustion score ( Figure 6B ), suggesting that the tumor microenvironment of mHCC1 was immune-hot but highly immune-exhausted type. Most suppressive ligands and receptors for immune checkpoint showed the highest expression levels in mHCC1 ( Figure 6C ), which is consistent with an immune-suppressive phenotype. Emerging experimental data indicated that the presence of a pre-existing intra-tumoral T cell infiltration, checkpoint molecules (PD-1, PD-L1 expression) could favor a clinical response (41), therefore, we hypothesized that this subtype could benefit much from immune therapy only if the immunosuppressive situation was relieved.

Different from the extreme low/high immune infiltration shown above for mHCC3 and mHCC1, mHCC2 demonstrated a more balanced and favorable immune profile. One major difference from mHCC1 is that the immune composition in mHCC2 was enriched with B cells, macrophages, mast cells, and neutrophils, suggesting an immune microenvironment tended to towards innate immunity and inflammation. The remaining mHCC4 were more diverse with intermediate levels of immune features. The mHCC4 had high infiltration levels of NKdim, Th, Tcm, and Tgd, but showed a low level of other lymphocytic cells, and extremely low levels of B cells, macrophages, mast cells, and neutrophils, which was opposite to mHCC2. Additionally, mHCC4 had the highest level of eosinophils. Together with the specificity of eosinophils infiltration as an index for good prognosis in mHCC4, the above findings supported the major anti-tumor role of eosinophils in the metabolic subtype mHCC4.

The prognostic role of immune cells shown by survival analysis between inter-and intra-cluster was in a high degree of heterogeneity ( Figure 6D ). For example, among the important cytotoxic T cells, such as CD8+ T cells, and Tgd cells, we found that increased rate of CD8+ T cell infiltration was associated with good prognosis in mHCC1, mHCC2, and mHCC3, but it was not linked with prognosis in mHCC4; in contrast, Tgd infiltration was associated with significantly poor prognosis in mHCC4, but showed no correlation with prognosis in other three clusters. Meanwhile, some immune cells had the opposite prognostic role in different clusters, such as that infiltration of Treg cells was a poor factor in the prognosis of mHCC4 but a good factor of prognosis in mHCC3. The prognostic heterogeneity of immune infiltration revealed different roles of immune cells in different HCC subtypes. Since the metabolic subtypes were produced based on metabolism genes that did not involve signature genes of immune cells, these results suggest a high heterogeneity in component and function of immune microenvironment across these four metabolic subtypes, and support a combination treatment strategy targeting metabolism and immune microenvironment, or the attempt to dig out new target of dual function that combined metabolism-regulation as well as immune-modulation.

Considering that mHCC1 had the poorest prognosis, we tried to identify new therapeutic targets for this subtype. The prognostic metabolic genes were identified, and their correlation with T cell exhaustion, the important feature of mHCC1, was investigated ( Figure 7A ). Of genes most closely related with T cell exhaustion, we identified palmitoyl-protein thioesterase 1 (PPT1), which was specifically upregulated in mHCC1 and was a specific prognosis factor for mHCC1 ( Figures 7B, C ). The positive correlation between PPT1 and Tex score ( Figure 7D , R = 0.3495, P<0.0001) suggested that PPT1 might contribute to T cell exhaustion. To determine the relationship between PPT1 and glucose metabolism, we analyzed the data of Transcriptomics and metabolomics in the study of Chaisaingmongkol et al. (42). Data LIHC_GSE76297 was used to confirm the high expression of PPT1 in HCC tumor tissues. Through the analysis of glucose metabolism-related metabolites, the content of glucose, maltose and maltotriose decreased in the high expression group of PPT1, suggested that the level of glycolysis was exuberant in tissues with high expression of PPT1 ( Figures 7E, F ). Then, we analyzed the data LIHC_GSE166635 in GEO database using TISCH platform, and the results showed that the expression of PPT1 in malignant cell was strikingly higher than that in epithelial cells ( Figures 7G, H ). To verify the biological function and mechanism of PPT1, functional annotation analysis was performed on PPT1 in malignant cell. The KEGG pathways and GO biological processes enrichment suggested that PPT1 was closely related to tumor metabolism, immunity, tumor microenvironment and extracellular matrix ( Figures 7I, J ), which was consistent with the above results that PPT1 is related to tumor metabolism and immunosuppression (43, 44).

In a very recent study, Rebecca et al. identified PPT1 is the molecular target of a novel chloroquine derivative, the DC661; they found that knockout of PPT1 in several cancer cell lines using CRISPR-Cas9 editing resulted in significant impairment of tumor growth similar to that observed with DC661 treatment, which supported the tumor driver role of PPT1 in cancer (45, 46). Therefore, a further question was raised whether PPT1 as a therapeutic target had subtype-specific property and immunomodulating function. In vitro study was carried out to verify the above questions. We enhanced the glycolysis activity through overexpression of PKM2 in Huh-7 and JHH-7 cells, cell lines with the intermediate glycolysis level. Downregulation of PPT1 significantly neutralized the proliferative advantage of PKM2 overexpression cells ( Figures 8A–C ). In addition, PKM2 overexpression markedly increased ECAR ( Figure 8D ), and we found downregulation of PPT1 in PKM2-overexpressed cells induced a larger drop in ECAR than it did in normal cells, and the same was observed for the lactate production ability ( Figure 8E ), which partly explained the positive correlation between PPT1 expression and Tex score. We next divided TCGA samples into two groups based on the expression of PKM2, HK2, or FBP1, respectively. Interestingly, the prognostic value of PPT1 expression was only significant in TCGA samples with higher expression of PKM2 ( Figure 8F ) or HK2 ( Figure 8G ), two enzymes promote glycolysis, or in TCGA samples with lower expression of FBP1 ( Figure 8H ), an enzyme inhibits glycolysis.

Recently, it has been reported that PPT1 is essential for the function of lysosomes, which, as the center of cellular energy sensing and metabolic regulation, play a driving role in the malignant progression of tumor cells (45). In HCC, the increase of glycolysis level and the further enhancement of lysosome and autophagy activity lead to poor prognosis of patients (47). PKM2 is a key rate-limiting enzyme in glycolysis and an important factor in HCC metabolism (48). Therefore, we decided to find out whether PPT1 depends on PKM2 to promote tumor progression through autophagy. Chloroquine, a lysosomal autophagy inhibitor, was used to inhibit autophagy. In Colony formation assay ( Figures 9A–D and Supplementary Figures S1A–D ), Transwell assay ( Figures 9E–H and Supplementary Figures S1E–H ) and EdU assay ( Figures 9I–L and Supplementary Figures S1I–L ), compared with the stably knockdown PKM2 cells, PPT1 could significantly promote the proliferation, migration and invasion of HCC tumor cells in stable overexpression of PKM2 between Huh-7 and JHH-7 cells, interestingly, the proliferation, migration and invasion of tumor cells decreased strikingly after the treatment of chloroquine. Furthermore, Western blotting analysis showed that under the condition of high expression of PKM2, PPT1 could significantly affect the expression of autophagosome protein marker LC3B, while the accumulation of LC3B -II expression increased after chloroquine was added, suggesting that the lysosomal function might be impaired ( Figures 9M, N ). Collectively, PPT1 may be dependent on PKM2 to promote the progress of HCC through the autophagy.

To further verify the role of PPT1, we constructed Mouse subcutaneous xenograft models. Compared with the control group, the injection of Huh-7 cells stably overexpressing PPT1 and PKM2 significantly promoted tumor growth, while after treatment with chloroquine, tumor growth was inhibited to some extent ( Figures 10A–D ). In addition, IHC staining showed that Ki-67 levels were notably increased in the tissues of mice treated with stable overexpression of PPT1 and PKM2, suggesting that tumor proliferation was obvious ( Figures 10E, F ). Taken together, PPT1 relies on PKM2 to promote the malignant progression of HCC by mediating autophagy.