Introduction

Hepatocellular carcinoma (HCC) is the most prevalent form of primary liver cancer. The global incidence of HCC is rising, primarily due to factors such as hepatitis virus infections, excessive alcohol consumption, and metabolic syndrome1. Recent meta-analyses estimate that there are approximately 747,000 cases of HCC worldwide, with global mortality rates ranging from 0.8–0.92%2,3.

The prevalence of HCC has shifted from regions with low to moderate sociodemographic index (SDI) to those with high SDI, with viral hepatitis identified as the primary cause of HCC in Asia and much of Africa. In contrast, alcohol-related liver disease (ALD) and non-alcoholic steatohepatitis (NASH) are the leading contributors to HCC in North America, Europe, Latin America, and Australasia4,5.

As HCC progresses, tumors often exhibit invasive and metastatic behavior, largely due to immune escape mechanisms6. This significantly worsens the prognosis of patients with HCC and presents a considerable challenge for clinicians7. Currently, treatment options for HCC include hepatectomy (the primary therapeutic approach), liver transplantation, and percutaneous ablation8,9. However, these treatments are associated with a high recurrence rate10, and liver transplantation is limited by donor availability and extended waiting periods11. Once metastasis occurs, the effectiveness of these treatments diminishes or may become unfeasible, severely compromising patient prognosis. Therefore, there is an urgent need to identify novel prognostic biomarkers and elucidate the underlying mechanisms, which could enhance our understanding of HCC pathogenesis and provide new targets for clinical diagnosis and therapy.

Polyamines, including spermine (SPM), spermidine (SPD), and putrescine12, play pivotal roles in normal growth, development, and specific pathological conditions, such as tumorigenesis13. The regulation of polyamines encompasses their synthesis, transport, and catabolism, and is intricately involved in processes like cell proliferation, apoptosis, and gene regulation14. Polyamine metabolism is not only a key metabolic pathway in the development of the immunosuppressive phenotype of tumor-associated macrophages (TAMs) but also plays a critical role in mitochondrial quality control, signal transduction, and energy metabolism regulation11,. These processes are directly linked to cancer initiation and progression. In the context of HCC, polyamine metabolism not only accelerates tumor progression but also mediates interactions between HCC cells and immune cells11. Exosomes derived from HCC can serve as diagnostic and prognostic biomarkers by modulating immune responses11. Thus, polyamine metabolism is an attractive target for HCC therapy. Polyamine metabolism-targeting drugs or analogs have demonstrated efficacy against cancer in experimental animal models and even in clinical trials11,15.

Despite these advances, the precise molecular mechanisms underlying polyamine metabolism in the prognosis of HCC remain unclear.

In this study, three prognostic genes associated with polyamine metabolism-related genes (PMRGs) in HCC were identified using transcriptome data from patients with HCC and normal control samples sourced from the GEO and UCSC Xena databases. A risk model was then constructed based on these three prognostic genes, and the expression of these markers was validated through qRT-PCR, ELISA, and in the TCGA-LIHC cohort. Further analyses, including functional enrichment, immune cell infiltration, and regulatory network exploration, were conducted. This represents investigation into the role of PMRGs in HCC, providing novel insights into their potential function and molecular mechanisms in HCC. These findings could serve as a theoretical foundation for the development of new therapeutic strategies for HCC.

Results

The 8 candidate genes significantly enriched for polyamine metabolic-related signaling pathways

Following differential expression analysis, 5,492 differentially expressed genes (DEGs, they were all protein-coding transcripts) were identified, consisting of 3,468 upregulated and 2,024 downregulated genes (Fig. 1A and B). These DEGs were intersected with 19 PMRGs, resulting in the identification of 8 candidate genes (Fig. 1C). Functional analyses were then performed to explore the biological functions and pathways associated with these candidate genes. Gene Ontology (GO) enrichment analysis revealed that the candidate genes were involved in 129 functional entries, including 91 biological processes (BPs), 8 cellular components (CCs), and 30 molecular functions (MFs), with key processes such as polyamine metabolic process and polyamine biosynthetic process (Fig. 1D). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis highlighted significant pathways such as arginine and proline metabolism and ferroptosis (Fig. 1E). The Protein-Protein Interaction (PPI) network revealed 10 pairs of interactions among the candidate genes, such as AGMAT-SAT1 and SRM-OAZ3 (Fig. 1F).

Fig. 1
figure 1

Differential expression analysis and candidate gene identification analysis. (A) A volcano plot illustrating differentially expressed genes (DEGs), highlighting top 5 up-regulated genes and top 5 down-regulated genes. (B) A heatmap displaying the expression patterns of the top 20 up-regulated and down-regulated DEGs. (C) Venn diagram between DEGs and PMRGs. (D) GO enrichment analysis of candidate genes, to elucidate their potential biological roles. (E) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of candidate genes. (F) A Protein-Protein Interaction (PPI) network of candidate genes.

Full size image

SMOX, SRM, and SAT1 were used to construct a risk model

Based on these 8 candidate genes, univariate Cox regression analysis identified 4 genes (SMOX, SRM, SAT1, and OAZ3) significantly associated with the survival of patients with HCC (Hazard Ratio (HR) ≠ 1, p < 0.05) (Fig. 2A). All 4 genes satisfied the PH hypothesis test (p > 0.05) (Table 1). Multivariate Cox regression further confirmed SMOX, SRM, and SAT1 as prognostic genes for constructing the risk model (HR ≠ 1, p < 0.2) (Fig. 2B). Among these, SAT1 acted as a protective factor, while SMOX and SRM were identified as risk factors for HCC, consistent with the univariate analysis results.

Fig. 2
figure 2

Prognostic risk model construction. (A) A forest plot of univariate Cox regression analysis. (B) A forest plot of multivariate Cox regression analysis. (C) Kaplan-Meier (K-M) survival curves comparing the survival outcomes of high-risk and low-risk groups in the training set. (D)The distribution of risk score among patients in the training set. (E) The Receiver Operating Characteristic (ROC) curve of the risk model. (F) K-M survival curves comparing the survival outcomes of high-risk and low-risk groups in the validation set. (G) A curve showing the distribution of risk score among patients in the validation set. (H) ROC curve demonstrating the predictive performance of the risk model in the validation set.

Full size image
Table 1 The results of the proportional hazards (PH) hypothesis testing for the genes (SMOX, SRM, SAT1, OAZ3) identified through univariate Cox regression analysis indicate that the null hypothesis cannot be rejected (P > 0.05).
Full size table

The risk score was calculated as follows: risk score = (−0.2865) × SMOX + (0.2319) × SRM + (−0.2882) × SAT1. Using a minimum proportion threshold of 0.45, TCGA-LIHC samples were categorized into high- and low-risk groups (high: low = 169:194). Patients in the high-risk group exhibited significantly lower survival rates on the Kaplan-Meier (K-M) survival curve (p < 0.0001) (Fig. 2C), with similar results observed in the risk curve (Fig. 2D). Furthermore, the 1-, 3-, and 5-year Area Under the Curve (AUC) values from the Receiver Operating Characteristic (ROC) curve were 0.70, 0.69, and 0.64, respectively, all exceeding 0.6, indicating good diagnostic performance of the risk model for HCC (Fig. 2E).

The risk model was validated in the GSE14520 dataset. In this cohort, patients with HCC in the high-risk group again demonstrated significantly lower survival rates, both in the K-M curve and risk curve (p = 0.016) (Fig. 2F-G). The AUC values for 1- and 5-year survival were above 0.6, and the 3-year AUC was close to 0.6 in the ROC curve, further confirming the strong diagnostic performance of the risk model (Fig. 2H).

Constructing a good nomogram

An independent prognostic analysis was performed by combining the risk score from Sect. 3.2 with clinical characteristics. Univariate Cox analysis revealed that pathologic T stage, stage, risk score, and pathologic M stage were significantly associated with the survival of patients with HCC (Fig. 3A), with p-values from the Proportional Hazards (PH) hypothesis test greater than 0.05 for all variables (Table 2). Multivariate Cox regression identified stage and risk score as independent prognostic factors for HCC (Fig. 3B).

Fig. 3
figure 3

Construction of an independent prognostic model. (A) Forest plot of univariate Cox analysis. (B) Forest plot of multivariate Cox analysis. (C) Nomogram model. (D) Correction curve of nomogram. (E) ROC curve of nomogram. (F 1–3) The 1, 3, and 5 years decision curves of nomogram.

Full size image
Table 2 The results of the proportional hazards (PH) hypothesis testing for the factors (T stage, stage, risk score, and M stage) identified through univariate Cox regression analysis indicate that the null hypothesis cannot be rejected (P > 0.05).
Full size table

A nomogram model was then constructed based on these independent prognostic factors to predict 1-, 3-, and 5-year survival rates in patients with HCC (Fig. 3C). The nomogram’s predictive ability was further assessed using various methods. The calibration curve demonstrated a slope close to 1 (Fig. 3D), and the AUC values for the nomogram model in the ROC curve at 1, 3, and 5 years were all greater than 0.7 (Fig. 3E). Additionally, Decision curve analysis (DCA) showed that the nomogram provided significant additional benefit at 1, 3, and 5 years (Fig. 3F), confirming its excellent diagnostic performance.

Next, the relationship between clinical features and the risk score was explored. Significant differences were observed between risk groups for subgroups of stage, grade, and T stage (Fig. 4A). Specifically, Wilcoxon tests revealed that the proportions of stages 2–4 in stage (p < 0.001), grade (p < 0.001), and T stage (p < 0.001) were higher in the high-risk group (Fig. 4B). Additionally, risk scores differed significantly among the subtypes of stage, grade, and T (p < 0.05), with higher risk scores associated with more advanced stages (Fig. 4C).

Fig. 4
figure 4

Correlation analysis of clinical features. (A) Heatmap showing the correlation between prognostic genes and various clinical features. (B) Comparison of subgroup counts with different clinical characteristics between high- and low-risk groups (Wilcoxon test). (C) Violin plots displaying RiskScore distributions across different clinical feature subgroups (Wilcoxon test). **** represented p < 0.0001, *** represented p < 0.001, ** represented p < 0.01, * represented p < 0.05, and ns represented p > 0.05.

Full size image

Pathway enrichment and immune cell expression differed between risk groups

Gene Set Enrichment Analysis (GSEA) was conducted to investigate biological functions across different risk groups. A total of 46 signaling pathways were significantly enriched, with the top 5 pathways visualized, including drug metabolism cytochrome P450 and retinol metabolism (Fig. 5A).

Fig. 5
figure 5

Pathway enrichment and immune infiltration analysis of different risk groups. (A) Gene Set Enrichment Analysis (GSEA) between the high and low-risk groups. (B) Heatmap displaying immune cell infiltration scores across different risk groups. (C) Box plots showing differences in immune cell infiltration levels between high- and low-risk groups. (D) Heatmap illustrating the correlation coefficients between prognostic biomarkers and differentially infiltrated immune cells. (E) Comparison of immune checkpoint gene expression levels between high- and low-risk groups. (F) The relationship between risk score and expression levels of differential immune checkpoints. (G) Boxline plot displaying differences in Tumor Immune Dysfunction and Exclusion (TIDE) scores between high- and low-risk groups. (H) Scatter plot showing the correlation between individual risk score and TIDE scores, with a regression line indicating the trend. **** represented p < 0.0001, *** represented p < 0.001, ** represented p < 0.01, * represented p < 0.05, and ns represented p > 0.05.

Full size image

Immune infiltration analysis was performed by calculating immune scores for 28 immune cell types across the different risk groups (Fig. 5B). Significant differences were found for 11 immune cell types between the risk groups. For instance, Eosinophil expression was significantly lower in the high-risk group, while the remaining 10 cell types, including activated CD4 T cells and activated dendritic cells, were elevated in the high-risk group (Fig. 5C). Correlation analysis revealed that Eosinophils exhibited the strongest positive correlation with SAT1 (r = 0.41, p < 0.001) and the strongest negative correlation with SRM (r = −0.15, p < 0.001) (Fig. 5D). The expression of 8 immune checkpoints was found to be downregulated in the high-risk group (Fig. 5E), with a significant positive correlation to the risk score. LGALS9 showed the highest correlation with the risk score (r = 0.43, p < 0.001) (Fig. 5F). Furthermore, the Tumor Immune Dysfunction and Exclusion (TIDE) score was significantly higher in the high-risk group (Fig. 5G), suggesting increased immune evasion potential, which could be linked to a poorer response to immune checkpoint inhibitor (ICI) therapy. Additionally, the TIDE score exhibited a weak positive correlation with the risk score (r = 0.37, p = 6.4 × 10–3), indicating that a higher risk score may be associated with a less favorable response to immunotherapy (Fig. 5H).

TF-mRNA-miRNA and miRNA-SNP-mRNA networks of prognostic genes were generated

In this study, the regulatory mechanisms of prognostic genes were explored. Initially, 90 miRNAs predicted by the starBase database were intersected with 59 miRNAs predicted by the miRTarBase database, resulting in 8 targeted miRNAs for the 3 prognostic genes (Fig. 6A). Additionally, 66 Transcription factors (TFs) corresponding to the prognostic genes were predicted. Subsequently, a TF-mRNA-miRNA network was constructed, consisting of 85 relationship pairs, such as POU2F2-SAT1-hsa-miR-335-5p, TAL1-SMOX-hsa-let-7d-5p, and others (Fig. 6B). Furthermore, 11 Single Nucleotide Polymorphisms (SNPs) corresponding to the prognostic genes were predicted using the miRNASNP database. The miRNA-SNP-mRNA network was composed of 19 nodes and 22 relationship pairs, including pairs like hsa-miR-186-5p-rs1445781653-SAT1 and hsa-miR-760-rs569643928-SMOX (Fig. 6C).

Fig. 6
figure 6

Construction of prognostic gene regulatory network. (A) A Venn diagram showing the intersection of miRNAs identified through different databases. (B) LncRNA-miRNA-mRNA network. (C) MiRNAs-single nucleotide polymorphisms (SNPs-mRNAs network.

Full size image

HLA family genes, tumor immune cycle, and IC50 of sensitive drugs were different among different risk groups

The expression differences of Human leukocyte antigen (HLA) family genes across risk groups were also examined. The results revealed that the expression of 15 HLA family genes, such as HLA-A, HLA-DMA, and HLA-DMB, was significantly upregulated in the high-risk group (Fig. 7A). Additionally, tumor immune cycle analysis identified 10 categories of the immune cycle, including Step 1, Step 4.Th22.cell.recruiting, and Step 4.Neutrophil.recruiting, with significant differences between risk groups (Fig. 7B). Furthermore, 117 drugs exhibited significant differences in Half maximal inhibitory concentration (IC50) values between the risk groups, indicating varying therapeutic effects in different stages of HCC progression. Among the top drugs, Doramapimod_1042 was more effective in the high-risk group, while five other drugs, such as Sepantronium bromide_1941 and Paclitaxel_1080, showed greater efficacy in the low-risk group (Fig. 7C).

Fig. 7
figure 7

HLA family gene expression, tumor immune circulation and drug sensitivity analysis. (A) Boxline plot showing the differences in HLA family gene expression between high- and low-risk groups. (B) Comparison of Tumor Immune Cycle (TIP) differences between high- and low-risk groups. (C) Boxline plot illustrating the differences in IC50 values of chemotherapeutic drugs between high- and low-risk groups. **** represented p < 0.0001, *** represented p < 0.001, ** represented p < 0.01, * represented p < 0.05, and ns represented p > 0.05.

Full size image

Validation of expression levels of SMOX, SRM, and SAT1

In the TCGA-LIHC dataset, the expression of prognostic genes was further investigated. SMOX and SRM were significantly overexpressed in the HCC group, while SAT1 displayed the opposite pattern (Fig. 8A). In the GSE14520 dataset, SRM showed significantly higher expression in HCC samples (p < 0.0001), while SAT1 was significantly more highly expressed in control samples (p < 0.0001), consistent with the findings from TCGA-LIHC. However, the expression of SMOX did not differ significantly between HCC and control samples (p > 0.05) (Supplementary FigureS1).

Fig. 8
figure 8

Verification of the expression levels of three prognostic genes. (A) The expression levels of prognostic genes in the TCGA-LIHC dataset. (B) Prognostic gene expression in quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR) experiments. (C) Prognostic gene expression in serum of HCC patients and normal people by Elisa experiment. (D) The expression of prognostic genes in the supernatant of HepG2 cells in the control group and the intervention group was detected by Elisa. (E) The expression of prognostic genes in the supernatant of Hep3B cells in the control group and the intervention group was detected by Elisa. **** represented p < 0.0001, *** represented p < 0.001, ** represented p < 0.01, * represented p < 0.05, and ns represented p > 0.05.

Full size image

Similarly, Quantitative real-time polymerase chain reaction (qRT-PCR) experiments showed that SRM and SAT1 were expressed at higher levels in the HCC group, while SMOX expression was significantly higher in control samples (Fig. 8B, Supplementary Figure S2). Variability across datasets and qRT-PCR results can be attributed to differences in sample types and experimental methods. Further validation was performed by measuring the serum levels of SMOX, SRM, and SAT1 in patients with HCC and healthy individuals, which also aligned with the qRT-PCR findings (Fig. 8C). Finally, HepG2 and Hep3B cells were treated with Doramapimod, and changes in the levels of SMOX, SRM, and SAT1 in the cell supernatant were consistent with the above results (Fig. 8D-E).

Discussion

HCC is the sixth most common cancer in the world, with a high mortality rate16. From 2000 to 2016, the 5-year survival rate of HCC patients in the United States was only 18%17. Recent studies have found that polyamine metabolism plays a vital role in changing the tumor immune microenvironment18. The exploration for potential regulatory mechanism of polyamine metabolism can provide an important means to understand the process of tumor evolution19. Therefore, the polyamine metabolic pathway is a new way for anti-HCC treatment.

In this study, prognostic gene related to PMRGs (SMOX, SRM, SAT1) were identified for constructing risk model, and a nomogram was constructed to predict the prognosis of patients. The identification of three prognostic genes was based on a series of analyses, making these genes relatively more reliable compared to other markers, such as AFP.

At the same time, we used TCGA-LIHC, qRT-PCR and ELISA experiments to verify the expression levels of SMOX, SRM, SAT1 in HCC patients, and explored their prognostic value in HCC patients. We also analyzed the correlation between polyamine metabolism-related subgroups and immune cells and immune checkpoints in HCC patients, and screened IC50-sensitive drugs in different subgroups.

Spermine oxidase (SMOX) is an inducible enzyme in the polyamine catabolism pathway, which specifically oxidizes spermine to spermidine and produces hydrogen peroxide and 3-aminopropanal20,21. Hu et al.found that SMOX is highly expressed in HCC cell lines and clinical samples, and can promote the tumor growth of HCC22. Targeted SMOX can reduce the progression of liver cancer22. At the same time, studies have found that SMOX can promote the inflammatory response and fibrosis of hepatocytes through the β-catenin pathway23, suggesting that SMOX may indirectly promote the occurrence and development of HCC by affecting the inflammatory response and fibrosis process. In addition, SMOX catalyzes the oxidation of polyamines to produce reactive oxygen species (ROS)22. Excessive accumulation of ROS may lead to DNA damage, cell apoptosis and gene mutation, thereby promoting the occurrence and development of HCC. In conclusion, SMOX may synergistically drive the occurrence and development of HCC by mediating ROS production, activating β-catenin pathway and promoting fibrosis, and targeting SMOX may become a potential therapeutic strategy.

Spermidine synthase (SRM, also known as SPDS) is an enzyme that is widely present in organisms and is responsible for catalyzing the synthesis of spermidine24. Studies have shown that spermidine secreted by tumor cells can inhibit the aggregation of T cell receptor (TCR), thereby inhibiting the activation of CD8+ T cells in tumors25. This suggests that SRM may inhibit the activation of CD8+ T cells in HCC, thereby affecting the progression of HCC. In addition, SRM is enriched in the extracellular vesicles of HCC cells under hypoxic conditions. It is an effective regulator of polyamines and energy metabolism in tumor-associated macrophages, which can promote the progression of liver cancer26. In other words, SRM may inhibit CD8+ T cell activation by regulating spermidine synthesis and mediate tumor-macrophage metabolic interaction, thereby promoting HCC progression.

Spermidine/spermine N1-acetyltransferase 1 (SAT1) encodes a rate-limiting enzyme in polyamine metabolism27,28. SAT1 regulates intracellular polyamine levels by acetylating spermine and spermidine, thereby affecting various biological processes of cells28,29. SAT1 can regulate the growth and metastasis of HCC cells through AKT/β-catenin signaling pathway29. SAT1 can be used as a gene-specific transcriptional regulator to drive the invasiveness of brain tumors28, suggesting that SAT1 may also drive the invasiveness of HCC through the same mechanism in HCC. In addition, SAT1 activation is associated with p53-mediated ferroptosis30. Ferroptosis is an iron-dependent way of cell death. In low-grade gliomas, high expression of SAT1 is associated with poor prognosis and immune cell infiltration31. This means that SAT1 may affect the progression of HCC by affecting ferroptosis. In summary, SAT1 may affect the growth, metastasis and invasion of HCC through polyamine metabolism, signaling pathways and ferroptosis.

In this study, GSEA found that pathways such as drug metabolism cytochrome P450 and Steroid hormone biosynthesis were significantly enriched in the low-risk group of HCC. Cytochrome P450 enzymes play a vital role in drug metabolism, especially in the liver, which is responsible for the oxidative metabolism and clearance of many clinical drugs. These enzymes affect the body’s response to drugs, safety, bioavailability and toxicity32. Cytochrome P450 3A4, a subtype of cytochrome P450, is involved in the metabolism of various chemotherapeutic drugs. Changes in the expression or activity of cytochrome P450 3A4 in tumor cells can affect the efficacy of chemotherapeutic drugs33. In other words, by affecting the metabolism of drugs, cytochrome P450 may affect the effective concentration of HCC therapeutic drugs, thus affecting the therapeutic effect. Steroid hormones, such as estrogen and androgen, are important regulators of cell cycle progression34. These hormones affect the transcription of target genes by binding to intracellular receptors, thereby regulating cell proliferation and differentiation35,36. Therefore, in HCC, abnormal steroid hormone signal transduction may lead to excessive proliferation and abnormal differentiation of cancer cells. In summary, the enrichment of drug metabolism cytochrome P450 and Steroid hormone biosynthesis pathways in the low-risk group of HCC may jointly regulate the treatment response and malignant progression of HCC by affecting drug metabolism and cell proliferation.

In summary, this study identified three prognostic genes associated with polyamine metabolism in HCC, and constructed a risk model based on these prognostic genes, and explored their functions and potential regulatory mechanisms. It provides a new theoretical reference for the in-depth study of the pathogenesis and new treatment strategies of HCC. In the future, we will focus on studying these prognostic genes, further elucidating their specific mechanisms in HCC, and optimizing the treatment strategies of HCC.

However, this study has some limitations. First, the clinical application of the prognostic genes we screened requires more sample data support. Secondly, due to differences in the datasets and gene expression profiles, the validation results of the risk model in the GSE14520 dataset exhibited lower AUC values. Moving forward, we plan to expand the sample size and perform additional validations using diverse datasets to enhance the robustness of the model. Thirdly, given the limited number of clinical samples used in the qRT-PCR and Elisa validation, it is essential to expand the sample size to enhance the robustness and generalizability of our findings. Meanwhile, a large number of experimental mechanisms (such as gene knockout, etc.) need to be further studied to explore the specific functions and mechanisms of these genes related to the prognosis of HCC. We will continue to pay attention to the role of these genes.

Methods

Data acquisition

The HCC-related dataset TCGA-LIHC with survival information was selected from the UCSC Xena database (https://xenabrowser.net/datapages/) for analysis. The TCGA-LIHC dataset included 419 tissue samples (Adjacent control: HCC tumor = 50:369). Of these, 363 samples had available overall survival (OS) data. To process the TCGA-LIHC Counts matrix, several normalization and filtering steps were performed to ensure high-quality data and reliable results. Initially, a logarithmic transformation was applied to stabilize variance and mitigate skewness. The data were then normalized using Counts Per Million (CPM), and low-expression genes were excluded, retaining only those with a CPM value exceeding 0.5 in at least half of the samples. Additionally, genes with expression levels above zero in more than 80% of the samples were selected. These steps refined the TCGA-LIHC Counts matrix, ensuring robust and accurate analysis.

Similarly, the HCC-related dataset GSE14520 (GPL3921) with survival information in the GEO database (https://www.ncbi.nlm.nih.gov/gds) was selected to validate the risk model. GSE14520 dataset included 445 tissue samples (Control: HCC tumor = 220:225), with 221 HCC samples possessing OS data. Additionally, according to the references37, 19 PMRGs (Supplementary Table 1) were derived from the GOBP_POLYAMINE_METABOLIC_PROCESS functional entry in the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb) (accessed on November 27, 2023).

Identification and function exploration of candidate genes

To identify genes with significant expression changes under different conditions, differential expression analysis was performed on the TCGA-LIHC dataset using the DESeq2 package (v 1.36.0). DEGs between the HCC and control groups were identified with the following criteria: |log2 fold change (FC)| > 0.5 and p.adjust < 0.0538. Candidate genes were then identified by intersecting DEGs with PMRGs using the ggvenn package (v 0.1.9)39. GO and KEGG enrichment analyses were conducted using the clusterProfiler package (v 4.7.1.001) to further explore the molecular functions of the candidate genes (p < 0.05)40,41,42,43. Subsequently, a PPI network for the candidate genes was constructed using the STRING database (https://string-db.org/), which provided insights into their protein-level interactions (interaction confidence = 0.4).

Risk model construction and verification

Univariate Cox regression analysis was conducted on the candidate genes using the survival package (v 3.4-0) in the TCGA-LIHC dataset to identify genes associated with survival (HR ≠ 1, p < 0.05)44. To screen for prognostic genes for the risk model, genes that met the PH hypothesis test (p > 0.5) were subjected to multivariate Cox regression analysis (HR ≠ 1, p < 0.2). The resulting risk scores were ranked from low to high, and the optimal threshold was determined using the surv_cutpoint function in the survminer package (v 0.4.9) to categorize patients with HCC into high and low-risk groups45. K-M survival curves and risk curves were plotted to visualize survival differences between the groups. ROC curves were also generated using the survival ROC package (v 1.0.3) to assess the performance of the risk model46. To further evaluate the generalizability of the risk model, validation was performed using the GSE14520 dataset. A similar risk model was constructed for GSE14520, and K-M curves, risk curves, and ROC curves were plotted to evaluate the model’s effectiveness.

Construction of the independent prognostic model

To enhance the accuracy of the prognostic model, clinical features (age, stage, grade, pathologic M, N, T, and gender) along with the risk score were incorporated into the univariate Cox analysis (HR ≠ 1, p < 0.05). Notably, in the TNM staging system for HCC, the T (Tumor), N (Nodes), and M (Metastasis) categories are critical indicators used to assess tumor progression and guide staging. Following this, the PH hypothesis test was applied to univariate Cox results (p > 0.5). Independent prognostic factors (HR ≠ 1, p < 0.2) were then identified using multivariate Cox regression analysis. A nomogram was constructed for independent prognostic factors based on the rms package (v 6.3-0)47. The predictive power and validity of the nomogram were further evaluated through calibration curves, ROC curves, and DCA. Additionally, relationships between risk score, prognostic genes, and clinical features in TCGA-LIHC were examined by generating heat maps using the pheatmap package (v 1.0.12)48. Differences in risk scores across various clinical features were analyzed using the Wilcoxon rank sum test.

Functional enrichment in high- and low-risk groups

To investigate relevant signaling pathways and potential biological mechanisms between the two risk groups, differential analysis was performed using the DESeq2 package (v 1.36.0), with log2FC serving as the sorting criterion. GSEA was then conducted based on the KEGG pathway gene set (c2.cp.kegg.v7.4.symbols.gmt) from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb), using the clusterProfiler package (v 4.7.1.001). Enrichment was considered significant with a normalized enrichment score (|NES| > 1, p-adjust < 0.05)40.

Immune-related analysis

The scores of 28 immune cell types for each sample were calculated using the single sample GESA (ssGSEA) algorithm from the GSVA package (v 1.44.5)49, and the results were visualized using heat maps generated by the ggplot2 package (v 3.3.6)50. Differences in immune cell distributions and 12 immune checkpoint molecules51 between the high- and low-risk groups were assessed using the Wilcoxon rank-sum test. Spearman correlation analysis was conducted between prognostic genes and differential immune cells as well as differential immune checkpoints using the psych package (v 2.2.9)52. To explore the association between risk score and tumor immune escape, TIDE scores for HCC samples were calculated using the TIDE website, with differences in TIDE scores across risk groups tested by the Wilcoxon rank-sum test. Additionally, Spearman correlation was applied to investigate the relationship between risk score and TIDE score.

Construction of regulatory network for prognostic gene

A regulatory network was constructed to explore the potential regulatory mechanisms of the prognostic genes. Initially, miRNAs targeting the prognostic genes were predicted using the starBase and miRTarBase databases. The results from these two databases were intersected to identify common target miRNAs. TFs for the prognostic genes were predicted using the ChEA3 database (https://maayanlab.cloud/chea3/), with TFs supported by ChIP-seq data from the ENCODE database (https://www.encodeproject.org/) being identified. A TF-mRNA-miRNA network was then constructed. To further investigate the impact of Single Nucleotide Polymorphisms (SNPs) on prognostic gene expression, SNPs within the seed regions of target miRNAs were predicted using the miRNASNP database (http://www.biobigdata.com/miRNASNP/). The TF-mRNA-miRNA and miRNA-SNP-mRNA networks were visualized and analyzed using Cytoscape (v 3.8.2)53.

HLA family gene, tumor immune cycle, and drug sensitivity analysis

The HLA genes, which are crucial for initiating cellular and humoral immunity, were compared between high- and low-risk groups in HCC samples using the Wilcoxon rank-sum test54. Differences in the scores of the 23 subcategories of the seven-step tumor immune cycle55 between risk groups were also visualized and analyzed using the Wilcoxon test. To assess potential therapeutic options, the IC50 for various targeted agents in the GDSC database (https://www.cancerrxgene.org/) was calculated using the oncoPredict package (v 0.2)56, and IC50 differences were compared between risk groups using the Wilcoxon test.

Expression verification of prognostic genes

To validate the expression trends of the prognostic genes, expression differences were analyzed in the TCGA-LIHC dataset using the Wilcoxon rank-sum test. Additionally, the expression of prognostic genes was also evaluated in the GSE14520 dataset (p < 0.05).

qRT-PCR

This study utilized five paired tissue samples from patients at Zhejiang Provincial People’s Hospital for qRT-PCR analysis, including five HCC samples and five control samples. The study received approval from the Ethics Committee of Zhejiang Provincial People’s Hospital (No. QT2024070), and informed consent was obtained from all patients.

Total RNA from frozen HCC and control samples was extracted using TRIzol reagent (Ambion, USA) (Supplementary Table 2), followed by RNA quantification using the NanoPhotometer N50. Equal amounts of mRNA were reverse transcribed into cDNA using the SureScript First-Strand cDNA Synthesis Kit. The resulting cDNA was diluted 5–20 times with RNase/DNase-free ddH2O (Sevier). qPCR was performed using the 2 x Universal Blue SYBR Green qPCR Master Mix Kit under the following conditions: initial denaturation at 95 °C for 1 min, followed by 40 cycles of denaturation at 95 °C for 20 s, annealing at 55 °C for 20 s, and extension at 72 °C for 30 s. Gene expression levels were calculated using the 2–ΔΔCT method. Primer sequences are provided in Supplementary Table 3.

cell culture

HepG2 and Hep3B human hepatoma cell lines were purchased from Wuhan Pricella Biotechnology Co., Ltd. Both cell types were cultured in Minimum Essential Medium (MEM) supplemented with 10% fetal bovine serum (FBS, Gibco, USA) and maintained in a humidified atmosphere at 37 °C with 5% CO2. CCK8 experiment was used to explore Doramapimod, and it was found that 100 nM Doramapimod had a significant inhibitory effect on proliferation at 24 h. This concentration and time were selected for subsequent studies (Supplementary Figure S3-S4). The cells were divided into two groups: a control group, which received only medium, and an intervention group, which was treated with 100 nM Doramapimod (HY-10320, MedChemExpress) for 24 h.

Enzyme-linked immunosorbent assay

Serum and cell supernatant levels of SMOX (ED-201717), SRM (ED-17237), and SAT1 (ED-15124) were measured using enzyme-linked immunosorbent assay (ELISA) kits according to the manufacturer’s instructions (AMOY LUNCHANGSHUOBIOTECH, CO., LTD, Xiamen, China).

Statistical analysis

Bioinformatics analyses were conducted using R software (v 4.2.3). Differences were considered statistically significant at p < 0.05.