Transcriptomic data and clinical information of HCC patients were obtained from TCGA-LIHC (The Cancer Genome Atlas database-Liver Hepatocellular Carcinoma, n = 370; https://www.cancer.gov/ccg/research/genome-sequencing/tcga), the International Cancer Genome Consortium (ICGC-LIHC, n = 232), which were downloaded from UCSC XENA (https://xenabrowser.net/). GSE148355 (n = 52) was derived from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Hypoxia and glycolysis-related gene sets were obtained from the Molecular Signature Database (GSEA $|$ MSigDB (https://gsea-msigdb.org)), and lactylation-related gene sets were obtained from a previously published study [17]. To investigate the spatial expression of each gene in the signature, spatial transcriptomic data was downloaded from a previous study [22] (https://ngdc.cncb.ac.cn/gsa-human/browse/HRA000437), and subsequent analysis was performed using the ‘Seurat’ R package (https://cran.r-project.org/web/packages/Seurat/index.html, v5.1.0).

The hypoxia, glycolysis and lactylation genesets were derived from a previously published study [17] and the Molecular Signature Database (GSEA $|$ MSigDB (https://gsea-msigdb.org)). We obtained the differentially expressed genes (DEGs) in HCC patients versus normal individuals with the ‘DEseq2’ R package (https://github.com/thelovelab/DESeq2, v1.38.3), by the criteria of FDR 0.05 and $|$log2fc$|$ $\ge$1.

Employing univariate COX analysis, prognostic HGL genes were obtained, then Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/), Gene Ontology (GO, https://www.geneontology.org/), Gene Set Enrichment Analysis (GSEA, https://www.gsea-msigdb.org/gsea/index.jsp) enrichment analysis were performed.

The TCGA-LIHC patients were divided into two clusters with prognostic HGL genes using the ‘CancerSubtypes’ R package (https://github.com/taoshengxu/CancerSubtypes, v.1.24.0).

To compare the immune microenvironment between the two clusters, we predicted differences in infiltration of 6 immune cells using the Tumor Immune Estimation Resource (TIMER) algorithm (http://cistrome.org/TIMER/). The difference of immune function was predicted by the ‘GSVA’ R package (https://github.com/rcastelo/GSVA, v1.46.0) [23]. The immune evasion ability of the two clusters was evaluated by the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/).

To assess the guiding potential of HGL clustering for clinical drug use, we measured the ‘PRROPHETIC’ R package (https://github.com/paulgeeleher/pRRophetic, v0.5) to predict IC50 values for each drug and compared differences between the two clusters for commonly used drugs in HCC.

We obtained DEGs between 2 clusters with the ‘Limma’ R package (https://www.bioconductor.org/packages/release/bioc/html/limma.html, v3.54.2), the criteria were $|$log2fc$|$ $\ge$1, p 0.05; then measured COX LASSO regression screening for genes to construct the optimal model. The ‘Survminer’ R package (https://cran.r-project.org/web/packages/survminer/index.html, v0.5.0) was used for Kaplan-Meier (K-M) survival analysis to compare the difference of Overall Survival (OS) between high-risk and low-risk groups in different data sets. The predictive ability of this signature was further evaluated using ‘timeROC’ R package (https://cran.r-project.org/web/packages/timeROC/, v0.4.0) to calculate the receiver operating characteristic (ROC) curve and the area under the curve (AUC).

We explored the association between the HGL signature (HGLs) and clinical features, finding that HGLs had an independent ability to predict patients’ outcomes by K-M analysis.

A total of 8 pairs of HCC adjacent and tumor tissues were obtained from the Xiangya Hospital, Central South University. The Central South University Xiangya Hospital Ethics Committee approved this study. All patients or their families/legal guardians obtained their informed consent to use their materials. The patient’s clinical information is presented in Supplementary Table 1.

The tissue samples from patients were cryogenically ground, followed by RNA extraction using the TRIZOL reagent (ET101-01-V2, Transcript®, Beijing, China). The extracted RNA was then reverse transcribed into cDNA using the reverse transcription kit (E047-01B, NovoScript®, Suzhou, China). Finally, the PCR reaction system mixed with primers and cDNA was prepared. Primers used are listed in Supplementary Table 2. To calculate the relative expression of HOMER1, we used $\beta$-actin as a reference gene, the 2^{-Δ⁢Δ⁢Ct} method is used: $\mathrm{\Delta }$Ct = Ct_HOMER1 – Ct_β-actin, $\mathrm{\Delta }\mathrm{\Delta }$Ct = $\mathrm{\Delta }$Ct_Tumor – $\mathrm{\Delta }$Ct_Normal, Relative Expression = 2^{-Δ⁢Δ⁢Ct}.

Tissue or cell lysates were prepared using lysis buffer and protein concentrations were measured. The protein was separated by SDS-PAGE and transferred to a nitrocellulose membrane. The membrane was blocked with 5% milk, followed by incubation with the primary antibody overnight at 4 °C. After washing with Tris-Borate-Sodium Tween-20 (TBST) buffer, the membrane was incubated with the secondary antibody for 2 hours. Protein was visualized using chemiluminescent detection. Antibodies used were as follows: HOMER1 (1:1000, db12743, Diagbio, Hangzhou, China), $\beta$-Tubulin (1:10,000, 66240-1-Ig, Proteintech, Wuhan, China), Peroxidase-conjugated AffiniPure Goat Anti-Mouse IgG (H+L) (1:5000, 115-035-003, Jackson, West Grove, PA, USA), Peroxidase-conjugated AffiniPure Goat Anti-Rabbit IgG (H+L) (1:5000, 111-035-003, Jackson, West Grove, PA, USA).

LM3 and Hep3B in this study were purchased from BDBIO (Zhejiang, China). The two types of cells have been successfully authenticated by STR profiling and have passed mycoplasma testing. They were cultured in DMEM medium (10% FBS) under 5% CO₂ at 37 °C.

According to the knocking down of HOMER1, Hep3B and LM3 cell lines were transfected with lipofectamine 2000 following the manufacturer’s instructions. All small interfering RNA (siRNA) fragments were designed and synthesized by Gene Pharma (Shanghai, China); and are listed in Supplementary Table 3.

After knocking down HOMER1 in Hep3B and LM3 cells as mentioned above, they were incubated for 72 hours. Then the cells were digested and collected in a cell suspension, and washed with PBS. The cells were fixed by incubating with 70% ethanol for 12 hours. After removing the ethanol, the cells were permeabilized using 0.1% Nonidet P-40 (NP-40) (N8030, Solarbio®, Beijing, China). Finally, the cells were stained with propidium iodide for 30 minutes (The propidium iodide solution is prepared by combining 3 mL of staining buffer, 150 µL of 20$\times$ propidium iodide staining solution, and 60 µL of 50$\times$ RNase A), and analyzed at different cell cycle stages using flow cytometry based on the fluorescence intensity.

Hep3B and LM3 cells were seeded in 6 cm dishes, and transfected with a mixture of 2 µL of lipofectamine 2000, 200 µL of Opti-MEM, and 3 µL of siHOMER1 to knock down the expression of HOMER1 in both cell lines; subsequently 1 $\times$ 10⁵ cells were seeded into each six-well plate; scratches were made with 10 µL pipette until cells were 100% the following day. The images were taken at 0 h and 24 h. The statistical results were analyzed by ImageJ (V1.54g, National Institutes of Health, Bethesda, MD, USA).

After transfected with HOMER1 siRNA or random control, 500 cells were seeded into six-well plates, and 10 days later fixed with 4% paraformaldehyde (G1101, Servicebio, Wuhan, China) and stained with crystal violet solution. The statistical results were analyzed by ImageJ.

We used a 24-hole Transwell chamber (Transwell, Corning Costar, Corning, NY, USA). For migration assays, 2 $\times$ 10⁴ cells were seeded in the upper chamber, cultured in serum-free DMEM, and 20% serum DMEM was added in the lower chamber. For the invasion assay, 30 µL Matrigel was used to coat the superior compartment, and then 2 $\times$ 10⁵ cells were seeded into the superior compartment. After 48 hours of culture, the cells were fixed with 4% paraformaldehyde, stained with crystal violet, then swabbed to remove the upper cells. The statistical results were analyzed by ImageJ.

We obtained 149 HGL-DEGs between tumor and normal tissues (Fig. 1A,B). Univariate COX regression showed that 72 genes were significantly associated with the prognosis of HCC. The GO and KEGG enrichment analyses revealed that prognostic-associated HGL genes were associated with metabolism-related pathways, including glycolysis, the HIF-1 signaling pathway, and organic acid catabolic and carboxylic acid catabolic processes (Fig. 1C,D). It reveals that both hypoxia, glycolysis and lactylation effect the prognosis of HCC.

Fig. 1.

HGL associate genes differentially expressed in TCGA-LIHC. (A) The volcano plot showing the HGL-DEGs between tumor and normal samples in TCGA dataset. (B) Heatmap showing top 30 DEGs. (C) The bar plot of KEGG enrichment pathway of prognostic-associated HGL genes. (D) GO enrichment analysis bubble plot. HGL, hypoxia-glycolysis-lactylation; TCGA-LIHC, The Cancer Genome Atlas database-Liver Hepatocellular Carcinoma; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; NS, no significance.

Based on the prognostic-associated HGL genes, we conducted a consensus clustering algorithm on the TCGA-LIHC dataset to further investigate the effect of HGL on HCC. When K = 2, we obtained the appropriate clustering result (Fig. 2A). Principal component analysis (PCA) also displayed that the two clusters were separated visibly (Fig. 2B). K-M curves revealed that patients in cluster1 had better median survival time than those in cluster2 (Fig. 2C, p 0.001). GSEA analysis revealed that pathways enriched in cluster1 were mainly metabolic-related pathways, while cluster2’s were related to the malignant process of HCC. These pathways include: E2F_TARGETS, G2M_CHECKPOINT, EPITHELIA_MESENCHYMAL_TRANSITION (Supplementary Fig. 1A,B).

Fig. 2.

HGL clusters based on consensus clustering. (A) Consensus clustering divides 370 patients in TCGA dataset into two distinct clusters. (B) PCA analysis shows that the two clusters of patients are significantly segmented. (C) Kaplan-Meier (K-M) curve analysis of cluster1 and cluster2. (D) Infiltration differences of six immune cells between cluster1 and cluster2. (E) The levels of immune function of the two HGL clusters. (F) The immune checkpoint receptor genes are highly expressed in cluster2. (G) Cluster2 has higher TIDE score. Statistical significance was denoted as follows. *p 0.05, **p 0.01, ***p 0.001. PCA, principal component analysis; TIDE, Tumor Immune Dysfunction and Exclusion.

There were variances in immune cells between the two clusters (Fig. 2D). The check-point function of cluster2 was significantly higher than that of cluster1 (Fig. 2E, p 0.05). PDCD1, CLTLA4, LAG3, and Havcr2, the major inhibitory immune checkpoint receptors in HCC [24], were highly expressed in cluster2 (Fig. 2F, p 0.001). A TIDE score of cluster2 was also higher than that of cluster1 (Fig. 2G, p 0.01), which suggests that cluster1 may have a better response to immunotherapy.

Chemotherapy and immunotherapy are both important therapeutic methods for HCC. After IC50 prediction of chemotherapeutic drugs, we found that patients in cluster1 were more sensitive to sorafenib than patients in cluster2 (p 0.001). The same is true for other drugs, such as: 5-Fluorouracil, etoposide, cytarabine, docetaxel, cisplatin, doxorubicin, nilotinib, paclitaxel, gemcitabine, gefitinib, sunitinib, vinblastine (p 0.001). Overall, the drug sensitivity of cluster1 was superior to cluster2 (Fig. 3).

Fig. 3.

The drug sensitivity of cluster1 is better than cluster2. **p 0.01, ***p 0.001.

We performed differential analysis between the two clusters, and after conducting univariate COX regression analysis on the differentially expressed genes, we finally identified these 123 survival-related differential genes (Fig. 4A). Using LASSO and multivariate COX regression, we acquired 13 genes to construct the signature (Supplementary Fig. 2A). The riskScore calculation was as follows: riskScore = (0.20233 $\times$ CDCA8) + (0.140306 $\times$ GAGE2D) + (0.0913041 $\times$ S100A9) + (–0.08057605 $\times$ SLC16A11) + (–0.1422534 $\times$ SOCS2) + (0.0711202 $\times$ FAM163A) + (–0.08746885 $\times$ LGSN) + (0.09041832 $\times$ FRMD1) + (0.09532562 $\times$ HOMER1) + (–0.1636347 $\times$ CD8B) + (–0.1596561 $\times$ CD27) + (0.1255294 $\times$ E2F5) + (01694574 $\times$ PGF).

Fig. 4.

Construction of HGL signature in LIHC. (A) Parameter selection for LASSO. (B) Kaplan-Meier survival curve analysis between high risk and low risk groups in TCGA-LIHC datasets. (C) Risk curves: overall survival status, and the heatmap of the expression of 13 HGL signature associate genes. The dotted line represents the division of the median. (D) ROC curves of the 1–5 years survival rates of TCGA-LIHC dataset. (E) Kaplan-Meier survival curve analysis between high risk and low risk groups in ICGC-LIHC datasets. (F) Kaplan-Meier survival curve analysis between high risk and low risk groups in GSE148355 datasets. LASSO, Least Absolute Shrinkage and Selection Operator; ROC, receiver operating characteristic.

The patients in the TCGA training set were divided into a high-risk group (N = 185) and low-risk group (N = 185) based on the riskScore, while the OS of the low-risk group was significantly better than the high-risk group (Fig. 4B, p 0.001). The risk curve also showed that as the risk score increases, the survival rate decreases (Fig. 4C). At 1–5 years, the AUC of the ROC curve for risk signature was 0.818, 0.824, 0822, 0.829, and 0.821, respectively, which indicates a great predictive accuracy for HCC patients (Fig. 4D).

To further validate the predictive effectiveness of the signature, we grouped the ICGC and the GSE148355 together, since the validation set was similar: ICGC high-risk group (N = 105), ICGC low-risk group (N = 127), and GSE148355 high-risk group (N = 27), low-risk group (N = 25). Patients in the low-risk group still had significantly better OS than those in the high-risk group (Fig. 4E,F, p 0.001/p = 0.01). At 1–4 years, the AUC of the ROC curve for ICGC-LIHC were 0.700, 0.720, 0.738, 0.750 (Supplementary Fig. 2B), respectively; the AUC for GSE148355 were 0.758, 0.689, 0.707, 0.710 (Supplementary Fig. 2C), respectively. The signature also exhibited good prediction ability for these external validation cohorts.

HGL signature exhibits predictive power, thus we desired to explore whether it can be used as an independent prognostic factor for HCC patients. Age, gender, alpha-fetoprotein level (AFP), microvascular invasion (MVI), and pathological stage were crucial factors affecting the prognosis of HCC patients, riskScore could also accurately predict the prognosis. In the K-M curve, patients in different ages, different sex, and in different pathological stages, with different expression of AFP, whether with microvascular invasion or not, showed a significant variance in OS between the high and low risk group (Fig. 5A, p 0.01). Univariate and multivariate COX analysis showed that the riskScore and pathologic stages were related with the OS of patients according to the univariate COX analysis (Fig. 5B); while the multivariate COX analysis revealed they are also the independent predictors (Fig. 5C).

Fig. 5.

HGL riskScore is an independent prognostic factor. (A) The Kaplan-Meier survival analysis was performed based on HGL risk signature for patients with different gender, with different ages, different AFP values, different MVI and pathological stages. (B) Forest plot showing Univariate-COX analysis of the TCGA-LIHC cohort including age, AFP, MVI, pathologic stage, and riskScore as variables. (C) Forest plot showing Multivariate-COX analysis including Pathologic stage, and riskScore as variables in the TCGA-LIHC cohort. AFP, alpha-fetoprotein level; MVI, microvascular invasion.

To obtain a better prediction of prognosis, we incorporated the above clinical factors together with the riskScore to construct a nomogram. Patients with a higher score represent worse survival (Fig. 6A, p 0.001). The C-index of each factor suggests that the riskScore still contributes the most to the nomogram (Fig. 6B). The survival ROC analysis shows that the AUC were 0.822, 0.676 for riskScore and pathologic stage, which were superior to other clinical factors. The prognostic ability of gender and age is relatively poor, and there are also fewer patients with MVI data in the TCGA-LIHC cohort. As a result, the AUC of these clinical prognostic factors is not ideal (Fig. 6C). The calibration curve revealed that the predicted OS and observed OS have great consistency, and that nomogram has good predictive ability (Fig. 6D).

Fig. 6.

Construction and validation of the nomogram. (A) Nomogram for predicting 1-, 3-, and 5-year OS. (B) Time-dependent C-index of riskScore and other clinical factors. (C) ROC curves of riskScore and other clinical factors predicting 3-year OS of patients. (D) Calibration curve for nomogram. *p 0.05, ***p 0.001. OS, Overall Survival.

To explore the expression patterns in HCC, we incorporated a set of spatial transcriptomic data to investigate genes in HGL signature from both the expression level and spatial orientation [22]. The organization of four sets of spatial transcriptomes were taken from the junction of tumor and normal tissue, and cells obtained from these samples were divided into 14 cell clusters (Fig. 7A). Based on the annotations that come with the spatial transcriptome data, these cells could be defined as five kinds of localized cell populations, namely: Hepatocytes, Immune infiltration region (IIR), Transition, Tumor, Others (Fig. 7B).

Fig. 7.

Spatial transcriptome data revealing HOMER1 significantly enriched in hepatocellular carcinoma cell. (A) Four spatial transcriptome samples taken from the junction of hepatocellular carcinoma and para-carcinoma were analyzed by dimensionality reduction and cluster analysis, 14 different clusters were identified. (B) We identified clusters derived from Hepatocytes, Tumor, Immune infiltrate region (IIR), Transition, Others these five regions. (C) Bubble diagram showing genes expression in each region. (D) Spatial distribution of different cell clusters. (E) Spatial distribution of HOMER1 is concentrated in tumor cells. UMAP, uniform manifold approximation and projection.

As shown in the bubble diagram, among the genes in the signature, HOMER1 is mainly expressed in tumor-derived cells (Fig. 7C). To further validate our findings, we performed spatial transcriptomic analysis on these four sets of samples. From the spatial distribution of the cell clusters and the spatial distribution of HOMER1, we found that HOMER1 highly coincides with the distribution of tumor cells (Fig. 7D,E). The differential expression of HOMER1 indicates that this protein may play a key role in the malignant process of hepatocellular carcinoma.

The mechanism of HOMER1 in HCC has not been explored, thus we sought to explore the actual expression of HOMER1 in HCC. Eight pairs of adjacent and tumor tissues from HCC patients were obtained. Western blot revealed HOMER1 was elevated in tumor tissues (Fig. 8A,B, p 0.01); RT-qPCR also confirmed its presence at the transcription level (Supplementary Fig. 2D).

Fig. 8.

HOMER1 is highly expressed in tumor tissues and affects the malignant proliferation of hepatocellular carcinoma cell lines. (A) Western Blot revealed that the expression of HOMER1 was evaluated in tumors and in adjacent tumors at the protein level. (B) The comparison of HOMER1 protein expression levels between normal and tumor tissues. (C) HOMER1 was successfully knocked down in LM3 and Hep3B cell lines. (D) The proliferation ability of LM3 and Hep3B decreased after knock down of HOMER1. (E) Knock down of HOMER1 blocked the cell cycle of LM3 and Hep3B. **p 0.01, ***p 0.001, ****p 0.0001.

To explore its function in hepatocellular carcinoma cell lines, we knocked down HOMER1 in LM3 as well as in Hep3B cell lines, with Western-Blot demonstrating that HOMER1 knock down was very successful (Fig. 8C; Supplementary Fig. 2E,F). We observed that knocking down HOMER1 significantly inhibited the proliferation of LM3 and Hep3B. The colony formation assay clearly exhibited the difference of cell proliferation between untreated and si-HOMER1 cells (Fig. 8D, p 0.0001). Cell cycle analysis further revealed that knocking down HOMER1 increased the proportion of cells arrested in the G2/M phase (Fig. 8E). These results suggested that HOMER1 may play an important role in the cell cycle of hepatocellular carcinoma cell lines. The wound healing assay showed that HOMER1 knock down also inhibited the migration of LM3 and Hep3B cells (Fig. 9A, p 0.01), while the Transwell assay confirmed both migration and invasion of LM3 and Hep3B cells were decreased after knocking down of HOMER1 (Fig. 9B,C, p 0.01). These experiments revealed knocking down HOMER1 could inhibit the proliferation of hepatocellular carcinoma cell lines by blocking the cell cycle, and also affect the invasion and metastasis of hepatocellular carcinoma cell lines.

Fig. 9.

The knocking down of HOMER1 affected the ability of metastasis and invasion of hepatocellular carcinoma cell lines. (A) The wound healing assay showed that knocking down of HOMER1 significantly decreased the migration ability of Hep3B and LM3. Scale bars = 100 µm. (B) The Transwell migration assay showed that knocking down of HOMER1 decreased the migration ability of Hep3B and LM3. Scale bars = 50 µm. (C) The Transwell invasion assay showed that knocking down of HOMER1 decreased the invasion ability of Hep3B and LM3. Scale bars = 50 µm. **p 0.01, ***p 0.001, ****p 0.0001.