Development of a Starvation Response-Based Model and Its Application in Prognostic Assessment of Liver Hepatocellular Carcinoma

Abstract

Background: Hepatocellular carcinoma (LIHC) is a highly prevalent and poorly prognostic malignancy worldwide, and nutrient deprivation in the tumor microenvironment activates the starvation response in tumor cells. Starvation response-related genes (SRRGs) play critical roles in maintaining energy metabolism and promoting tumor development, but their value in prognostic prediction of LIHC has not been clarified. Methods: We based on public databases to obtain transcriptome and single-cell RNA sequencing (scRNA-seq) data for LIHC and SRRG from previous studies. Key modules relevant to SRRGs were identified by weighted gene co-expression network analysis (WGCNA). Functional enrichment analysis was conducted using clusterProfiler R package. Independent prognostic genes were screened to build a RiskScore model and its performance was further verified. The immune microenvironmental profile of patients in different risk groups was assessed using the single-sample gene set enrichment analysis (ssGSEA), MCP-Counter, ESTIMATE, and TIMER algorithms. Seurat package for single-cell profiling and validation of key gene expression based on Huh7 and transformed human liver epithelial-2 (THLE-2) cell lines. The LIHC cell migration and invasion were measured by conducting wound healing and transwell assays. Results: The key module identified by WGCNA showed the strongest correlation with SRRGs and the glycolysis-related SRRGs were mainly enriched in metabolism-correlated pathways. Two protective genes (FBXL5 and PON1) and three risk genes (TFF2, TBC1D30, and SLC2A1) were discovered as the independent prognostic genes for LIHC. Activation of cytokine-cytokine receptor interaction and IL-17 signaling pathway and higher infiltration of immune cells in high-risk group was observed. The five independent prognostic genes were mainly expressed in cancer stem cells and epithelial cells, in particular, SLC2A1 and TFF2 were significantly high-expressed in epithelial cells in the tumor group than in nontumor group. FBXL5 and PON1 were downregulated, while TFF2, TBC1D30, and SLC2A1 were upregulated in LIHC cells. Silencing SLC2A1 significantly inhibited LIHC cell migration and invasion. Conclusion: In this study, we constructed the first risk model based on SRRGs to accurately predict the prognosis of LIHC, which provides a new idea for individualized treatment and targeted intervention.

Keywords: RiskScore model; glycolysis; liver hepatocellular carcinoma; prognosis; single-cell RNA sequencing; starvation response; transcriptome.

Figure 1

Identification of key gene module associated with starvation response-related genes (SRRG) in LIHC by WGCNA. (A) Scale-free fitting index analysis of soft threshold (power). (B) Cluster dendrogram of hierarchical clustering. (C) Module-trait relationships between each module and SRRG score. (D) Number of genes of each module.

Figure 2

Functional enrichment analysis of the glycolysis-related SRRG. (A) Venn diagram of glycolysis-related genes (GRG) and blue module genes. (B) KEGG enrichment pathways of glycolysis-related SRRG. (C) GO enrichment terms in biological process (BP). (D) GO enrichment terms in cellular component (CC). (E) GO enrichment terms in molecular function (MF).

Figure 3

A RiskScore model developed based on independent prognostic genes and validation of the model. (A) LASSO Cox regression analysis to reduce gene range. (B) Randomforest plot of multivariate Cox regression analysis. (C) ROC curve of RiskScore model in TCGA-LIHC cohort. (D) Kaplan–Meier (K–M) curve of overall survival (OS) between different risk groups in TCGA-LIHC cohort. (E) Percentage of survival status between different risk groups in TCGA-LIHC cohort. (F) ROC curve of RiskScore model in ICGC-LIRI-JP dataset. (G) K–M curve of survival probability between different risk groups in ICGC-LIRI-JP dataset. (H) Percentage of survival status between different risk groups in ICGC-LIRI-JP dataset. (I, J) Expression levels of independent genes for LIHC prognosis in TCGA-LIHC and ICGC-LIRI-JP datasets.

Figure 4

Immune cell infiltration between different risk groups. (A–C) Infiltration levels of immune cells between low- and high-risk groups calculated by ESTIMATE, TIMER, and ssGSEA algorithms. (D) Correlation between independent prognostic genes, RiskScore and infiltration levels of immune cells calculated by MCP-counter. (E) Expressions of immune checkpoint genes in the risk groups. ⁣^∗∗∗∗ means p < 0.0001; ⁣^∗∗∗ means p < 0.001; ⁣^∗∗ means p < 0.01; ⁣^∗ means p < 0.05; ns means not significant.

Figure 5

GSEA of different risk groups in TCGA-LIHC cohort. (A) KEGG enrichment pathways in high-risk group. (B) KEGG enrichment pathways in low-risk group.

Figure 6

Construction of nomogram and validation. (A) Distribution of RiskScore in different clinical characteristics such as gender, AJCC stage, grade, and age. (B, C) Univariate and multivariate Cox regression analysis of RiskScore and clinical features. (D) Nomogram model constructed by combining RiskScore and AJCC stage; ⁣^∗∗∗ means p < 0.001. (E) Calibration curved of the nomogram for 1, 3, and 5 year(s). (F) Decision curve analysis (DCA) of nomogram.

Figure 7

Single-cell atlas of LIHC. (A) UMAP plot of 18 primary tumor and adjacent nontumor liver samples. (B) UMAP dimensionality reduction plot of single-cell clustering and annotation. (C) Bubble plot of marker genes in each cell cluster. (D) Percentage of each cell cluster in tumor and nontumor groups. (E) Expression levels of independent prognostic genes in each cell cluster of tumor group. (F) Expressions of SLC2A1 and TFF2 in epithelial cells.

Figure 8

The in vitro validation by LIHC cells. (A) qRT-PCR detecting the relative mRNA expressions of five prognostic SRRGs. (B) CCK-8–based assay to assess the effect of silencing SLC2A1 on LIHC cell viability. (C) Transwell assay assessing the impact of SLC2A1 silencing on cell invasion of Huh7 cells. (D) Wound healing assay evaluating the effect of SLC2A1 silencing on cell migration of Huh7 cells. All procedures were repeated three times independently. ⁣^∗∗∗∗ means p < 0.0001; ⁣^∗∗∗ means p < 0.001; ⁣^∗∗ means p < 0.01; and ⁣^∗ means p < 0.05.