Using the WGCNA algorithm, these genes were assigned to different modules by clustering dendrograms, and finally, we obtained three modules, among which the turquoise module consisting of 86 senescence genes showed the highest correlation and significance, therefore, we selected the senescence genes in the turquoise module for subsequent analysis (Fig. S1).

Construction of Senescence-Related Gene Signature

First, we took the intersection of differential expressed genes (DEGs), prognostic genes and turquoise module genes and obtained 33 intersecting genes. Using these 33 genes to construct a signature of senescence-related genes. First, 342 samples were randomly divided 1:1 into training set and Test set. These 33 intersecting genes were analyzed by univariate Cox regression to derive genes with prognostic features in the training set for the next step of analysis. Then, the signature genes and their regression coefficients were obtained by Lasso and multivariate cox regression analysis (Fig. 3A-D). Finally, we obtained 5 signature genes: CBX2, CDKN2B, ETS2, HMGA1 and UBE2S. Therefore, we created a risk score system based on the signature genes and regression coefficients to calculate the risk score for each sample. In each data set, the sample was divided into high-risk and low-risk groups based on the median value of the risk scores. The risk score is calculated as follows:

Fig. 3

LASSO and cox regression analysis to identify signature genes in the training set. A We intersected the differentially expressed genes (DEGs), prognostic genes and modular genes obtained from WGCNA analysis and obtained 33 shared genes. B In the training set, we performed Univariate cox regression analysis on these 33 shared genes to screen for prognosis-related genes. C Cross-validation of the LASSO regression. D Coefficient value of prognostic genes

$$\mathrm= (\mathrm2*0.385) + (\mathrm2\mathrm*0.417) + (\mathrm2*-0.322) + (\mathrm1*0.197) + (\mathrm2\mathrm*0.301)$$

The expression of signature genes was demonstrated with a heat map (set 1) (Fig. 4A), and by analysis it was found that patients in the high-risk group had a lower survival rate (p < 0.05) (Fig. 4B). As the risk score gets higher, the survival time gets shorter and the number of deaths gets higher (Fig. 4C). The AUC values at 1, 2 and 3 years were 0.867, 0.755 and 0.711, respectively (Fig. 4D).

Fig. 4

Prognostic signature based on 5 senescence-related genes. A The expression of 5 senescence-related genes in the TCGA train set. B Survival analysis between the high-risk and low-risk groups in the TCGA train set. C Distribution of risk scores and survival outcomes. D Receiver operating characteristic curve (ROC) of risk score in the TCGA train set

Expression Validation of Signature Genes

First, the expression levels of the signature genes in cells (normal liver cells vs hepatocellular carcinoma cells) and tissues (hepatocellular carcinoma tissues vs paraneoplastic tissues) were verified by q-PCR assay (Fig. 5A-J).

Fig. 5

Validation of differential expression of 5 signature genes in cells and tissues by q-PCR. A-E Differential mRNA expression of 5 signature genes in 7702, LM3, 97H, HepG2 and 7721 cells. The results showed that the expression of CBX2, CDKN2B, HMGA1 and UBE2S in HCC cells was significantly higher than that in normal liver cells. In contrast, the expression of ETS2 was higher in normal liver cell. F-J Differential mRNA expression of 5 signature genes in HCC tissue and paraneoplastic tissue (30 pairs), The expression of CBX2, CDKN2B, HMGA1 and UBE2S were significantly higher in HCC and paraneoplastic tissues. In contrast, ETS2 was lower expressed in HCC tissues

Then, we examined the differences in expression of the signature genes in eight pairs of hepatocellular carcinoma and paraneoplastic tissues by Western blot (Fig. 6A). At the cellular level, we compared the protein expression differences of the signature genes between normal liver cell (7702) and HCC cells (LM3, 97H, HepG2, 7721 and Hu7) (Fig. 6B).

Fig. 6

Western Blot showed the protein expression difference of 5 signature genes in HCC tissues and cells. A Differential protein expression of 5 signature genes in 8 pairs of HCC tissues and paraneoplastic tissues. B Differential protein expression of 5 signature genes in normal liver cell (7702) and HCC cells (LM3, 97H, HepG2, 7721 and Hu7)

Finally, we carried out an immunohistochemical experiment. We showed representative images of five genes, and then compared the expression differences between HCC tissues and paraneoplastic tissues using relative optical density scores (Fig. 7).

Fig. 7

Staining images of five signature genes in HCC tissues and paraneoplastic tissues. Relative optical density scores were used to compare the differences between the two groups. A CBX2; B CDKN2B; C ETS2; D HMGA1; E UBE2S

Molecular Interaction Networks of Signature Genes and Signature Validation.

We used the GeneMANIA (http://genemania.org/) online website to analyze the molecular interaction network between the five signature genes and found that the functions of these five genes and interacting genes (e.g. ID1, ERF, CDK4, PIAS2, CDK6,CDKN2C, etc.) were mainly related to regulation of cellular senescence, cyclin-dependent protein serine/threonine kinase regulator activity, regulation of G1/S transition of mitotic cell cycle, ubiquitin ligase complex, negative regulation of cell cycle phase transition, protein kinase inhibitor activity and nuclear ubiquitin ligase complex and CDKN2B seems to be enriched with even more features (Figure S2).

In the test group (set2), survival was worse in the high-risk group (p < 0.05). The AUC values at 1, 2 and 3 years were 0.745, 0.734 and 0.719, respectively (Figure S3A-C).

In the total TCGA set (set3), to validate the extent of classification between risk groups, we used t-SNE and PCA downscaling, and we found that the samples between risk groups could be well differentiated between HCC patients (Figure S4A-B). Again, the low-risk group had better survival (p < 0.05). The AUC values at 1, 2, and 3 years were 0.810, 0.748, and 0.719, respectively (Figure S3D-F).

The signature was next tested using ICGC data, and we also used t-distributed Stochastic Neighbor Embedding (t-SNE) and Principal Component Analysis (PCA) downscaling analysis, and found that the samples between risk groups also distinguished HCC patients well (Figure S4C-D), In addition, the survival analysis was consistent with the results of the other validation sets (Figure S3G-I).

In addition, to verify the accuracy and reliability of our signature, our signature was compared with four signatures from previous studies, and it was found that the consistency index (C-index) of our signature was higher. Meanwhile, we constructed signatures with our data using the signature genes of four other studies, and then obtained ROC curves for 1, 2, and 3 years, and found that the predictive power of our signature has higher accuracy (Figure S5).

Independent Prognostic Analysis

Both univariate and multifactorial analyses found that Stage staging and riskScore were significantly associated with patients (P < 0.05) (Fig. 8A, B). For further evaluation for individual patients, we simplified the statistical prediction signature using Nomograms. The calibration chart also shows good accuracy (Fig. 8C, D).

Fig. 8

Independent prognostic analysis and Construction and validation of the Nomograms. A Univariate cox regression analysis for the TCGA cohort. B Multivariate cox regression analysis for the TCGA cohort. C Construction of the Nomograms. D The calibration curves displayed the accuracy of the nomogram in the 1-, 2- and 3 years

Clinical Correlation Analysis

We first presented the signature genes with clinically relevant indicators in heat map form, and found significant differences in T-stage, Stage staging and Grade staging between risk groups (P < 0.05) (Fig. 9A) (Figure S6). We present the clinical characteristics between risk groups in the form of box plots (Figure S7). Also, we further validated the reliability of the signature by analyzing the survival rates of patients at different clinical stages, which were low in the high-risk group in all stages and grades (p < 0.05) (Fig. 9B-E).

Fig. 9

Relationship between risk score and clinical characteristics. A The correlations between the senescence-related genes and clinicopathologic characters of the high-risk group and low-risk group were shown as a heatmap. B Survival analysis between the high-risk and low-risk groups in the patients with Stage I-II. C Survival analysis between the high-risk and low-risk groups in the patients with Stage III-IV. D Survival analysis between the high-risk and low-risk groups in the patients with Grade I-II. E Survival analysis between the high-risk and low-risk groups in the patients with Grade III-IV. *P < 0.05, **P < 0.01, ***P < 0.001

GO and KEGG Enrichment Analysis

GO enrichment analysis suggested that the high risk group was mainly associated with HUMORAL_IMMUNE_RESPONSE_MEDIATED_BY_CIRCULATING_IMMUNE, PHAGOCYTOSIS_RECOGNITION RECOGNITION and IMMUNOGLOBULIN_COMPLEX. The low-risk group is mainly associated with ALPHA_AMINO_ACID_CATABOLIC_PROCESS, CELLULAR_AMINO_ACID_CATABOLIC_PROCESS and FATTY_ACID_BETA_OXIDATION (Fig. 10A,B).

Fig. 10

GO and KEGG enrichment analysis of senescence-related genes. A GO enrichment analysis of high-risk group. B GO enrichment analysis of low-risk group. C KEGG enrichment analysis of high-risk group. D KEGG enrichment analysis of low-risk group

KEGG enrichment analysis suggested that the high-risk group was mainly associated with CELL_CYCLE, DNA_REPLICATION and ECM_RECEPTOR_INTERACTION. The low-risk group is mainly associated with COMPLEMENT_AND_COAGULATION_CASCADES, DRUG_METABOLISM_CYTOCHROME_P450 and FATTY_ACID_METABOLISM (Fig. 10C,D).

Immune Cell Infiltration Analysis

In the high-risk group, patients had higher levels of T-cell follicular helpers, T-cell regulation (Tregs), T-cell CD4 memory activation, and macrophage M0, B-cell memory (P < 0.05). However, the level of resting T-cell CD4 memory, monocytes, macrophage M1 and mast cells were lower (P < 0.05) (Table S4) (Fig. 11A). Also, Analysis of immune-related functions between the risk groups revealed significant differences in type II_IFN_response, MHC_class_I and type I_IFN_response between the two groups (Fig. 11B). We also show immune cells that have significant correlation with signature genes (Figure S8).

Fig. 11

Immune infiltration analysis and Analysis of the tumor mutation burden (TMB) and Microsatellite instability (MSI). A The comparison of CIBERSORT scores derived from 22 different immune cells. B Analysis of immune-related pathways between high-risk group and low-risk group. C The expression of TMB between high-risk group and low-risk group. D Survival analysis between the high-TMB and low-TMB groups. E Survival analysis between the risk score and TMB levels. F The expression of MSI between high-risk group and low-risk group

Tumor Microenvironment Analysis

Tumor microenvironment analysis showed significant differences in TMB levels between risk groups, with higher risk groups having higher levels of TMB expression, while higher TMB levels were associated with lower survival rates (P < 0.05) (Fig. 11C, D). The TMB levels between different risk groups affected the survival rate of patients (P < 0.05) (Fig. 11E). Microsatellite instability analysis showed lower MSI levels in the low-risk group (P < 0.05) (Fig. 11F).

Immune Checkpoint and Drug Sensitivity Analysis

Differential analysis results showed that HDAC2, PD-1, CTLA4, CD86, HHLA2, SOAT1, ICOS, CD40, CD27, CD28, IDO1, CDK1, CD276 and MMP9 were more expressed in the high-risk group (Fig. 12A). The results of correlation analysis showed significant correlations between 14 immune checkpoints and risk scores (Fig. 12B). In addition, the high-risk group had lower TIDE scores, less likelihood of immune escape, and better efficacy during immunotherapy (Fig. 12C). In recent years, tumor progression, metastasis, recurrence and tumor resistance to cytotoxic therapy play a key role, We analyzed the correlation between risk score and tumor stem cells and found that the higher the risk score, the higher the tumor stem cell score (P < 0.05) (Fig. 12D).

Fig. 12

Evaluation of immune checkpoint profiles and immunotherapy between risk groups. A The expression of immune checkpoints between high-risk group and low-risk group. B Correlation analysis between immune checkpoints and risk scores. C Comparison of the scores of TIDE between the high and low risk group. D The relationship between risk score and RNAss. E The relationship between risk score and sorafenib senstivity (IC50). F Comparison of the sorafenib senstivity (IC50) between the high and low risk group

Meanwhile, In order to analyse the differences in immunotherapy between the high and low risk groups, we used the "pRRophetic" R package to predict the gene expression and drug sensitivity of the cell lines. The analysis of riskScore and drug sensitivity showed that the higher the risk score, the higher the sensitivity to sorafenib and the better the treatment outcome is likely to be (Fig. 12E). And the high-risk group was more sensitive to sorafenib treatment (Fig. 12F).