Development of a prognostic histamine-based risk score signature for HCC

The differential analysis of HCC and paracancerous tissue samples identified 1873 genes (Fig. 1A). Integration of histamine-related signaling pathways from the database revealed 1541 specific histamine-related genes. After intersecting these datasets, we identified 143 histamine-related DEGs for further analysis (Fig. 1B). After employing GO and KEGG analysis for functional pathway enrichment, the primary enriched pathways were the neuroactive ligand–receptor interaction, PI3K–Akt signaling pathway, and cytokine–cytokine receptor interaction (Fig. 1C). Histamine-related genes were predominantly enriched in the biological processes such as cellular divalent inorganic cation homeostasis, calcium ion homeostasis, and cellular calcium ion homeostasis signaling pathway (Fig. 1D). The Cox univariate regression and LASSO analyses (p < 0.05) revealed six histamine-related genes that have a significant impact on the prognosis of patients with HCC (Fig. 1E, F). EZH2 and FLVCR1 were determined to be independent risk factors affecting the prognosis of HCC through Cox multivariate regression analysis (p < 0.05; Fig. 1G). Risk scores were calculated for all the patients in TCGA training cohorts using a defined formula: RiskScore = 3.599e – 04 × FLVCR1 expression + 5.473e – 04 × EZH2 expression. The median risk score was used as the cutoff value to classify patients into high- and low-risk groups (Fig. 1H). Survival analyses revealed that patients in the high-risk group had worse overall survival (OS) than those in the low-risk group in the TCGA training cohort (p < 0.001; Fig. 1I). The area under the curve (AUC) for OS at respective time points was as follows: one year, 0.77; three years, 0.71; and five years, 0.69 (Fig. 1J). The heatmap demonstrated the overexpression of most histamine-related genes in the high-risk group (Fig. 1K), suggesting that high expression levels of these genes characterize this group.

Fig. 1

Construction of a prognostic signature based on the differentially expressed histamine-related genes. A Volcano map showing differentially expressed genes in HCC and paracancerous tissue samples. B Identification of differentially expressed histamine-related genes. C Results of KEGG functional enrichment analysis of 1541 histamine-related genes. D Results of GO functional enrichment analysis. E The coefficients of genes were calculated by multivariate Cox regression using LASSO. F The partial likelihood of deviance of genes. G Results of Cox multivariate analysis of six differentially expressed histamine-related genes. H The association of risk scores with survival status and gene expression in patients with HCC in the training set. I Kaplan–Meier curves were used to compare the overall survival of patients with HCC between the high- and low-risk groups in the training set. J ROC curves of the prognostic signature for predicting the risk of death at 1, 3, and 5 years in the training set. K The heatmap shows 1541 histamine-related differentially expressed genes in the high-risk and low-risk groups

Relationship between HRGPS and hypoxic microenvironment

The hypoxic microenvironment is critical in immune cell suppression and tumor cell invasion in the TME. GSEA enrichment analysis revealed that the high-risk group was predominantly enriched in the hypoxia signaling pathway (Fig. 2A). The majority of genes associated with hypoxia were found to be overexpressed in the high-risk group (Fig. 2B), including VEGFA and ANKZF1 in the hypoxia pathway (Fig. 2C, D) [28, 29]. Correlation analysis indicated a strong relationship between mRNA expression levels of VEGFA and ANKZF1 and histamine-related genes (Fig. 2E, F). These findings suggest that overexpression of histamine-related genes may result from the hypoxic microenvironment in tumor tissues.

Fig. 2

Relationship between the prognostic signature and hypoxic microenvironment. A The GSEA functional enrichment analysis results show that the hypoxia signaling pathway enriched the high-risk group. B Heatmap demonstrating that most hypoxia-related genes are overexpressed in high-risk groups. C–D The mRNA expression levels of two key genes (VEGFA and ANKZF1) in the hypoxic microenvironment were significantly higher in the high-risk group. E The mRNA levels of the EZH2 were positively associated with those of the VEGFA. F The mRNA levels of the EZH2 were positively associated with those of the ANKZF1. G The mRNAsi levels were significantly higher in the high-risk group than in the low-risk group. H–I A significant positive correlation was observed between the mRNA expression levels of histamine-related genes (EZH2 and FLVCR1) and mRNAsi scores

Previous studies have suggested that the hypoxic tumor microenvironment could affect cancer stem-like cell phenotype [30, 31], prompting us to investigate the association between histamine-related genes and mRNAsi score. Notably, the high-risk group exhibited higher mRNAsi scores than the low-risk group (Fig. 2G). Correlation analysis also demonstrated a positive correlation between the mRNA expression levels of both EZH2 and FLVCR1 with mRNAsi scores (Fig. 2H, I). These findings suggest that in the hypoxic environment of HCC, overexpression of genes related to histamine may affect the phenotype of cancer stem cell-like cells.

Association between HRGPS and tumor immunity

We investigated the correlation between histamine and the immune status of patients in the TCGA cohort and observed significant alterations in the immune cell populations. Interestingly, we noted a substantial infiltration of immunosuppressive cells in the TME, including Tregs, macrophages, dendritic cells, and neutrophils, particularly within the high-risk group, indicating pronounced immunosuppression (Fig. 3A). Similar outcomes were derived from ssGSEA enrichment analysis (Fig. 3B). The findings of the study of immune cell infiltration indicate that histamine-related genes could play a role in developing an immunosuppressive microenvironment by recruiting immunosuppressive cells and helping tumor cells avoid immune surveillance. To further validate this result, we applied the same methodology to calculate immune cell infiltration abundance in a validation set comprising 220 patients with HCC (Fig. 3C, D). The results were consistent with those obtained in the training set. Correlation analysis revealed a positive association between EZH2 and FLVCR1 expression levels with infiltrating Treg cells and macrophages (Fig. 3E–G). Utilizing the ESTIMATE algorithm, it was determined that individuals within the high-risk group exhibited lower StromalScore and ImmuneScore but higher TumorPurity (Fig. 3H–J).

Fig. 3

Relationship between the prognostic signature and immune microenvironment. A Differences in the proportion of infiltrating immune cells between the high- and low-risk groups in the TCGA database. B Differences in the proportion of 28 immune cells infiltrates in the TME between high- and low-risk groups in the TCGA database. C Differences in the proportion of infiltrating immune cells between the high- and low-risk groups in the ICGC database. D Differences in the proportion of 28 immune cells infiltrates in the TME between high- and low-risk groups in the ICGC database. E The mRNA expression level of FLVCR1 was positively correlated with the proportion of Tregs infiltration. F The mRNA expression level of EZH2 was positively correlated with the proportion of Tregs infiltration. G The mRNA expression level of EZH2 was positively correlated with the proportion of macrophage infiltration. H–J Differences in StromalScore, ImmuneScore, and TumorPurity between the high- and low-risk groups

Correlations between HRGPS and clinical characteristics

We explored the association between HRGPS and clinical characteristics. Notably, patients with liver fibrosis, more advanced tumor stage (Stage III + IV), and higher tumor grade (G3 + G4) were primarily distributed in the high-risk group (Fig. 4A–C). In addition, microvascular invasion (MVI) is an important risk factor for early recurrence and metastasis of HCC after surgery. We found many patients with MVI-positive HCC in the high-risk group (Fig. 4D), suggesting that the risk score can predict early recurrence and metastasis of HCC. The frequency of TP53 somatic mutations in patients with HCC in the high-risk group was 40.98%, significantly higher than the 14.48% in the low-risk group (Fig. 4E). Conversely, mutations in CTNNB1 showed the opposite pattern (Fig. 4F). It was demonstrated that there was a strong correlation between the TP53 mutation and the TME of HCC, wherein FOXP3 + Tregs were more prevalent, and CD8 + T cells were less invasive, leading to a downregulated immune response. Our study showed patients with TP53 mutations overexpressed immune checkpoint molecules (such as PD1, PD-L1, TIGIT, and CTLA4) (Fig. 4G). These results suggest that TP53 mutations enhanced the expression of the tumor cell surface molecule PD1 and immune cell surface molecules PD-L1, TIGIT, and CTLA4. This indicates that patients with HCC in the high-risk group are more likely to be responsive to immunotherapy. Survival analysis showed that patients with HCC with MVI-negative, TP53 wild-type, and CTNNB1 wild-type in the low-risk group had a better prognosis (Fig. 4H).

Fig. 4

The relationship between risk scores and clinical characteristics. A–C The relationship between risk score and clinical characteristics such as background of liver fibrosis, tumor stage, and tumor grade. D–F The relationship between risk score and MVI, TP53, and CTNNB1 mutations. G The relationship between immune checkpoint gene expression and TP53 mutations. H Impact of risk scores and MVI, TP53, and CTNNB1 mutations on the prognosis of patients with HCC

Both univariate and multivariate Cox analyses demonstrated that the risk score independently influenced the prognosis of patients with HCC (Table 1). Furthermore, the ROC curves indicated that the risk score had superior predictive power compared to other clinical characteristics (Figure S2A). Subsequently, we constructed a nomogram based on the risk score and tumor stage, which could effectively predict the prognosis of patients with HCC (Fig. 5A). Calibration and ROC curves showed that the nomogram accurately predicted patients with HCC survival at 1, 3, and 5 years (Fig. 5B, Figure S2B). To further assess the robustness of the risk scores in predicting the prognosis of patients with HCC, we first stratified the patients according to age (Fig. 5C, D), gender (Fig. 5E, F), MVI (Fig. 5G, H), tumor grade (Fig. 5I, J), and tumor stage (Fig. 5K, L). After stratification, OS was significantly lower in all high-risk groups except for the female subgroup.

Table 1 Results of univariate and multivariate Cox analyses of clinical characteristics and risk scores of HCC patients in the training setFig. 5

Construction of nomogram and internal validation of prognostic signature. A A nomogram was constructed based on the tumor stage and risk score. B Calibration curves showing strong predictive power of nomogram for 1-, 3-, and 5-year survival in patients with HCC in TCGA database. C–L The validity of prognostic signature was validated in different clinical subgroups, including young (C), old (D), male (E), female (F), MVI-negative (G), MVI-positive (H), low tumor grade (I), high tumor grade (J), early stage (K), and advanced stage(L)

Validation of the prognostic signature

To assess the predictive robustness of HRGPS, we computed the risk scores for two independent cohorts using the same formula as the training set. Subsequently, patients were stratified into low- and high-risk groups according to the median risk score. According to the findings, patients in the low-risk group had a better OS than those in the high-risk group (Figure S3A–S3B). However, HRGPS exhibited high AUC values in predicting the probabilities of 1-, 3-, and 5-year survival (Figure S3C–S3D), suggesting its excellent performance in predicting the prognosis of patients with HCC. In the two external independent datasets, expression levels of EZH2 and FLVCR1 were found to increase concomitantly with higher risk scores (Figure S3E–S3F).

Prediction of immunotherapy efficacy by the HRGPS

Based on the aforementioned study of immune cells, we further investigated whether histamine regulates the immune checkpoints, which are important regulators of TME, and their connection to HRGPS-based risk score. The expression levels of immune checkpoint-related genes play a crucial role in HCC immunotherapy, and our findings indicate that YTHDF1, CD86, and CD80 were significantly overexpressed in the high-risk group (Fig. 6A). Furthermore, all six immune checkpoint genes showed higher expression in the high-risk group compared to the low-risk group, as observed from our initial comparison of PDCD1, TIGIT, TIM-3, and CTLA4 expression patterns (Fig. 6B–G). These findings indicate that immunotherapy shows promise in reducing the immunosuppressive conditions observed in the high-risk group. In addition, rearrangement analysis of TCR and BCR was used to stratify and monitor the patients undergoing immunotherapy by identifying MHC-presented antigens. The cancer testicular antigen (CTA) score is an indirect measure of the tumor immunogenicity and thus indicates a potential response to immunotherapy. We found that the high-risk group exhibited higher TCR, BCR, and CTA scores than the low-risk group (Fig. 6H–J), suggesting a possible positive response to immunotherapy. Furthermore, TIDE was utilized to predict responses to immunotherapy. Interestingly, the high-risk individuals displayed lower TIDE scores (Fig. 7A). We also calculated risk scores for IMvigor210 patients. The complete/partial response (CR/PR) group had higher risk scores than the stable/progressive disease (SD/PD) group (Fig. 7B), and greater high-risk patients responded positively to immunotherapy (Fig. 7C).

Fig. 6

Prediction of immunotherapy response by prognostic signature. A Relationship between the mRNA expression levels of immune checkpoint-related genes and prognostic signature. B–G Differences in mRNA expression levels of immune checkpoint-related genes commonly used in HCC immunotherapy in high- and low-risk groups. H–J Differences in TCR richness, BCR richness, and CTA scores between the high- and low-risk groups

Fig. 7

Prediction of immunotherapy response and chemotherapeutic drugs by prognostic signature. A Differences in the distribution of TIDE scores between the high- and low-risk groups. B In the IMvigor210 dataset, patients who responded to immunotherapy had significantly higher risk scores than the low-risk group. C A higher proportion of patients in the high-risk group responded to immunotherapy compared to the low-risk group in the IMvigor210 dataset. D–I Higher risk scores were associated with lower IC50 values for anticancer drugs such as Embelin, Salubrinal, Shikonin, Lenalidomide, Doxorubicin, and Tipifarnib. J EZH2 was overexpressed in HCC tissues in clinical samples compared to paraneoplastic tissues. K FLVCR1 was overexpressed in HCC tissues in clinical samples compared to paraneoplastic tissues

Analysis of response to chemotherapy

The calculated IC50 varied significantly between the two risk populations. The fact that Embelin, Salubrinal, Shikonin, Lenalidomide, Doxorubicin, and Tipifarnib IC50 values were lower in the high-risk group than in the low-risk group (Fig. 7D–I) suggests that patients with high-risk scores may be more sensitive to these chemotherapeutic agents.

Verification of signature genes

We assessed the EZH2 and FLVCR1 mRNA expression levels in 10 pairs of HCC and paraneoplastic tissue samples. The results showed that the expression levels of EZH2 and FLVCR1 were higher in HCC samples compared to paracancerous tissues (p < 0.01) (Fig. 7J, K), suggesting that these genes may be involved in the progression of HCC.