Single-Cell RNA Sequencing Revealed That the Enrichment of TPI1⁺ Malignant Hepatocytes Was Linked to HCC Metastasis and Immunosuppressive Microenvironment

Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignancy and primary liver cancer, with an estimated incidence of more than 1 million by 2025.^1,2 Despite the adoption of promising treatment strategies, including curative resection, transarterial chemoembolization (TACE), molecular targeted therapy, and immunotherapy, the prognosis of HCC remains unsatisfactory.^3,4 The main reasons for poor prognosis in HCC include its propensity for metastasis and recurrence, coupled with the lack of effective predictive and therapeutic strategies.^5,6 Exploring possible factors affecting HCC metastasis and potential treatment options will effectively improve the prognosis of the patients.^5,7 As a typical inflammation-related process in tumorigenesis, immune evasion is an important feature of HCC occurrence and development.⁸ Moreover, the tumor microenvironment (TME) in HCC has a strong dependence on the number and status of immune cells, which leads to poor clinical effectiveness of HCC treatment.⁹ Therefore, elucidating the role of metastasis-related genes in the immune microenvironment is crucial for the prognosis and treatment strategies of HCC.

Due to the highly heterogeneous nature of HCC, its current efficacy is limited.¹⁰ Single-cell RNA sequencing (scRNA-seq) technology could observe the unique characteristics of individual cells and gain insight into the tumor biology, thereby resolving the expression patterns of unique cell populations and enabling studies of the relationship between unique subtypes and HCC.^11,12 The unique ecosystem of HCC tissue contains a variety of cell lines, including hepatocytes and immune cells. The formation of tumor requires mutual regulatory functions and crosstalk between specific cell types.¹³ A detailed understanding of specific cell clusters in the HCC microenvironment could facilitate the identification of novel biomarkers and improve the prognosis of HCC patients.

In the present study, we investigated the biological functions of malignant hepatocyte clusters and metastasis-related genes using scRNA-seq data retrieved from the Gene Expression Omnibus (GEO) database. The correlations between key gene expression and HCC prognosis, immune cell infiltration, and immune checkpoints were also analyzed using public databases and validated with clinical samples. Assessing the role of metastasis-related gene in HCC microenvironment may offer novel insights for the surveillance and treatment of HCC.

Materials and Methods

Data Collection and Tissue Samples

The single-cell RNA-sequencing data, including 10 HCC tissue samples, were downloaded from the GEO database GSE149614. Using the Tumor-Node-Metastasis (TNM) staging system, 6 patients were divided into the non-metastasis group (TNM stage I–IIIA), and 4 patients were divided into the metastasis group (TNM stage IIIB-IV). In addition, the transcriptome data and corresponding clinicopathological and survival information of the HCC patients were also collected from The Cancer Genome Atlas liver hepatocellular carcinoma (TCGA LIHC) database. Furthermore, 46 HCC tissues and matched tumor-adjacent and normal tissue samples were obtained by surgical resection at Sichuan Cancer Hospital between June 2018 and July 2021. The follow-up ended in June 2023. Postoperative pathological results of all tissue samples confirmed HCC.

The study was approved by the Clinical Research Ethics Committee of Sichuan Cancer Hospital, and written informed consent was received from all participants. All procedures were carried out in accordance with the ethical guidelines of the Helsinki Declaration.

Identification of Hepatocyte Subsets and Differential Expressed Metastasis-Related Genes

To determine the marker genes for each malignant hepatocyte cluster, the “FindAllMarkers” function in the R package “Seurat” was used to compare cells of one cluster with all other cells. In addition, different malignant hepatocyte clusters were annotated by the R package “SingleR” based on the composition pattern of the marker genes. And the transcriptome dataset Human Primary Cell Atlas was referenced to independently infer the origin cells of each single cell and identify the cell type. Hepatocytes were characterized by high expression of ALB, APOC3, APOB, APOA1, HPR, SEC16B, BCHE, etc. Moreover, difference analysis was conducted using HCC tissues with and without metastasis. The Wilcox test was performed to determine differentially expressed genes (DEGs) according to the criteria of P value < 0.05 and log₂ fold change > 0.58. All codes used in this study were available in the R package “Seurat” and “SingleR”.

Go and KEGG Analysis

To further explore the potential biological mechanisms, 9 malignant hepatocyte clusters and 465 differential expressed metastasis-related genes were analyzed through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, respectively.

Prognostic Model Development for HCC

Based on the TCGA-LIHC database, univariate Cox regression analysis was utilized to determine the correlation between metastasis-related genes and overall survival (OS). With p < 0.05, prognosis-related genes were identified. Based on these genes, the Least Absolute Shrinkage and Selection Operator (LASSO)-penalized Cox regression analysis were then conducted to construct a prognostic model and reduce the complexity of solving overfitting problems.

Immunohistochemistry (IHC) and Immunofluorescent (IF) Staining

IHC staining was performed on the tumor tissues derived from 46 HCC patients in our center. First, the tumor tissues were fixed in 10% formalin for 48 hours and embedded in paraffin to prepare the specimens. Then, the specimens were incubated with 1:150 rabbit monoclonal to TPI1 (Abcam, USA) overnight at 4°C. All experiments were performed three times. The intensity of immunohistochemical staining of TPI1 was evaluated by calculating the percentage of positive cells and the staining concentration in 5 random fields of views. All IHC results were assessed by two independent pathologists.

For IF staining of TPI1 and CD8, the cells were fixed in 4% paraformaldehyde for 25 minutes, permeated in 0.5% Triton X-100 for 20 minutes, sealed in 5% BSA for 1 hour, and then incubated overnight with the primary antibodies. After washing the cells in poly (butylene succinate-co-butylene terephthalate) (PBST) three times, the cells were incubated with the secondary antibody for 1 hour. Then, the nucleus was re-stained with 4’,6-diamidino-2-phenylindole (DAPI). Finally, the staining results were observed by the confocal laser scanning microscope (Olympus Corporation, Tokyo, Japan). All experiments were performed three times.

Statistical Analysis

All statistical analyses were performed using GraphPad Prism version 8.0 and R software version 4.0.1. Comparisons between the two groups were evaluated using the Student’s t–test or the Wilcoxon rank-sum test, as appropriate. Differences with p < 0.05 were considered statistically significant. Univariate analysis was conducted by the Cox proportional hazards model. Survival curves were plotted by Kaplan-Meier method and the differences were compared by log rank test.

Results

Hepatocyte Clustering and Subtype Analysis

To reveal tumor heterogeneity in HCC, we explored the malignant hepatocyte clusters in HCC via analyzing single-cell RNA-seq data of 10 HCC tissues obtained from GEO database GSE149614. A total of 14,779 single cells were analyzed to obtain these scRNA-seq data. Subcluster analysis identified 9 malignant hepatocyte clusters with different marker genes and biological functions (Figure 1A and B). In GO enrichment analysis and KEGG analysis, the top 10 GO and KEGG pathways revealed that C1_Hepatocyte-CRP, C2_Hepatocyte-TFDP3, C4_Hepatocyte-SPP1, C5_Hepatocyte-HHIPL2, C7_Hepatocyte- GTSE1, and C9_Hepatocyte-CT45A10 were mainly involved in the activation of oncogenes. C8_Hepatocyte-SLCO1B3 was mainly involved in regulating the metabolic processes. C3_Hepatocyte-SERF2 and C6_Hepatocyte-IL13RA2 were mainly involved in regulating the immune microenvironment via regulating the JAK/STAT signaling pathway and primary immunodeficiency, respectively (Figure 1C and D). The top three marker genes of these malignant hepatocyte clusters were as follows: CRP, FGF19, and SLC25A47 (C1_Hepatocyte-CRP); TFDP3, XX-CR54.1, GREM2 (C2_Hepatocyte-TFDP3); SERF2, TMSB4X, RPS19 (C3_Hepatocyte-SERF2); SPP1, ADH4, ALDH3A1 (C4_Hepatocyte-SPP1); HHIPL2, SLC22A11, NTM (C5_Hepatocyte-HHIPL2); IL13RA2, C11orf53, C6orf141 (C6_Hepatocyte-IL13RA2); GTSE1, DLGAP5, BUB1 (C7_Hepatocyte- GTSE1); SLCO1B3, TAT-AS1, LINC01554 (C8_Hepatocyte-SLCO1B3); CT45A10, ADAD1, MSC-AS1 (C9_Hepatocyte-CT45A10) (Figure 2A). In conclusion, we identified 9 malignant hepatocyte clusters with different marker genes and biological functions, uncovering their relationship with the immune microenvironment.

Figure 1 Heterogeneity of hepatocytes in HCC tissues. (A) malignant hepatocyte clusters were identified; (B) 9 major malignant hepatocyte clusters were further determined; (C) The top 10 GO enrichment pathways of 9 malignant hepatocyte clusters; (D) The top 10 KEGG pathways of 9 malignant hepatocyte clusters. Abbreviations: HCC, hepatocellular carcinoma; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 2 Landscape of differentially expressed genes associated with HCC metastasis. (A) Heatmap of top 10 DEGs between the metastasis and non-metastasis group in 9 malignant hepatocyte clusters; (B and C) Go enrichment analysis of total DEGs between the metastasis and non-metastasis group; (D and E) KEGG pathway analysis of total DEGs between the metastasis and non-metastasis group. Abbreviations: HCC, hepatocellular carcinoma; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of Differentially Expressed Genes Associated with HCC Metastasis

Among the 10 HCC patients analyzed in single-cell sequencing data, 6 patients were at the early/intermediate stage without HCC metastasis. Therefore, they were divided into non-metastasis group. Specifically, three of them were at TNM stage I, one patient was at TNM stage II, and 2 of them were at TNM stage IIIA. Four patients were at the advanced stage with HCC metastasis. And they were divided into metastasis group. Specifically, 2 of them were at TNM stage IIIB and the other 2 patients were at TNM stage IV. In C3_Hepatocyte-SERF2, C6_Hepatocyte-IL13RA2, and C9_Hepatocyte-CT45A10, the gene expression in metastasis group was higher than that in non-metastasis group (Supplemental Table 1). And in C1_Hepatocyte-CRP, C5_Hepatocyte-HHIPL2, and C8_Hepatocyte-SLCO1B3, the gene expression in metastasis group was lower than that in non-metastasis group (Supplemental Table 1).

In 9 malignant hepatocyte clusters, four hundred and sixty-five differentially expressed genes (DEGs) between metastasis and non-metastasis group were identified (Supplementary Table 2). To explore the biological functions of the DEGs, we further conducted GO enrichment analysis and KEGG analysis on these genes. In the biological processes (BP) of GO enrichment analysis, the DEGs mainly participated in immune system process and metabolic process (Figure 2B and C). According to KEGG analysis, these genes were mainly involved in cancers and immune system pathways (Figure 2D and E). Therefore, we believed that the immune microenvironment was a significant pathway involved in regulating HCC metastasis.

TPI1 Was Verified to Be the Key Gene of C6_Hepatocyte-IL13RA2 and HCC Metastasis

In order to verify the DEGs associated with metastasis, further bioinformatics and IHC analysis were performed. Based on the TCGA-LIHC database, univariate Cox regression analysis was performed to identify metastasis-related genes associated with the patient prognosis. The results showed that a total of 79 metastasis-related genes were significantly associated with the OS of HCC patients (Supplemental Table 3, P < 0.05). LOSSO regression analysis was further conducted to construct a prognostic model with these genes (Figure 3A–F, Supplemental Table 4). Finally, four key metastasis-related genes (ITM2A, HSP90AA1, RSU1, and TPI1) were selected (Coefficient >0.1). We found that these four genes were all marker genes of C6_Hepatocyte-IL13RA2 and preliminarily validated their expression in HCC tissues with and without the evidence of metastasis (Supplemental Figure 1). The results showed that only the expression of TPI1 presented an association with the prognosis and metastasis of HCC. Therefore, the mRNA level and prognosis results of TPI1 were further analyzed in the TCGA-LIHC database. According to the results, the expression of TPI1 was much higher in HCC tissues than in normal controls based on patient age, gender, pathologic stage, pathologic T stage, histologic grade and vascular invasion (Figure 4A). Consistent with the previous results, the expression of TPI1 was significantly higher in patients with advanced pathologic T stage (Figure 4A). The prognosis results also demonstrated that patients with higher TPI1 expression presented poorer overall survival (P < 0.001) and disease-specific survival (P = 0.010) (Figure 4B). The results suggested that TPI1, a key gene in C6_Hepatocyte-IL13RA2, was significantly associated with the prognosis and metastasis of HCC.

Figure 3 LOSSO regression analysis of DEGs associated with HCC metastasis and prognosis in the TCGA database. (A) The results of LASSO regression for the 15 DEGs; (B) Risk score of each HCC patient; (C) Patient survival based on the risk score; (D) Heatmap of the 15 DEGs in the low- and high-risk group; (E) Kaplan-Meier survival curves of OS for the patients in the low- and high-risk group; (F) ROC curves for the prognostic performance of the model for 1-, 3-, and 5-year OS. Abbreviations: LASSO, Least Absolute Shrinkage and Selection Operator; DEGs, differentially expressed genes; HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; OS, overall survival; ROC, receiver operating characteristic.

Figure 4 The correlation between TPI1 expression and clinicopathological data, prognosis, and immune cell infiltration in HCC patients. (A) The correlation between TPI1 expression and the clinicopathological data in the TCGA database; (B) The Kaplan-Meier survival curves of OS and DSS for the HCC patients with different TPI1 expression in the TCGA database; (C) Relationship between TPI1 expression and immune cell infiltration in the TCGA database; (D) Relationship between TPI1 expression and immune cell infiltration in the TIMER database; (E) Relationship between TPI1 expression and the expression of PDCD1 and CTLA4 in the TIMER database. *P < 0.05; **P < 0.01; ***P < 0.001. Abbreviations: HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; OS, overall survival; DSS, disease specific survival; TIMER, Tumor Immune Estimation Resource.

The Role of TPI1 in Immune Microenvironment

According to the previous results, TPI1 was one of the marker genes of C6_Hepatocyte-IL13RA2, which was mainly involved in regulating the immune microenvironment of HCC (Figure 1C and D). The immune microenvironment was also a main pathway for the enrichment of metastasis-related genes (Figure 2B-E). Therefore, the relationships between TPI1 expression and immune cell infiltration (B cells, Th2 cells, CD8⁺ T cells, dendritic cells (DCs), neutrophils, and macrophages) were analyzed based on the TCGA-LIHC database. As shown in Figure 4C, the level of TPI1 expression in HCC was negatively correlated with the levels of infiltrating CD8⁺ T cells (R = −0.199, P < 0.001), DCs (R= −0.172, P < 0.001), and B cells (R = −0.139, P = 0.007), while it was positively correlated with the levels of infiltrating Th2 cells (R= 0.356, P < 0.001). Moreover, the expression of TPI1 was positively correlated with the levels of infiltrating regulatory T cells (Tregs) and myeloid-derived suppressor cells (MDSCs) in the Tumor Immune Estimation Resource (TIMER) database (R= 0.265, P < 0.001; R= 0.460, P < 0.001; Figure 4D).

Programmed cell death protein 1 (PD-1) and cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), two prominent immune checkpoints, have been found to be the most reliable targets for attenuating anti-tumor immune responses in numerous malignancies, including HCC.^14,15 They play a pivotal role in the regulation of immune responses within the tumor microenvironment. To preliminarily explore the potential mechanisms underlying the effect of TPI1 on the immune microenvironment, the relationship between TPI1 expression and the expression of PDCD1 and CTLA4 was evaluated in the TIMER database. In HCC, the expression of TPI1 presented a significant positive correlation with the expression of both PDCD1 and CTLA4 (R= 0.224, P < 0.001; R= 0.204, P < 0.001; Figure 4E).

These results preliminarily demonstrated that TPI1 could participate in regulating the HCC immune microenvironment and presented a significant positive correlation with the expression of PDCD1 and CTLA4.

The Association of TPI1 Expression with CD8⁺ T Cell Infiltration in HCC

To validate the expression levels and prognostic value of TPI1 in HCC, 46 HCC specimens were analyzed by IHC and prognosis analysis (Figure 5A). The results showed that the expression of TPI1 was much higher in HCC patients with evidence of metastasis than in those without HCC metastasis (Figure 5B). In addition, patients with higher TPI1 expression presented poorer overall survival (P = 0.011, Figure 5C). The expression of CD8 was also assessed in these HCC specimens (including tumors, tumor-adjacent tissues, and normal tissues) (Figure 5D). Interestingly, the expression of CD8 was negatively correlated with the expression of TPI1 in tumor-adjacent tissues (P < 0.01, Figure 5E), but not in tumor or normal tissues samples (P > 0.05, Figure 5E). Moreover, the expression of CD8 in tumor-adjacent tissues was much higher in HCC patients without HCC metastasis than in those with evidence of metastasis (P < 0.05, Figure 5E). Immunofluorescence staining was further performed in our cohort. The results also showed that the expression of TPI1 was negatively correlated with the expression of CD8 in the tumor tissues (Figure 5F). In conclusion, these findings validated the positive correlation between TPI1 expression and HCC metastasis and unfavorable prognosis, while indicating the negative correlation with the infiltration of CD8⁺ T cells.

Figure 5 The expression of TPI1 and CD8 in HCC tissue samples. (A) Representative IHC staining of TPI1 in HCC tissues with and without metastasis; (B) The correlation between TPI1 expression and metastasis in HCC tissues; (C) The Kaplan-Meier survival curves of OS for HCC patients with different TPI1 expression; (D) Representative IHC staining of CD8 in HCC tissues and matched tumor-adjacent and normal tissue samples; (E) The correlation between the expression of CD8 and TPI1 expression and HCC metastasis in HCC tissue, matched tumor-adjacent tissue and normal tissue samples; (F) Representative IF images for TPI1 (red) and CD8 (green) staining in HCC tissues with and without metastasis. *P < 0.05; **P < 0.01. Abbreviations: HCC, hepatocellular carcinoma; IHC, immunohistochemistry; OS, overall survival; IF, immunofluorescent.

Discussion

General analysis of tumor samples may obscure critical information about changes across cells and within cell groups.¹⁶ ScRNA-seq analysis is a novel technique to better understand the collective behavior and mutual regulated communication patterns of cells in the liver ecosystem.¹³ For example, a pioneering study has revealed that a cluster of LAMP3⁺ DCs could regulate lymphocytes in the TME through crosstalk and migrate to the lymph nodes, which is associated with dysfunctional T cells.¹⁷ Another novel study conducted deep scRNA-seq on T cells from HCC patients, suggesting that specific subpopulations, such as exhausted CD8⁺ T cells and Tregs, might be preferentially enriched and clonally expanded in HCC.¹⁸ Currently, little is known about the relationship between malignant hepatocyte clusters and HCC metastasis at the single-cell level. In this study, we collected single-cell transcript level sequencing data of 14,779 single cells from HCC patients and analyzed metastasis-related genes in 9 malignant hepatocyte clusters. C3_Hepatocyte-SERF2 and C6_Hepatocyte-IL13RA2 were found to significantly up-regulated in the metastasis group, and were mainly involved in the regulation of immune microenvironment. Furthermore, we investigated the key gene in malignant hepatocyte clusters and its relationship with immune cell infiltration. The results revealed that TPI1, a key gene in C6_Hepatocyte-IL13RA2 and HCC metastasis, could participate in regulating the HCC immune microenvironment. Specifically, it was negatively correlated with the infiltration of CD8⁺ T cells, which may contribute to HCC metastasis. This provides a new perspective to explore potential biomarkers for immunotherapy and the underlying mechanisms between HCC metastasis and tumor immune microenvironment (TIME).

HCC tumors are highly heterogeneous, so patients may respond differently to the same treatment, resulting in poor treatment outcomes.^10,19 Therefore, we explored the intratumor heterogeneity of HCC to further understand the HCC microenvironment, and thus improving the prognosis of the patients. HCC tissue mainly contains nine major cell types, such as immune cells, endothelial cells, and hepatocytes.²⁰ It has been found that HCC originates from hepatocytes, and changes in the characteristics of hepatocytes not only affect their own proliferation, but also affect the TME through intercellular communication and other regulatory ways.¹⁹ According to our results, C3_Hepatocyte-SERF2 and C6_Hepatocyte-IL13RA2 were mainly involved in the regulation of immune microenvironment, which was also a significant pathway in regulating HCC metastasis.

The tumor microenvironment of HCC is a complex and dynamic ecosystem that generates an immunosuppressive milieu, resulting in HCC initiation, invasion, and metastasis.²¹ It could coexist and interact with immune cells and their products, which is called TIME.²² To further explore the possible mechanism of HCC metastasis, we investigated the key gene in the malignant hepatocyte cluster and its role in the TIME. In the present study, the results showed that TPI1, a key gene in C6_Hepatocyte-IL13RA2, was associated with the poor prognosis and HCC metastasis. Accumulating research has uncovered that TPI1 could be a key enzyme in glycolysis and a candidate oncogene.²³ It is overexpressed in a variety of cancers, including intrahepatic cholangiocarcinoma, stomach cancer, and breast cancer.^24–26 However, the specific role of TPI1 in HCC remains to be identified. According to the results, the expression of TPI1 was negatively correlated with the infiltration of CD8⁺ T cells, DC cells and B cells, while it was positively correlated with the infiltration of Tregs and MDSCs. Recently, researches have revealed that various immune cells in HCC microenvironment, including cytotoxic CD8⁺ T cells and DC cells, play a crucial role in maintaining immune surveillance.^22,27 Moreover, several immune cells such as MDSCs and Tregs could contribute to immune evasion and tumor progression in HCC.^28,29 Therefore, we believed that TPI1 might contribute to poor prognosis and HCC metastasis via inhibiting the tumor immune microenvironment. Furthermore, cancer immunotherapy in recent years has greatly changed the treatment paradigm of HCC, especially the application of immune checkpoint inhibitors (ICIs).²² The efficacy of ICIs depends not only on the TIME, but also on the adequate expression of immune checkpoints.³⁰ PD-1 and CTLA-4 are the two most common immune checkpoints in HCC, mainly providing immunosuppressive signals. In the present study, TPI1 expression was found to be positively correlated with the expression of PDCD1 and CTLA4, suggesting that targeting TPI1 could facilitate the efficacy of immune checkpoint inhibitors in HCC.

The effect of TPI1 on the immune microenvironment of HCC might be related to glycolysis. Upregulation of TPI1 has been proved to cause abnormally upregulated glycolysis in various cancers, including HCC.^26,31,32 Increased glycolysis is often hijacked by cancer cells to promote high energy production, thus becoming a key mechanism contributing to the high rates of tumor growth.^33,34 In addition, glycolysis could regulate the TIME by inhibiting immune cell function and promoting tumor immune escape.^35,36 Specifically, lactate is the end product of glycolytic metabolism, which is essential for tumor cells to maintain stemness and evade immune surveillance.³⁷ Accumulated lactate in the TME could lead to defective antigen presentation, restricted immune cell function, and insufficient immune cell infiltration.^38–40 In addition, a pioneering study has revealed that tumor-derived lactate could rewire CD8⁺ T cell metabolism.⁴¹ Low pH and high lactate concentration has also been found to impair the metabolic activity, proliferation and cytokine production of CD8⁺ T cells.⁴² Consistent with these studies and the bioinformatics results, our study verified that TPI1 was negatively correlated with the infiltration of CD8⁺ T cells, which may contribute to HCC metastasis.

There are still some limitations in the present study. The potential mechanisms of TPI1 on HCC metastasis remain to be elucidated with further functional studies in vitro and in vivo. Moreover, our results were mainly based on the TCGA, GEO, and TIMER databases, and although they have been validated in clinical samples, the sample size should be expanded.

Conclusion

In summary, our study revealed that a cluster of TPI1⁺ malignant hepatocytes was associated with the suppression of CD8⁺ T cell infiltration and HCC metastasis, providing novel insights into potential biomarkers for immunotherapy in HCC.