Cellular senescence throws new insights into patient classification and pharmacological interventions for clinical management of hepatocellular carcinoma

Abstract

BACKGROUND

Cellular senescence, a state of stable growth arrest, is intertwined with human cancers. However, characterization of cellular senescence-associated phenotypes in hepatocellular carcinoma (HCC) remains unexplored.

AIM

To address this issue, we delineated cellular senescence landscape across HCC.

METHODS

We enrolled two HCC datasets, TCGA-LIHC and International Cancer Genome Consortium (ICGC). Unsupervised clustering was executed to probe tumor heterogeneity based upon cellular senescence genes. Least absolute shrinkage and selection operator algorithm were utilized to define a cellular senescence-relevant scoring system. TRNP1 expression was measured in HCCs and normal tissues through immunohistochemistry, immunoblotting and quantitative real-time polymerase chain reaction. The influence of TMF-regulated nuclear protein (TRNP)1 on HCC senescence and growth was proven via a series of experiments.

RESULTS

TCGA-LIHC patients were classified as three cellular senescence subtypes, named C1–3. The robustness and reproducibility of these subtypes were proven in the ICGC cohort. C2 had the worst overall survival, C1 the next, and C3 the best. C2 presented the highest levels of immune checkpoints, abundance of immune cells, and immunogenetic indicators. Thus, C2 might possibly respond to immunotherapy. C2 had the lowest somatic mutation rate, while C1 presented the highest copy number variations. A cellular senescence-relevant gene signature was generated, which can predict patient survival, and chemo- or immunotherapeutic response. Experimentally, it was proven that TRNP1 presented the remarkable upregulation in HCCs. TRNP1 knockdown induced apoptosis and senescence of HCC cells and attenuated tumor growth.

CONCLUSION

These findings provide a systematic framework for assessing cellular senescence in HCC, which decode the tumor heterogeneity and tailor the pharmacological interventions to improve clinical management.

Key Words: Cellular senescence, Hepatocellular carcinoma, Prognosis, Subtypes, Tumor microenvironment, Gene signature, Pharmacological interventions

Core Tip: Cellular senescence, a state of stable growth arrest, is implicated in human cancers. Nevertheless, characterization of cellular senescence-associated phenotypes in hepatocellular carcinoma (HCC) is still indistinct. Here, we proposed a novel cellular senescence-based classification for HCC and identified TRNP1 as a novel therapeutic target.

INTRODUCTION

Cellular senescence is defined as an irreversible cessation of cellular division of cells with normal proliferation[1]. Human cells age due to progressive shortening of telomeres following cellular division, stress, oncogenes, etc[2]. Numerous genes have been implicated in cellular senescence as biomarkers and causal drivers[3,4]. Cellular senescence is a double-edged sword for cancer and its treatment[5,6]. The growth arrest and immunomodulatory features linked with senescence possess powerful antimalignant roles[7]. In addition, senescence bypass and secretory phenotype correlate to tumor progression and recurrence[8]. Research has unveiled the great potential for antiaging interventions as a novel antitumor strategy[9]. Nonetheless, the heterogeneity of senescence-related features makes the definition and targeting of treatment-induced senescent cells challenging[10].

Hepatocellular carcinoma (HCC) is a poorly managed malignancy with high mortality due to the lack of response to classical chemotherapy agents (doxorubicin, cisplatin, etc.) and targeted agents in the early stage[11]. For late-stage HCCs, single sorafenib or combination therapy remains the mainstay in first-line therapy, which improves overall survival by 3 mo[12]. The modest therapeutic success is largely attributable to sorafenib resistance[13]. Immunotherapy with checkpoint inhibitors (anti-PD-1/PD-L1) has displayed potent anti-HCC activity in a subset of patients[14]. The main unmet challenge in HCC immunotherapy is to discover and verify predictive biomarkers[15]. Accumulated evidence demonstrates that inducing tumor cells into senescence represents a potential anti-HCC therapy[16]. In HCCs, cellular senescence is primarily controlled by p53-dependent or -independent mechanisms[17]. Paradis et al[18] investigated replicative senescence in normal liver, chronic hepatitis C, and HCC, and demonstrated that chronic hepatitis C is a relevant model of accelerated replicative senescence and that accumulation of replicative senescent cells predispose to HCC progression[18]. Hepatic stellate cell activation and senescence also trigger the development of liver cirrhosis towards HCC[19]. Yildiz et al[20] found that cirrhosis and HCC exhibit expression patterns compatible with senescent and immortal phenotypes, respectively, while dysplasia is a transitional state. Senescence bypass exerts an essential role in hepatocellular carcinogenesis engendering systematic alteration in the transcription of genes modulating DNA repair, proliferation, differentiation, metabolism, etc[20]. Eggert et al[21] reported that while chemokines secreted by senescent hepatocytes inhibit liver cancer initiation, they enable to facilitate the growth of fully established HCC[21]. Due to the highly heterogeneous malignancy at the molecular and histological levels, characterization of cellular senescence-based classification might facilitate the personalized treatment of HCCs. Recently, cell senescence molecular subtypes have been conducted for predicting prognostic outcomes and immunotherapeutic responses of hepatitis B virus-related HCC patients[22]. A cellular-senescence-related classifier has been developed for inferring predicting prognosis, immunotherapeutic responses, and candidate agents in HCCs[23]. However, these findings are based upon retrospective analysis, and lack of experimental validation. To address these problems, our integrative analysis classified HCCs as three cellular senescence subtypes and defined a cellular senescence-relevant scoring system, which decoded the tumor heterogeneity as well as tailored the pharmacological interventions to boost clinical management of HCC.

MATERIALS AND METHODS

Acquisition of cellular senescence genes

Totally, 279 human cellular senescence genes were acquired from the CellAge database (https://genomics.senescence.info/cells/)[3,4]. Genes that induce cellular senescence present the overexpression with age in human tissue samples and are notably overrepresented in antilongevity and tumor-suppressor genes; meanwhile, genes that inhibit cellular senescence overlap with prolongevity genes and oncogenes. The detailed information is listed in Supplementary Table 1.

Public HCC datasets

HCC patients were acquired from three public datasets, covering the Cancer Genome Atlas (TCGA-LIHC) database (https://portal.gdc.cancer.gov/projects/TCGA-LIHC) (n = 368), the International Cancer Genome Consortium (ICGC) portal (https://dcc.icgc.org/projects/LIRI-JP/) (n = 232), and the GSE14520 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14520). Clinicopathological traits of above datasets are summarized in Supplementary Table 2. All expression data were transformed o transcripts per kilobase million, followed by log-2 conversion.

Differential expression analysis

Utilizing limma package, differentially expressed cellular senescence genes were screened in HCC relative to normal liver tissues[24]. To prevent high false-positive rate, P values were adjusted via Benjamini–Hochberg approach. Adjusted P < 0.01 and |log₂ fold change (FC)| > 0.58 were regarded as the criteria of differentially expressed genes.

Functional annotation analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were annotated by use of clusterProfiler package[25]. GO terms and KEGG pathways with adjusted P < 0.05 were significantly enriched. The activity of fifty hallmark pathways was computed through GSVA package[26] based upon the gene sets of Molecular Signatures Database[27].

Unsupervised clustering

Based upon the prognostic differentially expressed cellular senescence genes derived from univariate Cox regression (P < 0.05), unsupervised clustering was implemented for TCGA-LIHC patients utilizing ConsensusClusterPlus package[28]. This process was conducted with 1000 iterations through sampling 80% of all the data for each iteration, thus ensuring clustering stability. The optimal number of clusters was identified utilizing consensus heatmap together with cumulative distribution function (CDF) curves. Principal component analysis (PCA) was utilized for recognizing and visualizing distinct subtypes.

Nearest template prediction for subtype verification

Nearest template prediction (NTP) method is flexible for evaluating class prediction confidence for patients. Up-regulated genes were regarded as markers of each subtype with adjusted P < 0.05, which were adopted in the NTP method derived from CMScaller package[29], thus assessing the reliability and stability of subtypes.

Tumor microenvironment estimation

Single-cell gene set enrichment analysis (ssGSEA), a deconvolution algorithm from GSVA package, was executed for quantifying the compositions within the tumor microenvironment (TME), comprising 22 immune cells and two stromal components (fibroblasts and endothelial cells). The ssGSEA score denoted the abundance of these TME components. The abundance of the TME components was also inferred through TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, together with EPIC methods.

Genetic alteration analysis

The available mutation annotation format files from the TCGA were adopted for the analysis of somatic mutation utilizing maftools package[30]. Copy number variations (CNVs) of TCGA HCCs were stratified into three cellular senescence subtypes. Significant amplifications or deletions across the whole genome were assessed utilizing GISTIC2.0[31].

Cellular senescence subtype-relevant gene selection

Genes with differential expression between subtypes were selected with the thresholds of adjusted P < 0.05 together with |log₂FC| > 0.58. Cellular senescence subtype-relevant genes were determined following the intersection.

Gene signature establishment

Univariate Cox regression analysis was utilized for picking out cellular-senescence-relevant genes with prognostic implication based upon P < 0.05. Through least absolute shrinkage and selection operator (LASSO) algorithm, the best gene subset was found out via glmnet package[32]. The cellular senescence-relevant scoring system was defined following the formula: RiskScore = Σ(coefficient (β)*Expression β), where β denoted each selected prognostic cellular senescence-relevant gene. HCCs were stratified into low- and high-RiskScore groups with the median RiskScore.

Nomogram construction

Uni- and multivariate Cox regression analyses on the cellular-senescence-relevant gene signature and conventional clinicopathological variables were executed to select independent prognostic factors in the TCGA-LIHC cohort. A nomogram based upon independent factors was generated to predict the probability of overall survival through rms package. Decision curve analysis was conducted for validating the nomogram[33].

Therapeutic response prediction

The half-maximal inhibitory concentration (IC₅₀) value of commonly used chemotherapy or targeted therapy drugs was inferred utilizing pRRophetic package[34]. Immunotherapy response was inferred by use of Tumor Immune Dysfunction and Exclusion (TIDE)[35].

Patients and tissues

Thirty fresh HCC tumors together with adjacent normal tissues were harvested from The Affiliated Bozhou Hospital of Anhui Medical University. No patients experienced any preoperative adjuvant treatment. HCC diagnosis was confirmed pathologically. Written informed consent was provided by each patient. This project gained the approval of the Ethics Committee of The Affiliated Bozhou Hospital of Anhui Medical University (2022-17).

Immunohistochemistry

TRNP1 expression in HCC or normal tissues was tested through fixing tissue sections with 4% paraformaldehyde. The sections were sealed utilizing goat serum, followed by incubation with TRNP1 antibody (1:500; ab174303; Abcam, Cambridge, MA, USA) along with secondary antibody (1/1000; ab7090). After administration with diaminobenzidine tetrahydrochloride, images were acquired under a microscope (Zeiss, Germany).

Immunoblotting

Total protein extraction was analyzed utilizing immunoblotting. Protein content was measured utilizing BCA kit (Beyotime, Shanghai, China). Proteins were separated via 12% SDS-PAGE, and transferred onto PVDF membranes that were then probed with primary antibody against TRNP1 (1/5000; ab174303; Abcam), p16 (1/500; ab151303), p21 (1/1000; ab109199) or GAPDH (1/2500; ab9485) at 4°C overnight, and secondary antibody (1/1000; ab7090) at room temperature for 2 h. Proteins were developed using ECL reagent (Beyotime).

Quantitative real-time polymerase chain reaction

RNA extraction was achieved utilizing RNA easy mini kit (Invitrogen, Carlsbad, CA, USA), with cDNA preparation via PrimeScript RT Master Mix (Takara, Dalian, China). Quantitative real-time polymerase chain reaction (RT-qPCR) was conducted via ChamQ SYBR qPCR Master Mix (Vazyme, Nanjing, China). The relative expression value was estimated with 2^−ΔΔCt approach as well as normalized to endogenous GAPDH.

Cell culture and transfection

RPMI-1640 medium (Gibco) containing 10% fetal bovine serum (Gibco), and 1% penicillin–streptomycin (Gibco) was adopted for culturing SMMC-7721 and HepG2 HCC cells. All cells were maintained in an incubator at 37°C with 5% CO₂. For transfection, siRNAs of TRNP1 (si-TRNP1) and negative control (si-NC) were acquired from GenePharma. Cell transfection was conducted utilizing Lipofectamine 2000 (Thermo Fisher Scientific).

Flow cytometry

Apoptotic rate was tested through flow cytometry utilizing Annexin V–fluorescein isothiocyanate (apoptosis detection kit (BD Biosciences). Cells were harvested and the assay was performed. Next, samples were assessed instantly utilizing flow cytometry (Beckman Coulter).

SA-β-galactosidase staining

To investigate senescence, 10⁴ cells were seeded onto a six-well plate. After being fixed, they were stained with senescence-associated -galactosidase activity (SA-β-gal) (Gibco).

Tumor xenograft

Female BALB/c nude mice (5-wk-old, 16–18 g; Beijing Vital River Laboratory Animal Technology Co. Ltd., China) were fed under a 12-h light/dark cycle. They were divided into three groups (n = 5 per group). SMMC-7721 cells (n = 10⁵) with si-NC, si-TRNP1#1 or si-TRNP1#2 were inoculated into the armpit. After 36 d, they were killed, with subsequent tumor excision. Tumor volume was finally calculated. This experiment gained the approval the Animal Ethics Committee of The Affiliated Bozhou Hospital of Anhui Medical University (LLSC20232071).

Statistical analysis

For between-group comparisons, unpaired Student’s t-test was adopted; Mann–Whitney U test was utilized for variables with non-normal distribution. Kaplan–Meier curves were utilized to estimate the overall survival of groups, with log-rank for testing the difference significance between groups. Survival analysis was executed utilizing survival and survminer packages. Receiver operator characteristic curves were plotted to evaluate the prediction efficacy in overall survival. Statistical analyses were implemented utilizing R packages and GraphPad Prism software. A two-sided P < 0.05 indicated statistical significance.

RESULTS

Differentially expressed cellular senescence genes in HCC and biological significance

We identified 146 differentially expressed cellular senescence genes in TCGA HCCs relative to normal tissues with adjusted P < 0.01 and |log₂FC| > 0.58 (Figure 1A and B; Supplementary Table 3), which might participate in HCC initiation or progression. They were linked with metabolic process, cellular senescence, cell cycle, and immunity pathways (Figure 1C–F). Their prognostic value was then assessed. Ninety-seven differentially expressed cellular senescence genes were significantly linked with HCC prognosis (Table 1).

Figure 1 Differentially expressed cellular senescence genes in hepatocellular carcinoma and biological significance. A: Volcano diagram of the abnormal expression of cellular senescence genes in hepatocellular carcinoma (HCC) relative to normal liver tissues in the TCGA-LIHC dataset. Blue dots denote differentially expressed cellular senescence genes, with grey dots for the not differentially expressed genes; B: Heatmap of the transcript levels of differentially expressed cellular senescence genes in HCC and normal liver tissues; C–E: Bubble diagrams of the top 10 biological process, cellular component, molecular function terms enriched by differentially expressed genes. The bubble size indicates the count of genes enriched. The closer the color is to red, the smaller the adjusted p; F: Circle graph of the Kyoto Encyclopedia of Genes and Genomes pathways enriched by differentially expressed genes.

Table 1 Univariate-cox regression results of prognostic differentially expressed cellular senescence genes with P < 0.05 in TCGA-LIHC dataset.

Gene	HR	95%lower	95%upper	P value	Gene	HR	95%lower	95%upper	P value
MCRS1	1.6954	1.2504	2.2988	0.0007	TPR	1.2593	1.0274	1.5436	0.0264
FASTK	1.3714	1.0147	1.8536	0.0399	HDAC1	1.8639	1.4442	2.4055	< 0.0001
AURKA	1.2794	1.1083	1.4770	0.0008	SPOP	1.4186	1.0061	2.0004	0.0461
PTTG1	1.3361	1.1782	1.5151	< 0.0001	BTG3	1.3304	1.0766	1.6440	0.0082
PSMB5	1.8929	1.3286	2.6969	0.0004	IRF5	1.4232	1.1118	1.8217	0.0051
MAP2K2	1.5424	1.1846	2.0083	0.0013	IGFBP3	1.1741	1.0435	1.3210	0.0076
E2F1	1.2175	1.0804	1.3721	0.0012	CDKN2B	1.3951	1.1619	1.6751	0.0004
CDK1	1.3097	1.1553	1.4847	< 0.0001	AAK1	1.3018	1.0160	1.6682	0.0371
GRK6	1.7940	1.3455	2.3920	< 0.0001	ARPC1B	1.3007	1.0918	1.5495	0.0032
EZH2	1.5722	1.3117	1.8845	< 0.0001	LIMK1	1.3450	1.1517	1.5706	0.0002
DPY30	1.9266	1.3785	2.6928	0.0001	SFN	1.1541	1.0714	1.2432	0.0002
CBX8	1.4643	1.1477	1.8684	0.0022	CSNK2A1	1.3825	1.0896	1.7541	0.0077
SMARCA4	1.4396	1.1405	1.8170	0.0022	GLB1	1.4265	1.1097	1.8336	0.0056
CDKN2A	1.1748	1.0571	1.3055	0.0028	MOB3A	1.4527	1.1459	1.8416	0.0020
IRF3	1.4048	1.0834	1.8216	0.0103	BRD7	1.3897	1.0401	1.8567	0.0260
HRAS	1.4678	1.1958	1.8016	0.0002	PRKCD	1.4989	1.2458	1.8035	< 0.0001
ADCK5	1.3095	1.0375	1.6527	0.0232	CDK4	1.4798	1.2301	1.7800	< 0.0001
RUVBL2	1.7508	1.3029	2.3528	0.0002	BLVRA	1.1790	1.0296	1.3500	0.0172
ACLY	1.3734	1.1160	1.6901	0.0027	SENP1	1.5189	1.1628	1.9840	0.0022
TACC3	1.3618	1.1752	1.5780	< 0.0001	BMI1	1.6112	1.2516	2.0742	0.0002
SIRT6	1.6515	1.2586	2.1671	0.0003	DHX9	1.3908	1.1069	1.7477	0.0046
SUPT5H	1.5227	1.1449	2.0252	0.0039	RSL1D1	1.7134	1.2421	2.3634	0.0010
FOXM1	1.2772	1.1287	1.4454	0.0001	PAK4	1.2908	1.0419	1.5990	0.0195
PSMD14	1.9715	1.5043	2.5838	< 0.0001	PDCD10	1.6080	1.2369	2.0903	0.0004
HJURP	1.4529	1.2514	1.6869	< 0.0001	ASF1A	1.7774	1.3749	2.2977	< 0.0001
TRIM28	1.6005	1.2816	1.9988	< 0.0001	PNPT1	1.6790	1.2487	2.2576	0.0006
P3H1	1.9307	1.5257	2.4431	< 0.0001	MAPK12	1.2163	1.0478	1.4119	0.0101
MAGOHB	1.4601	1.0477	2.0346	0.0254	UBTD1	1.3794	1.0584	1.7977	0.0173
RBX1	1.6692	1.2410	2.2451	0.0007	KDM5B	1.3241	1.1029	1.5895	0.0026
MAGOH	1.6449	1.1977	2.2592	0.0021	USP1	1.4478	1.1790	1.7778	0.0004
CENPA	1.4947	1.2958	1.7241	< 0.0001	MAP3K7	1.4740	1.1327	1.9180	0.0039
EWSR1	1.9373	1.3937	2.693	< 0.0001	PKM	1.2079	1.0979	1.3290	0.0001
HSPA5	1.3234	1.0428	1.6796	0.0212	NINJ1	1.4743	1.0929	1.9888	0.0110
CHEK1	1.5484	1.2766	1.8780	< 0.0001	SERPINE1	1.1201	1.0226	1.2269	0.0147
PIAS4	1.4493	1.0683	1.9663	0.0171	BRCA1	1.4378	1.1653	1.7741	0.0007
PRPF19	2.5587	1.7707	3.6974	< 0.0001	DEK	1.2103	1.0125	1.4468	0.0360
MAPKAPK5	2.0718	1.4388	2.9834	< 0.0001	NDRG1	1.2946	1.1412	1.4686	< 0.0001
MAD2L1	1.4390	1.2065	1.7164	< 0.0001	SRC	1.2011	1.0507	1.3730	0.0073
TFAP4	1.9651	1.4205	2.7184	< 0.0001	ASPH	1.2007	1.0488	1.3746	0.0081
G6PD	1.3897	1.2500	1.5450	< 0.0001	STK32C	1.5149	1.2235	1.8759	0.0001
GAPDH	1.5984	1.2892	1.9818	< 0.0001	CDK2AP1	1.2327	1.0010	1.5181	0.0490
SMARCB1	1.4393	1.1297	1.8336	0.0032	KDM4A	1.3363	1.0476	1.7045	0.0196
LEO1	1.4005	1.0359	1.8934	0.0286	MMP9	1.1360	1.0384	1.2428	0.0054
TXN	1.3040	1.0667	1.5941	0.0096	HK3	1.2356	1.0293	1.4832	0.0232
HDAC4	1.5029	1.1470	1.9692	0.0031	VEGFA	1.2906	1.0657	1.5630	0.0090
FXR1	1.6935	1.2509	2.2927	0.0007	LGALS3	1.2041	1.0826	1.3392	0.0006
RAD21	1.3480	1.0928	1.6629	0.0053	AR	0.8662	0.7749	0.9682	0.0115
SRSF1	1.7610	1.2427	2.4955	0.0015	BAG3	1.2740	1.0326	1.5718	0.0238

Classification of HCC patients as three cellular senescence subtypes

Prognostic differentially expressed cellular senescence genes were utilized for probing HCC heterogeneity. Utilizing unsupervised clustering, TCGA-LIHC cases were initially assigned to 2–9 clusters. Combining consensus CDF and consensus matrix, the optimal number of clusters was generated when k = 3 (Figure 2A–C). Thus, HCCs were classified as three cellular senescence subtypes, named C1–3. Prognostic differentially expressed cellular senescence genes presented the highest transcript levels in C2, followed by C1 and C3 (Figure 2D). PCA proved the extensive discrepancy in transcript levels among three subtypes (Figure 2E). Additionally, we focused on the survival difference, with C2 having the worst overall survival, C1 the next, and C3 the best (Figure 2F). Based upon up-regulated markers of each subtype (Supplementary Table 4), the robustness and reproducibility of cellular senescence subtypes were verified utilizing NTP in the ICGC cohort (Figure 2G). The discrepancy in transcript levels and overall survival among subtypes was further proven in this cohort (Figure 2H and I).

Figure 2 Classification of TCGA-LIHC patients as three cellular senescence subtypes and external dataset validation. A–C: Consensus cumulative distribution function (CDF), relative alteration in area under CDF curve, and consensus matrix k = 3 based upon the transcriptome of prognostic differentially expressed cellular senescence genes across TCGA-LIHC patients; D: The transcript levels of differentially expressed cellular senescence genes with prognostic implications across three subtypes; E: Principal component analysis (PCA) plots of the discrepancy in transcript levels among subtypes; F: Kaplan–Meier (K-M) curves of overall survival in TCGA-LIHC; G: Nearest template prediction for verifying the subtypes in the International Cancer Genome Consortium (ICGC) cohort; H and I: PCA plots of the transcriptome difference and K-M curves of overall survival among subtypes in the ICGC cohort.

Responses to immunotherapy and chemotherapy of three cellular senescence subtypes

To elucidate the underlying mechanisms among the three cellular senescence subtypes, the activity of 50 hallmark pathways was inferred. Tumorigenic pathways (DNA repair, MYC, PI3K–AKT–mTOR, mTORC1, etc.) exhibited the highest activity in C2, with the lowest activity of metabolism pathways (Figure 3A). C3 presented the lowest activity of tumorigenic pathways, as well as the highest activity of metabolism pathways. Additionally, it was found that immune checkpoints displayed the highest transcript levels in C2, with the highest abundance of immune cells (Figure 3B). Immunogenetic indicators were then observed. Aneuploidy score, cancer-testicular antigen score, homologous recombination defects, and intratumor heterogeneity displayed the highest levels in C2, followed by C1 and C3 (Figure 3C–F). TIDE score was utilized to estimate the response to immune checkpoint inhibitors. Among three subtypes, C3 presented the lowest TIDE score, indicating that this subtype was most likely to respond to immune checkpoint inhibitors (Figure 3G). It was also found that cisplatin, doxorubicin and gemcitabine showed the lowest IC₅₀ values in C2 subtype (Figure 3H–J). Thus, C2 patients were most likely to benefit from above chemotherapeutic drugs.

Figure 3 Immunogenomic landscape of three cellular senescence subtypes across TCGA-LIHC. A: The activity of 50 hallmark pathways in three subtypes; B: Abundance of the tumor microenvironment components, transcript levels of immune checkpoints, and stromal and immune scores across subtypes; C–F: Differences in aneuploidy score, cancer-testicular antigen score, homologous recombination defects, and intratumor heterogeneity between subtypes; G: Difference in TIDE score between subtypes; H–J: Differences in IC₅₀ values of cisplatin, doxorubicin and gemcitabine between subtypes. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001. NS: Not significant; TIDE: Tumor Immune Dysfunction and Exclusion.

Genetic alterations of three cellular senescence subtypes

Overall, somatic mutation rate was the lowest in C2 among three cellular senescence subtypes (Figure 4A–C). Additionally, C1 presented the highest copy number amplified and deleted alterations (Figure 4D–I). Altogether, there was remarkable heterogeneity in genetic alterations among three cellular senescence subtypes.

Figure 4 Genetic alterations of three cellular senescence subtypes. A–C: Significant mutated genes in The Cancer Genome Atlas hepatocellular carcinomas stratified by cellular senescence subtypes; D–F: Copy number amplification plots in three subtypes. The green line denotes the significance threshold (q = 0.25); G–I: Copy number deletion plots in three subtypes.

Identification of cellular senescence subtype-relevant genes and functional implications

To select cellular senescence subtype-relevant genes, we assessed the genes with differential expression between cellular senescence subtypes based upon adjusted P < 0.05 together with |log₂FC| > 0.58. After the intersection, 666 cellular senescence subtype-relevant genes were eventually acquired (Supplementary Table 5 and Figure 5A). We elucidated the underlying functional implications. Consequently, these cellular senescence subtype-relevant genes were remarkably linked with cell cycle, DNA replication, oocyte meiosis, homologous recombination, cellular senescence, Fanconi anemia pathway, p53 pathway, progesterone-mediated oocyte maturation, mismatch repair, etc (Figure 5B–E).

Figure 5 Identification of cellular senescence subtype-relevant genes and functional implications. A: Venn diagram of the intersected genes from differentially expressed genes between cellular senescence subtypes in TCGA-LIHC cohort; B–D: The top 10 biological processes, cellular components, together with molecular functions of cellular senescence subtype-relevant genes; E: Kyoto Encyclopedia of Genes and Genomes pathways significantly enriched by these genes.

Definition and external verification of a cellular-senescence-relevant gene signature

To illustrate the relationships of the cellular-senescence-relevant genes and patient survival, univariate Cox regression method was adopted. A total of 511 cellular -enescence-relevant genes presented significant correlations to TCGA-LIHC prognosis (Supplementary Table 6). These prognostic genes were entered into LASSO analysis (Figure 6A and B). A 19-gene signature was generated in accordance with the optimal λ value. The cellular senescence-relevant scoring system was computed as follows: RiskScore = 0.0610156 * transcript level of SLC1A5 + 0.049731458 * transcript level of G6PD + 0.038762092 * transcript level of PSRC1 + 0.104396819 * transcript level of UCK2 + 0.004054037 * transcript level of TCOF1 + 0.03040522 * transcript level of CCT5 + 0.002669582 * transcript level of DTYMK + 0.053080689 * transcript level of NEIL3 + 0.031271078 * transcript level of TRNP1 + (-0.037524429) * transcript level of ADH4 + 0.018409162 * transcript level of HMMR + 0.001793118 * transcript level of SMG5 + (-0.016621762) * transcript level of CLEC3B + 0.083327624 * transcript level of PLOD2 + 0.035535231 * transcript level of SPP1 + (-0.023169364) * transcript level of CFHR3 + 0.020042596 * transcript level of TMEM106C + (-0.019959706) * transcript level of ANXA10 + (-0.039577478) * transcript level of LCAT. Based upon the median RiskScore, TCGA-LIHC cases were classified as low- and high-RiskScore groups (Figure 6C). Expression levels of these selected genes exhibited the notable differences between groups. Next, K-M curves illustrated that high-RiskScore patients’ overall survival was worse (Figure 6D). AUCs at 1-, 3- and 5-year survival all exceeded 0.75, demonstrating the excellent discrimination power of the gene signature (Figure 6E).

Figure 6 Definition and external verification of a cellular-senescence-relevant gene signature. A: Cross-validation for tuning the parameter selection in the least absolute shrinkage and selection operator (LASSO) analysis in TCGA-LIHC dataset; B: LASSO coefficient profiling; C: Distribution of RiskScore, survival status, and transcript levels of genes in the signature; D: Kaplan–Meier (curves of overall survival in low- and high-RiskScore hepatocellular carcinomas; E: Receiver operator characteristic curves (ROCs) at 1-, 3- and 5-year survival; F: External verification of distribution of RiskScore, survival status, and transcript levels of genes in the International Cancer Genome Consortium (ICGC) dataset; G: Kaplan–Meier curves of overall survival in two groups in the ICGC dataset; H: ROCs at 1- and 3-year survival based upon RiskScore in the ICGC dataset; I: External validation of distribution of RiskScore, survival status, and transcript levels of genes in the GSE14520 dataset; J: Kaplan–Meier curves of overall survival in two groups in the GSE14520 dataset; K: ROCs at 1- and 3-year survival in the GSE14520 dataset.

The ICGC and GSE14520 datasets were utilized to externally verify this signature. The current study stratified HCCs into low- and high-RiskScore groups based upon the median RiskScore in the ICGC dataset (Figure 6F). Overall survival rate of high-RiskScore group was prominently lower (Figure 6G). AUCs at 1- and 3-year survival were > 0.75 (Figure 6H). Above data proved the high reproducibility of the signature. The similar findings were also confirmed in the GSE14520 dataset (Figure 6I–K).

Generation of a prognostic nomogram for clinical practice

Univariate and multivariate Cox regression analyses were executed to select the independent prognostic parameters for HCCs. It was found that stage and the cellular-senescence-relevant RiskScore acted as independent risk factors of HCC prognosis (Figure 7A and B). As a visual representation of the prognostic model, a nomogram containing stage, and RiskScore was built to illustrate HCC patients’ survival more intuitively. The nomogram showed that RiskScore had the highest influence on 1-, 3- and 5-year survival of HCC patients, followed by stage (Figure 7C). Decision curve analysis demonstrated that the nomogram can accurately predict 1-, 3- and 5-year clinical outcomes (Figure 7D–F).

Figure 7 Establishment of a prognostic nomogram for clinical practice in the TCGA-LIHC dataset. A: Univariate-cox regression results of the cellular-senescence-relevant gene signature and conventional clinicopathological parameters with hepatocellular carcinoma prognosis; B: Multivariate Cox regression for selecting independent prognostic parameters; C: Generation of a nomogram based on stage and the cellular senescence-relevant RiskScore; D–F: Decision curve analysis at 1-, 1-, and 5-year survival.

Assessment of the cellular-senescence-relevant gene signature in predicting efficacy of pharmacological interventions

The IC₅₀ of some chemotherapy or targeted therapy agents was estimated in TCGA-LIHC dataset. High-RiskScore HCCs presented the notably lower IC₅₀ of cisplatin, doxorubicin, and gemcitabine relative to those with low RiskScore (Figure 8A–C). However, no significant difference in the IC₅₀ of sorafenib was found between the two groups (Figure 8D). Accordingly, high-RiskScore HCCs more possibly responded to cisplatin, doxorubicin or gemcitabine chemotherapy.

Figure 8 Assessment of the cellular senescence-relevant gene signature in predicting efficacy of pharmacological interventions in TCGA-LIHC dataset. A–D: IC₅₀ value of cisplatin, doxorubicin, gemcitabine, and sorafenib in low- and high-RiskScore hepatocellular carcinomas (HCCs); E: Abundance of the tumor microenvironment components inferred by multiple algorithms; F–K: Comparison of carcinoma-associated fibroblast (CAF), myeloid-derived suppressor cell (MDSC), interferon gamma (IFNG), dysfunction score, exclusion score and Tumor Immune Dysfunction and Exclusion levels in low- and high-RiskScore HCCs. ^cP < 0.001. ns: No significant difference.

Some reliable computational approaches were adopted to infer the abundance of the TME elements across TCGA-LIHC samples. Overall, most immune cells exhibited the higher infiltration in high-risk HCCs (Figure 8E). The TIDE method was used to predict immunotherapy response. We did not observe any difference in carcinoma-associated fibroblasts between the groups (Figure 8F). Lower myeloid-derived suppressor cells, interferon gamma, exclusion score and TIDE score as well as higher dysfunction score were found in low-RiskScore HCCs (Figure 8G–K). This indicated that low-RiskScore HCCs more possibly benefited from immunotherapy.

Experimental verification of expression of TRNP1 in HCC

Among the genes in the cellular-senescence-relevant gene signature, the role of TRNP1 in HCC remains unclear. Therefore, we focused on TRNP1. It was proven that TRNP1 presented remarkable upregulation in HCCs relative to normal tissues in accordance with immunohistochemistry (Figure 9A and B), immunoblotting (Figure 9C and D) and RT-qPCR (Figure E). Specific siRNAs of TRNP1 were transfected into SMMC-7721 and HepG2 cells. Immunoblotting demonstrated the notable decrease in TRNP1 expression induced by siRNAs (Figure 9F–9H).

Figure 9 Experimental verification of expression of TRNP1 in hepatocellular carcinomas. A and B: Representative immunohistochemistry of TRNP1 expression in human hepatocellular carcinomas (HCCs) and normal tissues. Bar, 50 μm; C and D: Representative immunoblotting of TRNP1 expression in human HCCs and normal tissues; E: Quantitative real-time polymerase chain reaction of TRNP1 expression in 30 pairs of human HCCs and normal tissues; F–H: Immunoblotting of TRNP1 expression in SMMC-7721 and HepG2 cells transfected with specific siRNAs of TRNP1. ^cP < 0.001; ^dP < 0.0001.

Suppression of TRNP1 induces apoptosis and senescence of HCC cells and attenuates tumor growth

Based upon flow cytometry, apoptotic rate of SMMC-7721 and HepG2 cells was prominently elevated by TRNP1 knockdown (Figure 10A–C). SA-β-Gal staining showed that TRNP1 knockdown notably induced cellular senescence of two HCC cells (Figure 10D–F). Additionally, two cellular senescence markers: p16 and p21 were overexpressed in HCC cells with TRNP1 knockdown (Figure 10G–K), further proving the role of TRNP1 in HCC senescence. In vivo tumorigenicity models were also developed for evaluating whether TRNP1 influenced tumor growth. As a result, TRNP1 knockdown was found to decrease in vivo tumor volume (Figure 10L and M).

Figure 10  Suppression of TRNP1 induces apoptosis and senescence of hepatocellular carcinoma cells and attenuates tumor growth. A–C: Flow cytometry for measuring the apoptotic rate of in SMMC-7721 and HepG2 cells with transfection of specific siRNAs of TRNP1; D-F: SA-β-galactosidase (SA-β-Gal) staining for evaluating senescence of transfected hepatocellular carcinoma cells. Bar, 10 μm; G–K: Immunoblotting of p16 and p21 expression in transfected cells; L: Representative photographs of tumors from mice of si-NC, si-TRNP1#1 and si-TRNP1#2 groups; M: Calculation of tumor volume in above groups. ^cP < 0.001; ^dP < 0.0001.

DISCUSSION

Cellular senescence is a permanent state of cell cycle arrest occurring in proliferating cells when face distinct stresses[36]. In cancers, senescence is usually an effective barrier against tumorigenesis because it prevents the division potential of cells[37]. Nonetheless, numerous research has demonstrated that senescent cells also have tumorigenic properties[38]. Thus, it is of significance to characterize key features of cellular senescence in HCC.

HCC is a typically fatal malignant tumor displaying genetic heterogeneity and limited therapy responses[39]. Based upon prognostic differentially expressed cellular senescence genes, we classified HCCs as three cellular senescence subtypes: C1–C3. The robustness and reproducibility of this classification were externally proven. C2 had the worst overall survival, C1 the next, and C3 the best, revealing the heterogeneity in prognostic outcomes among subtypes. Single-agent anti-PD-1 immune checkpoint blockade showed ponent efficacy in early-phase trials, but the findings were not confirmed in phase III studies[40]. In accordance with the lowest TIDE score, and immunogenetic indicators, C3 HCCs might possibly respond to immunotherapy. Additionally, C2 HCCs were most likely to benefit from chemotherapy. Thus, this classification might assist clinical decision-making. Genetic mutations associate with HCC initiation and progression[41-43]. For instance, mutant TP53 is the most frequent in HCC, affecting patient survival, and immune response[44]. CTNNB1 mutation occupies a large proportion of human HCCs, which correlates to high TMB and AFP in HCCs[45]. Among three cellular senescence subtypes, C2 presented the lowest somatic mutation rate, while C1 had the highest frequent CNVs. Accordingly, cellular senescence subtypes appear to associate with genetic mutations.

We defined a novel cellular-senescence-relevant gene signature comprising SLC1A5, G6PD, PSRC1, UCK2, TCOF1, CCT5, DTYMK, NEIL3, TRNP1, ADH4, HMMR, SMG5, CLEC3B, PLOD2, SPP1, CFHR3, TMEM106C, ANXA10, and LCAT, with the excellent power in survival prediction in HCCs. Previous research has proven the biological implications of the cellular-senescence-relevant genes in HCC. For example, SLC1A5 regulated by DDR1 contributes to HCC progression[46]. G6PD weakens ferroptosis in HCC through targeting cytochrome P450 oxidoreductase[47]. PSRC1, a hypoxia- and immune-associated gene, associates with HCC survival[48]. The nonmetabolic role of UCK2 facilitates HCC metastasis via EGFR–AKT signaling activation[49]. TCOF1 coordinates oncogenic activation and rRNA generation as well as results in HCC initiation[50]. Clinical features of HCC patients are notably associated with survival outcomes. To better optimize the cellular-senescence-relevant gene signature and improve the prediction accuracy, we incorporated stage in combination with the cellular-senescence-relevant gene signature to build a nomogram that enabled us to generate the individual survival probability in HCC patients. High-RiskScore HCCs might respond to cisplatin, doxorubicin or gemcitabine, while low-RiskScore HCCs more possibly benefit from immunotherapy, proving the potential of the cellular-senescence-relevant gene signature in inferring therapeutic efficacy. Among the cellular-senescence-relevant genes, the expression and role of TRNP1 in HCC remain indistinct. Only bioinformatics evidence demonstrated the prognostic significance of TRNP1 in HCC[51]. Our study experimentally proved that TRNP1 was upregulated in HCC, and TRNP1 knockdown induced apoptosis and senescence of HCC cells and attenuated tumor growth. Thus, TRNP1 potentially participates in HCC senescence and progression, which might be a promising therapeutic target.

Our study had some limitations. Due to the lack of HCC patients who received neoadjuvant immunotherapy, the relationship between cellular senescence subtypes and relevant gene signature with immunotherapeutic response requires further verification in the immunotherapy cohorts. Despite the external verification in the ICGC dataset, the predictive efficacy of cellular senescence-relevant gene signature needs to be proven in prospective cohorts.

CONCLUSION

Our findings showed the importance of cellular senescence in HCC classification and pharmacological interventions for clinical management. Additionally, we defined a cellular-senescence-relevant scoring system that can infer patient survival and therapeutic efficacy. Considering the clinically relevant parameters were closely linked with HCC survival, we incorporated stage in combination with the cellular-senescence-relevant gene signature to build a nomogram. Our integrated analysis provides a valuable framework for comprehending cellular senescence in HCC, which sheds light on the senescence-associated biomarker discovery as well as therapeutic targets.

ARTICLE HIGHLIGHTS

Research background

Cellular senescence, a state of stable growth arrest, is intertwined with human cancers. Due to the highly heterogeneous malignancy at the molecular and histological levels, characterization of cellular-senescence-based classification might facilitate the personalized treatment of hepatocellular carcinoma (HCC).

Research motivation

Nonetheless, the heterogeneity of cellular-senescence-related features makes the definition and targeting of treatment-induced senescent cells challenging.

Research objectives

This study aimed to characterize cellular-senescence-based phenotypes in HCC, and identify a novel cellular-senescence-related therapeutic target.

Research methods

We enrolled two HCC datasets, TCGA-LIHC and International Cancer Genome Consortium (ICGC). Unsupervised clustering was executed to probe tumor heterogeneity based upon cellular senescence genes. Least absolute shrinkage and selection operator algorithm was utilized to define a cellular-senescence-relevant scoring system. TRNP1 expression was measured in HCCs and normal tissues through immunohistochemistry, immunoblotting and quantitative real-time polymerase chain reaction. The influence of TRNP1 on HCC senescence and growth was proven via a series of experiments.

Research results

TCGA-LIHC patients were classified as three cellular senescence subtypes, named C1–3. The robustness and reproducibility of these subtypes were proven in the ICGC cohort. C2 had the worst overall survival, C1 the next, and C3 the best. C2 presented the highest levels of immune checkpoints, abundance of immune cells, and immunogenetic indicators. Thus, C2 might respond to immunotherapy. C2 had the lowest somatic mutation rate, while C1 presented the highest copy number variations. A cellular-senescence-relevant gene signature was generated, which can predict patient survival, and chemo- or immunotherapeutic response. Experimentally, it was proven that TRNP1 presented with remarkable upregulation in HCCs. TRNP1 knockdown induced apoptosis and senescence of HCC cells and attenuated tumor growth.

Research conclusions

These findings provide a systematic framework for assessing cellular senescence in HCC, which decode the tumor heterogeneity and tailor the pharmacological interventions to improve clinical management.

Research perspectives

Cellular senescence, a state of stable growth arrest, is implicated in human cancers. Nevertheless, characterization of cellular-senescence-associated phenotypes in HCC is still indistinct. Here, we proposed a novel cellular-senescence-based classification for HCC and identified TRNP1 as a novel therapeutic target.