Comprehensive analysis to develop a stromal senescence-associated gene signature for predicting hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC) is a common liver malignant tumor, which ranks sixth in terms of global incidence and second as a cause of global mortality (1). With the evolving knowledge of tumorigenesis, HCC is known to be dominantly induced by a series of chronic liver diseases, such as liver cirrhosis, hepatitis B virus/hepatitis C virus (HBV/HCV) infection, and fatty liver disease (2). Currently, more curative treatment approaches than ever have been proposed for HCC patients, including transcatheter arterial chemoembolization (TACE), radical surgery, immunotherapy, and molecular targeted agents (3). Amid the emergence of diverse treatment modalities, HCC patients now benefit from an expanded array of therapeutic options to enhance clinical outcomes. However, the disease’s aggressive biology and high recurrence rate persistently undermine patient survival, yielding a mere 18% 5-year survival rate (4). Furthermore, previous studies have demonstrated that patients at the same clinical stage exhibit significant variability in therapeutic effects and pronounced differences in prognosis (5). Multi-gene signatures have been reported to effectively improve survival outcomes and therapeutic responses to immunotherapy or chemotherapy in HCC patients (6,7). Thus, it is necessary to explore optimal predictive biomarkers to ascertain prognosis as well as therapeutic sensitivity, enabling personalized therapeutic selection and prolonged survival for HCC patients.

Cellular senescence signifies an irreversible halt in the cell cycle, specifically in diploid cells and curtailing their proliferative potential. This process functions as a defensive mechanism against genomic instability, thereby attenuating the amassment of DNA damage (8,9). Oncogene-driven stress resulting from dysregulated oncogene expression or tumor suppressor loss can also provoke senescence in cancer cells. Notably, senescence acts as a tumor-suppressive biology to impede cancer development, underscored by its application for its anti-tumor efficacy (10). Increasing evidence indicates that senescent cells can paradoxically promote tumorigenesis by stimulating cellular invasion, proliferation, and metastasis, facilitating vascularization, and disrupting antitumor immunity (11-13). Notably, tumor cells can adopt temporary and reversible senescent-like states, which mostly lead to treatment resistance (14,15). Given that cellular has been recently established as a defining feature of cancers (16), the detection and characterization of senescent cells emerge as critical strategies to prevent tumor-promoting consequences and hold significant potential for tailoring personalized cancer therapies.

The composition of stromal cells demonstrates marked variation across tumor types, encompassing endothelial cells (ECs), fibroblasts, mesothelial cells, mural cells (MCs), and glial cells. As fundamental physiological components of the tumor microenvironment (TME), stromal cells mediate tumor progression and therapeutic resistance (17). Accumulating research indicates that tumor-related stromal cells are pivotal participants within the TME, facilitating tumor growth through the secretion of growth factors, enhancement of angiogenesis, and suppression of immune responses. Building on these findings, targeted therapeutic strategies are being developed to modulate stromal cells as a means to disrupt tumor-supportive mechanisms (18).

Given the central role of stromal cells in tumorigenesis, we hypothesize that senescent stromal cells drive tumor progression and modulate therapeutic efficacy in HCC. Using machine-learning algorithms, we developed a stromal-senescence prognostic signature and externally validated it in Gene Expression Omnibus (GEO) cohorts. We systematically characterized high- and low-risk patients across multiple dimensions and profiled signature-gene expression by multi-omics analyses. Virtual screening identified three high-affinity DNASE1L3 ligands: STL167676, STK598132, and STK119225. This study provides a comprehensive framework for targeting stromal cell senescence to predict survival and optimize therapeutic responsiveness. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2894/rc).

Methods

Data source

In this investigation, eight independent publicly available datasets were acquired and processed. The single-cell RNA sequencing (scRNA-seq) cohorts, such as GSE149614 (n=13), GSE189903 (n=16), GSE156625 (n=43), and GSE151530 (n=32), were obtained from the GEO database, which collectively profiled transcriptomic landscapes of HCC tumor samples. After excluding samples with missing or unmeasurable data, we obtained The Cancer Genome Atlas (TCGA)-Liver Hepatocellular Carcinoma (LIHC) (n=369) and GSE14520 (n=221) datasets, along with RNA transcriptome profiles and comprehensive clinical annotations acquired from public databases. For the TCGA-LIHC dataset, somatic mutation profiles and copy number variation (CNV) data were extracted from the cBioPortal online platform. Protein expression matrices and relative clinical data of 159 paired HCC patients were acquired from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) cohort. Spatial transcriptome data were obtained from a previously published study (http://lifeome.net/supp/livercancer-st/data.htm), which includes 10 spatial transcriptomic samples. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Single-cell transcriptome analysis

Seurat package (version 5.2.1) was utilized for further computational analysis. Cells expressing between 500 and 6,000 unique features and containing less than 15% mitochondrial RNA content were retained for further processing. To mitigate batch effects and prevent individual patient-specific biases from dominating the analysis, the Harmony integration algorithm was implemented to harmonize the single-cell datasets. Dimensionality reduction was employed via screening out the top 2,000 highly variable genes via the FindVariableFeatures function, followed by data scaling. Principal component analysis (PCA) was employed by these variable genes. Nearest neighbors for graph clustering were achieved through the FindNeighbors and FindCluster functions based on the top 15 principal components. Cell type visualization was generated by the Uniform Manifold Approximation and Projection (UMAP) tool.

Cell annotation and differential expression analysis

Within the obtained clusters, scores were given to gene signatures specific to different cell subtypes. These signatures contained markers for T cells (CD3D, CD3E, KLRD1, and KLRC1), myeloid cells (CSF3R, C1QB, CD1C, and CLEC9A), B cells (CD79A and CD79B), stromal cells (COL1A1, COL1A2, PECAM1, and VWF), and epithelial cells (EPCAM, KRT8, and KRT18). Differential expression analysis was conducted employing the FindAllMarkers function across clusters, with parameters specified as logfc.threshold =0.5, min.pct =0.25, and only.pos = TRUE. The MAST analysis was applied to compute statistical significance, yielding adjusted P values for the pairwise comparisons.

Differentially expressed senescence genes in cells

A total of 866 senescence-related genes were retrieved from the CellAge database (https://genomics.senescence.info/cells). Detailed annotations for these 866 genes are summarized in Table S1. Overlapping gene sets between senescence genes and differentially expressed genes (DEGs) were screened out across cell types via the Venn visualization tool (19). The heatmap was expression comparison between HCC and normal tissues were generated to visualize the expressions of overlapping genes across distinct cell populations through the DoHeatmap function. Subsequently, the expression comparison between HCC and normal tissues was conducted within the TCGA cohort. Statistical significance was determined via two-tailed Student’s t-test, with P<0.05 indicating significant differential expression. We employed the AUCell_calcAUC (version 1.28.0) function to quantify cellular senescence scores for individual cells and assessed the mean senescence activity across diverse cell types. Subsequently, cells were selected into high- and low-senescence groups on the basis of the median senescence score. Differential expression analysis was employed utilizing the FindAllMarkers function, integrated with the MAST framework. Genes exclusively enriched in the high-senescence group were considered senescence-associated DEGs for downstream investigations.

Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) pathway analysis

ClusterProfiler R package (version 4.12.6) was employed to conduct KEGG pathway enrichment analysis of senescence-associated DEGs, aiming to identify significantly enriched biological processes (20). GO enrichment analysis explored the functional annotations of senescence-associated DEGs across multiple biological levels. Fisher’s exact test was utilized to ascertain the significance of each term. The threshold of significance was P<0.05.

Constructing prognostic features integrated with machine learning methods

A prognostic signature was developed utilizing the transcripts per million (TPM) profile of HCC samples from the TCGA-LIHC dataset with comprehensive survival records. Univariate Cox analysis was conducted to preliminarily obtain the OS of senescence-related genes in stromal cells from HCC patients. Subsequently, the least absolute shrinkage and selection operator (LASSO) analysis was used to screen out the prominent genes, and multivariate Cox analysis was further used to construct a prognostic signature.

Constructing the nomogram

The combination of risk score with pathological characteristics was performed via the RMS package (version 7.0.0) to develop a nomogram and corresponding calibration curves. It was assessed by aligning the predictive survival probabilities with the observed survival outcomes, where the ideal 45-degree diagonal line signifies the optimal predictive accuracy. Moreover, the survival package was applied to generate a forest plot, visually depicting the prognostic impacts of individual clinical variables, containing tumor-node-metastasis (TNM) and China Liver Cancer (CNLC) stages, and risk score on patient outcomes.

Immune checkpoint analysis

The expression of common immune checkpoints (ICG_genes_pmid_310434170) was compared between high- and low-risk groups using the Wilcoxon rank-sum test. Their association with the six stromal-cell senescence signature genes was evaluated by Pearson correlation analysis; gene pairs with R>0.1 and P<0.05 were considered statistically significant and selected for downstream visualization.

Chemotherapy response prediction

The “calcPhenotype” function of “oncoPredict” R (version 1.2) was used to predict the drug sensitivity for each patient based on bulk RNA sequencing (RNA-seq) results (21). The oncoPredict team has refreshed and shared the datasets used to fit the ridge regression model, incorporating data from the Genomics of Drug Sensitivity in Cancer (GDSC2) and Cancer Therapeutics Response Portal (CTRP) V2 databases (external, opens in a new tab. https://osf.io/c6tfx/).

Expression of the six-gene stromal senescence-associated gene signature

Gene-level expression values, represented as log₂(TPM +1), were extracted from the TCGA-LIHC dataset. Differential expression levels were assessed utilizing the two-sample Wilcoxon rank-sum test, with P<0.05 as the threshold for statistical significance. Protein-level log₂ tandem mass tag (TMT) intensities were obtained from the CPTAC-LIHC dataset, and comparisons between tumor and matched adjacent normal tissues were performed using the Wilcoxon rank-sum test, with P<0.05 as the threshold for statistical significance.

Analysis of spatial transcriptome databases

After loading the spatial objects with the R package Seurat (version 5.4.0), the data were normalized using the SCTransform function. The expression of the six stromal-senescence genes in each sample was plotted using the SpatialFeaturePlot function to visualize their spatial distribution across tissue regions.

Virtual screening

Protein binding pockets were predicted using PrankWeb (https://prankweb.cz/), and only those with a prediction score ≥0.7 were retained as principal potential binding sites. DrugCLIP (https://cn.bing.com/) employs contrastive learning to align representations of binding pockets and small molecules, enabling large-scale virtual screening. Screening was performed using the Vitas-M, ChemDiv, Enamine (SDF 202504), and Princeton BioMolecular Research libraries. Compounds were retained as potential drug-like candidates if they met the following criteria: DrugCLIP score >3, Vina docking score <−6.5, molecular weight 200–450 Da, logP <4.5, and quantitative estimate of drug-likeness (QED) ≥0.6.

Molecular docking of candidate drugs targeting DNASE1L3

The chemical compositions of candidate compounds were retrieved from the Vitas-M repository, while three-dimensional structural conformations of core proteins encoded by genes were obtained from the Protein Data Bank (PDB). Prior to analysis, protein structures underwent preprocessing protocols, including dehydration (removal of water molecules), hydrogenation (addition of hydrogens), and charge assignment. To elucidate binding mechanisms and interaction strengths, molecular docking was performed using AutoDock, aligning the candidate compounds with DNASE1L3 and other critical target proteins. Target protein preparation involved ligand and solvent elimination using PyMOL software. The docking outcomes were visualized through two-dimensional representations generated by AutoDock, depicting compound-protein interaction profiles.

Molecular dynamics (MD) simulation

MD simulations were conducted using GROMACS 2022, employing General AMBER Force Field (GAFF) parameters for ligands and AMBER14SB/TIP3P for protein systems. The solvated protein-ligand complex was simulated under periodic boundary constraints, with hydrogen atom-containing bonds constrained via LINCS (2.0 fs time step). Long-range electrostatics were calculated using Particle Mesh Ewald with a 1.2-nm cutoff, while non-bonded interactions were updated every 10 steps using a 1.0-nm cutoff. System equilibration was achieved at 298 K and 1 bar via a V-rescale thermostat and Berendsen barostat, followed by a 100-ns production phase. Simulation trajectories were captured at 10 ps intervals, and subsequent analyses were performed using VMD and PyMOL software suites. Two-dimensional interaction diagrams were generated through AutoDock to visualize compound-protein binding modes.

Correlation analysis

The gene sets were obtained from the MSigDB resource, including SAUL_SEN_MAYO, REACTOME_DNA_DAMAGE_TELOMERE_STRESS_INDUCED_SENESCENCE, REACTOME_ONCOGENE_INDUCED_SENESCENCE, REACTOME_OXIDATIVE_STRESS_INDUCED_SENESCENCE, and REACTOME_CELLULAR_SENESCENCE pathways. Their association with the senescence Area Under the Curve Cell (AUCell) scores was evaluated by Pearson correlation analysis to characterize the specific senescence mechanisms captured by our senescence-related score.

Statistical analyses

Data were analyzed via R (version 4.4.0) (https://www.r-project.org/). Kaplan-Meier estimator was employed to generate survival plots, and between-group differences were evaluated with the log-rank test. scRNA-seq and transcriptome data were deemed as statistically significance at a threshold of P or adjusted P value <0.05.

Results

Screening of senescence-related DEGs in HCC at the single-cell level

As shown in Figure 1, the flowchart summarizes the overall workflow. We integrated scRNA-seq data from four HCC datasets: GSE151530 (n=32), GSE149614 (n=13), GSE189903 (n=16), and GSE156625 (n=43). After PCA dimensionality reduction, batch effects were removed with Harmony, and the top 15 principal components were used for UMAP visualization, yielding 16 clusters (Figure 2A). These clusters were further categorized into five cell subtypes of epithelial cells, T cells, myeloid cells, B cells, and stromal cells via established HCC-specific marker genes (Figure 2B). Figure 2C displays the exact cell counts per subtype, and Figure 2D lists the top five marker genes for each. In total, 168,110 cells were successfully classified, and 1,150 HCC marker genes were identified. To assess the impact of senescence in HCC, we compiled 866 senescence genes. The Venn diagram revealed 93 overlapping DEGs between senescence genes and cell-type markers (Figure 2E). Heatmaps (Figure 2F) and bubble plots (Figure 2G) showed the expression of these 93 senescence-related DEGs across cell types. Stromal cells contained the largest set (33 genes), followed by myeloid cells (30 genes), epithelial cells (27 genes), T cells (3 genes), and none in B cells. These findings clarify the cell-type-specific roles of senescence genes in HCC (Table S2).

Figure 1 The workflow of this study. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. AUC, area under the curve; CPTAC, Clinical Proteomic Tumor Analysis Consortium; DEG, differentially expressed gene; GEO, Gene Expression Omnibus; GSEA, gene set enrichment analysis; K-M, Kaplan-Meier; LIHC, Liver Hepatocellular Carcinoma; M, metastasis; N, node; scRNA, single-cell RNA; T, tumor; TCGA, The Cancer Genome Atlas; UMAP, Uniform Manifold Approximation and Projection.

Figure 2 Differential expression profiling of senescence-related genes in HCC at the single-cell level. (A) UMAP visualization depicting 16 distinct cell clusters. (B) UMAP representation of annotated cell populations based on signature gene expression. (C) Quantification of identified cell type abundances determined by marker gene profiles. (D) Bubble plot illustrating the top 5 discriminatory marker genes for each characterized cell subtype. (E) Venn diagram depicting the intersection of senescence-related DEGs and cell type-specific marker genes. (F) Expression heatmap visualizing the transcriptional profiles of 93 senescence-associated genes across annotated cell populations. (G) Bubble plot quantifying the relative expression levels of 93 senescence-related genes in each characterized cell subtype. DEG, differentially expressed gene; HCC, hepatocellular carcinoma; UMAP, Uniform Manifold Approximation and Projection.

Scoring the cell subtype through 93 senescence genes

Due to the key effect of senescence in tumor progression, we quantified cellular senescence-related scores via expression profiles of 93 senescence-related genes using the AUCell algorithm. Correlation analysis revealed that stromal senescence AUCell scores were strongly correlated with SAUL_SEN_MAYO (SASP; Pearson r=0.67) but weakly correlated with other senescence pathways (Figure S1A). The UMAP visualization depicted the distribution of senescence scores across individual cells (Figure 3A). Comparative analysis revealed differential senescence activities among cell subtypes, with myeloid cells, epithelial cells, and stromal cells exhibiting elevated senescence activity, particularly stromal cells displaying the highest abundance, while T cells and B cells demonstrated relatively lower activity (Figure 3B and Figure S1B). Subsequently, cells were grouped into high-senescence and low-senescence cohorts through the median threshold of senescence scores. The UMAP plot vividly illustrated the single-cell-level distribution of these two groups, highlighting the preferential enrichment of high-senescence cells within myeloid, epithelial, and stromal populations (Figure 3C). The volcano plot revealed 5,608 upregulated and 1,598 downregulated TCGA-LIHC DEGs between HCC and normal tissues, as defined by the transcriptome profiles (Figure 3D). A total of 3,112 DEGs were identified within the high-senescence group. Table S3 shows the intersection of senescence-related genes with SASP genes, validated senescence regulators, and stromal markers. A total of 302 stromal-cell-specific marker genes were subsequently retrieved. Next, 63 overlapping DEGs were extracted by intersecting the DEGs from the TCGA dataset, the DEGs in the high-senescence group, and the stromal-cell marker genes (Figure 3E). Functional analyses were employed to enclose the biological functions. DEGs collectively depict a microenvironment-centric program through integrin-mediated regulation of the PI3K-AKT signaling axis, ultimately promoting tumor migration, invasion, and matrix remodeling (Figure 3F,3G). Consequently, these 63 shared DEGs were prioritized for further mechanistic investigations.

Figure 3 Quantification of cellular senescence scores based on senescence gene expression profiles. (A) UMAP visualization of cellular senescence enrichment scores calculated using the AUCell algorithm. (B) Violin plots depicting the distribution of senescence scores across characterized cell types. (C) UMAP clustering of cells stratified into high- and low-senescence subpopulations. (D) Volcano plot identifying DEGs between HCC and adjacent tissues. (E) Venn diagram delineating shared and unique DEGs between HCC and adjacent tissues, specifically enriched in the high-senescence group and stromal cells. (F,G) Functional annotation results of GO and KEGG enrichment analyses for DEGs in high-senescence stromal cells. AUCell, Area Under the Curve Cell; BP, biological process; CC, cellular component; DEG, differentially expressed gene; GO, Gene Ontology; HCC, hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; LIHC, Liver Hepatocellular Carcinoma; MF, molecular function; TCGA, The Cancer Genome Atlas; UMAP, Uniform Manifold Approximation and Projection.

Establishment of prognostic features integrated by machine learning tools

Firstly, univariate Cox regression was deployed to assess the prognostic significance of the 63 overlapping DEGs for HCC patient survival. Then, LASSO and multivariate Cox analysis were used to screen out the six prominent genes (DNASE1L3, NDRG2, ADAM15, IGFBP3, MMP14, and ANGPT2), and calculate the coefficient of each gene in the formula of risk score (Figure 5A,5B). Risk score = (−0.10385614 × DNASE1L3) + (−0.01711672 × NDRG2) + (0.10080776 × ADAM15) + (0.01860617 × IGFBP3) + (0.01461132 × MMP14) + (0.23162563 × ANGPT2). The forest plot (Figure 4C) revealed that, when all six genes were simultaneously entered into the model, DNASE1L3 and ANGPT2 maintained independent significance, suggesting these genes may serve as the primary effectors driving prognostic prediction. These weights serve to quantify the relative significance of each variable in contributing to the predictive capacity of the model. As plotted in Figure 4D, the high-risk cohort had a heightened mortality rate, and a heatmap visualized the expression profiles of DNASE1L3, NDRG2, ADAM15, IGFBP3, MMP14, and ANGPT2 across two subgroups. A nomogram was plotted by the integration of the risk signature and diverse clinicopathological factors and verified to predict survival outcomes in HCC patients (Figure 4E). The predictive performance for 1-, 3-, and 5-year survival probabilities of nomogram was rigorously assessed using calibration curves, which exhibited strong agreement between predicted and observed survival rates (Figure 4F). Kaplan-Meier curve revealed that high-risk patients exhibited significantly worse OS (P<0.0001) than those low-risk patients (Figure 5A). This prognostic signature was applied to predict the long-term survival probabilities of HCC patients, with the model demonstrating robust accuracy in estimating 1-, 3-, and 5-year survival outcomes, achieving areas under the curve (AUCs) of 0.78, 0.73, and 0.68, respectively, in HCC (Figure 5B). To enhance the study’s robustness, the GSE14520 dataset was leveraged as a validation set. In the validation cohort, patients in high-risk group experienced significantly worse outcomes than those in the low-risk group (P=0.003), and the AUCs for 1-, 3-, and 5-year survival were 0.69, 0.68, and 0.76, respectively (Figure 5C,5D). We further assessed the prognostic correlation of individual genes with HCC outcomes and obtained Kaplan-Meier curves, revealing the significant associations of all six genes with patient survival (Figure 5E-5J). Collectively, this study systematically elucidates the translational potential and clinical utility of utilizing stromal-cell senescence-related genes for prognostic prediction in HCC, providing a preliminary yet validated framework for precision oncology applications.

Figure 5 Survival analysis and prognostic performance of the stromal senescence risk signature. (A) Kaplan-Meier survival analysis comparing OS outcomes between high- and low-risk patient cohorts (n=343, median split, P values). (B) Time-dependent ROC curves assessing the discriminatory power of risk scores for 1-, 3-, and 5-year OS prediction. (C) OS disparity between risk groups validated using the GSE14520 external dataset (n=242, median split, P values). (D) ROC curve analysis of risk score performance for 1-, 3-, and 5-year OS prediction in the GSE14520 dataset. (E-J) The relationship between six senescence genes in stromal cells and survival time of HCC patients (n=343, optimal cut-off, P values). AUC, area under the curve; HCC, hepatocellular carcinoma; OS, overall survival; ROC, receiver operating characteristic.

Figure 4 Development and validation of a consensus stromal senescence-related signature through machine learning integration. (A) LASSO Cox regression analysis of the association between deviance and log(λ). (B) LASSO Cox regression analysis of the association between coefficients of genes and log(λ). (C) Forest plot of independent prognostic values for six senescence-related genes. (D) Distribution of survival status and time in HCC patients stratified by risk score, accompanied by a heatmap depicting differential expression of six senescence-related genes between high- and low-risk groups. (E) Detailed specification of the nomogram for prognostic prediction. (F) Calibration curves evaluating the accuracy of 1-, 3-, and 5-year OS predictions. *, P<0.05. AIC, Akaike information criterion; HCC, hepatocellular carcinoma; LASSO, least absolute shrinkage and selection operator; M, metastasis; N, node; OS, overall survival; T, tumor.

Enrichment analyses and mutation landscape

The pronounced differences in survival outcomes between groups underscored the existence of marked genomic heterogeneity. Gene set enrichment analysis (GSEA) was deployed to evaluate the mechanisms driving these differences, with results summarized in Figure 6A. As plotted in Figure 6B-6E, the enriched gene sets spanned a broad spectrum of biological processes, metabolic reprogramming (energy metabolism, lipid metabolism, and xenobiotic metabolism), immune response, cell cycle regulation, cell proliferation, and metastasis. In addition to the Hallmark gene sets, we incorporated two additional senescence-related gene sets into the analysis. GSEA of the high-risk cohort revealed significant enrichment of SASP pathways (P<0.05), whereas cellular senescence pathways were not enriched (Figure S1C). Furthermore, GO and KEGG analyses confirmed the roles of DEGs in regulating these biological processes (Figure 6F,6G), which finally mechanistically linked the distinct genomic landscapes to the observed prognostic disparity between the two HCC patient cohorts. Somatic mutation profiling was conducted across the entire sample cohort and compared between risk groups. Analysis indicated that TP53, TTN, CTNNB1, and MUC16 emerged as the top quartet of frequently mutated genes (FMGs) in HCC, accounting for the highest mutational prevalence within the tumor samples (Figure 6H,6I). In the high-risk group, TP53 emerged as the gene with the highest mutational frequency, accounting for 41% of cases analyzed. Conversely, CTNNB1 was identified as the top FMG in the low-risk group, with a mutation prevalence of approximately 34% (Figure 6H,6I).

Figure 6 Functional profiling and pathway enrichment analysis of DEGs between risk groups. (A) GSEA reveals significantly enriched biological processes and pathways associated with DEGs between high- and low-risk cohorts. (B-E) Enrichment plots illustrating the transcriptional signatures of metabolic reprogramming, inflammation & immunity, tumor proliferation, and tumor metastasis. (F) GO functional annotation unravels the biological functions of DEGs. (G) KEGG pathway enrichment analysis deciphering the regulatory networks of DEGs. (H) Somatic mutation landscape of the high-risk group. (I) Somatic mutation landscape of the low-risk group. BP, biological process; CC, cellular component; DEG, differentially expressed gene; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; TMB, tumor mutational burden.

Immune checkpoint analysis and drug sensitivity analysis

The association between immune checkpoint genes and risk signature was further explored. The TNFRSF4, CD276, TNFRSF18, ENTPD1, LAIR1, LGALS9, TNFSF15, HAVCR2, CTLA4, TNFSF4, TNFSF9, TNFRSF8, TNFRSF9, CD86, NRP1, CD200, TNFRSF14, TIGIT, ICOS, PDCD1, CD200R1, CD44, CD28, TNFRSF25, CD27, CD274, CD160, IDO1, CD40, CD48, and LAG3 were significantly overexpressed in high-risk group (Figure 7A). Furthermore, Figure 7B shows that the six stromal-senescence genes positively correlated with the 31 immune-checkpoint molecules listed above. Expression of immune checkpoint genes may induce a low response state. These findings underscore that the risk signature comprising six stromal-cell senescence genes was deemed as a predictive tool for optimizing immunotherapy selection in HCC patients. Notably, patients with low-risk scores commonly exhibit favorable responses to SB505124, IAP_5620, and JAK1_8709 while those in high-risk score demonstrate heightened sensitivity to sepantronium bromide, AZD6738, MK-1775, Tozasertib, Paclitaxel, YK-4-279, and BPD-00008900l treatment (Figure 7C).

Figure 7 Immune checkpoint and chemotherapy implications based on risk stratification. (A) Differential expression of ICI molecules between risk patient groups. (B) Correlation analysis between six senescence-related genes and infiltrating immune checkpoints. (C) Drug sensitivity prediction and therapeutic vulnerability assessment for risk groups. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. IC₅₀, half-maximal inhibitory concentration; ICI, immune checkpoint inhibitor.

Cell localization and expression analyses of 6 modeling genes

Aiming to decipher the expression landscapes of six modeling genes across various stromal cell subtypes, a thorough single-cell transcriptomic analysis was implemented. DNASE1L3, NDRG2, ADAM15, IGFBP3, MMP14, and ANGPT2 were primarily expressed in stromal cells (Figure S1D). Figure 8A and Figure S1E revealed differential gene expression patterns. Specifically, DNASE1L3, ADAM15, IGFBP3, and ANGPT2 were primarily expressed in ECs, while NDRG2 and MMP14 were mainly expressed in fibroblasts. Furthermore, differential expression profiles of these six prognostic genes were validated between HCC and normal hepatic tissues via transcriptomic data from the TCGA cohort and proteomic data from the CPTAC cohort. TCGA analysis revealed that ADAM15, MMP14, and ANGPT2 had significantly increased expression within HCC tissues relative to normal counterparts, while DNASE1L3, NDRG2, and IGFBP3 demonstrated downregulated expression in tumor samples (Figure 8B). Consistent with the TCGA database, CPTAC analysis revealed that ADAM15 and MMP14 exhibited significantly upregulated expression, while DNASE1L3, NDRG2, and IGFBP3 demonstrated downregulated expression at the protein level, except that ANGPT2 was not detected (Figure 8C). Spatial transcriptomics analyses in HCC were consistent with the TCGA findings (Figure 8D).

Figure 8 Expression of the six stromal senescence-related signature genes. (A) UMAP visualizations showing the differential expression of DNASE1L3, ANGPT2, ADAM15, IGFBP3, NDRG2, and MMP14 in stromal cells. (B) TCGA-LIHC of gene expression (DNASE1L3, ANGPT2, ADAM15, IGFBP3, NDRG2, and MMP14). (C) CPTAC showing protein expression (DNASE1L3, ADAM15, IGFBP3, NDRG2, and MMP14) from liver cancer. (D) Spatial transcriptome analysis reveals the spatial distribution and expression patterns of the six stromal senescence signature genes in HCC. ****, P<0.0001. CAF, cancer-associated fibroblast; CPTAC, Clinical Proteomic Tumor Analysis Consortium; EC, endothelial cell; HCC, hepatocellular carcinoma; LIHC, Liver Hepatocellular Carcinoma; MC, mural cell; TCGA, The Cancer Genome Atlas; TMT, tandem mass tag; TPM, transcripts per million; UMAP, Uniform Manifold Approximation and Projection.

Structure-based prediction and dynamic validation of DNASE1L3-targeted candidates using AutoDock Vina and MD simulations

Recently, HCC treatment has advanced to integrate targeted therapy alongside immunotherapy. However, a subset of certain patients remains refractory to this combination therapy. DNASE1L3, identified as a prominently upregulated protein in the treatment-responsive cohorts (22), emerges as a prospective target to potentiate the efficacy of this combination therapy in HCC.

Structure-based virtual screening was performed to ascertain compounds with high binding affinity to the DNASE1L3 protein. Protein pockets were predicted with PrankWeb (https://prankweb.cz/), and only those scoring >0.7 were retained as principal potential binding sites (Table 1). Using the DrugCLIP platform, we docked the Vitas-M, ChemDiv, Enamine Screening Collection SDF 202504, and Princeton BioMolecular Research libraries into the DNASE1L3 active pocket to identify potential binders. Among all the screened compounds, STL167676, STK598132, and STK119225 exhibited favorable predicted affinities and drug-like properties toward DNASE1L3 (Table 2).

STL167676 exhibited the highest predicted affinity for the DNASE1L3 chain A and chain B binding pocket (centroid coordinates: x=46.4445, y=53.0234, z=−4.0798), with DrugClip score of 3.01 and Vina score of −9.518 kcal/mol (Tables 1,2). Docking analysis revealed that STL167676 stably binds the DNASE1L3 protein through multiple stabilizing interactions. These included a hydrogen bond with residue ARG-132, in addition to π-anion and π-alkyl hydrophobic interactions involving GLU-159, ASP-273, LYS-226, ARG-95, and PRO-158. Additional van der Waals interactions were detected with TYR-98, THR-156, and THR-224 (Figure 9A). To probe the dynamic stability of the DNASE1L3-ligand complex, a 100-ns all-atom MD simulation was executed in an explicit TIP3P water environment. The protein was parameterized using the AMBER14SB force field, while the ligand was parameterized with GAFF. During the simulation trajectory, the ligand consistently retained a compact and embedded orientation within the binding cavity. Interaction mapping analyses revealed persistent interactions with pivotal residues such as ARG-132, GLU-159, ASP-273, and LYS-226. The trajectory analysis showed low root mean square deviation (RMSD) values (Figure 9D), a stable radius of gyration (Rg; Figure 9E), and limited flexibility at the binding interface as indicated by root mean square fluctuation (RMSF; Figure 9F). The solvent-accessible surface area (SASA; Figure 9G) and the center-of-mass distance between the ligand and the protein (Figure 9H) exhibited consistent behavior throughout the entire simulation, indicating a stable and prolonged interaction.

Figure 9 Structure-based molecular docking and dynamic simulations reveal high-affinity ligands targeting DNASE1L3 receptor. (A-C) Binding mode visualization of DNASE1L3 with STL167676, STK598132, and STK119225: molecular docking was conducted using AutoDock Vina, revealing favorable binding energies indicative of spontaneous interactions. (D-H) Validation through 100-ns all-atom MD simulation of DNASE1L3-STL167676 complex. (D) RMSD values maintained low magnitude, confirming structural equilibration. (E) Rg analysis suggests a compact and stable protein conformation. (F) RMSF profiles demonstrate limited flexibility at the ligand-binding interface. (G) SASA exhibits consistency over time, indicating unperturbed ligand-protein accessibility. (H) Stable center-of-mass separation between ligand and protein sustained throughout simulation, validating persistent binding. MD, molecular dynamics; Rg, radius of gyration; RMSD, root mean square deviation; RMSF, root mean square fluctuation; SASA, solvent-accessible surface area.

Two additional top-scoring ligands, STK598132 docked within pocket 2 (DrugClip score =3.33, Vina score =−7.010 kcal/mol; Table 2) and STK119225 docked within pocket 3 (DrugClip score =3.18, Vina score =−6.671 kcal/mol; Table 2), underwent the same molecular docking and MD protocol. Both compounds engaged with critical residues and exhibited comparable dynamic profiles, as evidenced by stable RMSD and SASA values and persistent binding-site contacts over 100 ns (Figure 9B,9C and Figure S2). They displayed favorable interaction profiles, highlighting the druggability potential of targeting DNASE1L3.

Discussion

Accurate prediction of patient survival outcomes and treatment response is paramount for developing tailored treatment strategies in precision oncology. Cellular senescence contributes to modulating cell function and is now recognized as a defining characteristic of solid tumors, colloquially termed a “hallmark of cancer” (16).

Notably, stromal cells undergo senescence-driven transformations that critically fuel the development of tumor. Senescent cancer-associated fibroblasts (CAFs) actively promote tumor progression and invasion by secreting cytokines and extracellular vesicles (23-25). In HCC, hepatic stellate cells (HSCs)—fibroblast-like cells within the TME—function analogously to CAFs in other malignancies. Diet-induced or genetic obesity alters gut microbiota, elevating deoxycholic acid (DCA) levels. Through enterohepatic circulation, DCA triggers DNA damage-induced SASP in HSCs. Senescent HSCs, characterized by SASP activation, secrete a diverse array of factors that collectively promote tumor proliferation, invasiveness, and immune evasion, thereby driving HCC aggressiveness (26). Furthermore, prior investigations have demonstrated that senescent ECs display upregulated expression of ICAM-1 and VCAM-1, alongside downregulation of vascular endothelial (VE)-cadherin, which collectively enhance vascular permeability and facilitate cancer cell dissemination (27-29). Additionally, chemokines secreted by senescent endothelium, such as IL-6 and CXCL11, directly promote tumor cell proliferation and aggressiveness (30). Notably, these changes in ECs also recruit immune cells into the TME, eliciting either antitumor or protumorigenic effects (29,31-33). Moreover, radiation-induced EC senescence has been demonstrated to stimulate cellular proliferation and accelerate the advancement of cancer (29,33). Therefore, targeting stromal cell senescence emerges as a compelling therapeutic strategy to improve patient survival outcomes and modulate tumor-immune interactions within the microenvironment.

scRNA-seq information from 104 HCC patients was retrieved from GSE149614, GSE151530, GSE156625, and GSE189903 datasets following stringent quality control screening. After comprehensive annotation, 16 cell clusters were successfully categorized into five distinct cell types. To assess the prognostic potential of senescence in HCC, we leveraged a well-established senescence-associated gene signature (34), which—unlike other available signatures (35,36)—comprehensively captures transcriptomic alterations specific to senescence. Then, 93 overlapping DEGs were selected between the senescence genes and marker genes for further cellular-level validation in HCC. Each cell was assigned a senescence score through these 93 genes. Our transcriptomic analysis unequivocally demonstrates that stromal cells exhibit the most pronounced senescence signatures across all cellular populations, highlighting their central role in the senescence landscape of HCC. To clarify whether the signature reflects the presence of senescent cells or the systemic effects of their secretome, we correlated stromal senescence AUCell scores with established gene sets. The results showed strong correlation with SAUL_SEN_MAYO (SASP; Pearson r=0.67) but weak correlation with other pathways. Consistently, GSEA of the high-risk cohort revealed significant enrichment of SASP pathways (P<0.05), whereas cellular senescence pathways were not enriched. Therefore, the model predicts survival based on the systemic effects of the secretome.

The 3,112 DEGs were determined by comparative analysis between senescence groups. Next, 63 DEGs were overlapped from 3,112 DEGs in the high-senescence group, 302 stromal cell marker genes, and 7,206 TCGA-LIHC DEGs. This 63-gene panel was then utilized to develop a stromal cell senescence-related signature. Multiple analyses screened out the six prominent genes (DNASE1L3, NDRG2, ADAM15, IGFBP3, MMP14, and ANGPT2) and formulated the risk score. Although ANGPT2 was not detected in the proteomic analysis (likely due to technical limitations of the detection platform), a study by Xie et al. demonstrated high protein expression of ANGPT2 in HCC tissues using immunohistochemistry (IHC), immunoblotting, and enzyme-linked immunosorbent assay (ELISA) (37). Consequently, we retained this gene in the prognostic model. The classification of HCC patients was applied for patients via the median risk score. Moreover, this signature effectively differentiated the survival outcomes between the two HCC patient cohorts and addressed the inherent limitations of conventional TNM staging in clinical practice. Remarkably, in the high-risk group, metabolic reprogramming shifts cells toward anaerobic glycolysis while suppressing lipid and xenobiotic metabolism; concurrently, it fuels immune-inflammatory activation, ultimately accelerating tumor proliferation and metastasis. The divergent mutation landscapes, immune-related gene signatures, and drug-sensitivity profiles between the two groups provide crucial clues for selecting personalized targeted therapies for individual HCC patients.

Our prognostic model integrates six functionally diverse genes whose roles in HCC and other cancers pathogenesis align well with current mechanistic understanding. Guo et al. demonstrated that DNASE1L3 inhibits tumor angiogenesis via impairing the senescence-associated secretory phenotype in response to DNA damage stress in HCC progression (38). According to Lee et al., NDRG2 can inhibit extracellular matrix-based, Rho-driven tumor cell invasion and migration and thereby plays important roles in suppressing tumor metastasis in HCC (39). Kuefer et al. demonstrated that ADAM15 is generally overexpressed in adenocarcinoma and is highly associated with metastatic progression of prostate and breast cancers (40). Li et al. revealed that activation of the IGFBP3-AKT/mTOR axis promotes HCC tumorigenesis and lung metastasis (41). According to Shi et al., MMP14 was identified as a novel cis-interacting partner for CD48. This interaction upregulates immunosuppressive genes in HCC (42). Zhu et al. revealed that upregulation of ANGPT2 induced angiogenesis and conferred sorafenib resistance to HCC cells (43).

The treatment of HCC has recently evolved to incorporate targeted therapy alongside immunotherapy (44). Nevertheless, certain patients continue to exhibit inherent resistance to this combination therapy (45). Wang et al.’s analysis of differentially expressed proteins in the combination therapy involving sorafenib and PD-1 with DNASE1L3 being identified as one of the most notably upregulated proteins in the treatment-responsive group (22). DNASE1L3 shows significant promise as a stromal senescence-related candidate for enhancing combination therapy efficacy in HCC. To determine its therapeutic potential, we deployed structure-based virtual screening, identifying STL167676, STK598132, and STK119225 as leading candidates. They showed robust binding affinity and conformational stability across simulations, validating their candidacy as DNASE1L3-targeted candidates.

The novel contributions of our study are multifaceted. First, we systematically characterized stromal senescence-associated genes in HCC at single-cell resolution and, by integrating scRNA-seq with bulk RNA-seq data, constructed a novel stromal-senescence signature. Second, we innovatively incorporated these findings with clinical parameters to develop a clinically actionable nomogram, which seems a practical and user-friendly tool enabling personalized prognosis assessment for HCC patients by stromal senescence-associated gene signature. This combination of molecular and clinical metrics significantly elevates the translational relevance and practical utility of our research. Thirdly, we performed structure-based virtual screening to find that STL167676, STK598132, and STK119225 hold promise as novel DNASE1L3-targeted candidates to overcome resistance to sorafenib plus anti-PD-1 monoclonal antibody (mAb) combination therapy.

However, some inherent limitations warrant acknowledgment. Firstly, the stromal cell senescence-driven risk prediction model was exclusively developed using retrospective data sourced from public repositories. Consequently, prospective validation across large-scale, multicenter HCC cohorts is imperative to establish its clinical applicability. Secondly, the public cohorts employed in this study (TCGA, GEO, and CPTAC) exhibit inherent heterogeneity in patient populations [e.g., diverse etiologies including HBV, HCV, alcohol, and non-alcoholic fatty liver disease (NAFLD)] and treatment regimens (surgical resection, transarterial chemoembolization, systemic therapy). These variations may introduce confounding factors that influence stromal senescence patterns and prognostic outcomes. Future studies should stratify analyses by specific etiologies and treatment modalities to validate the signature’s robustness across homogeneous subgroups. Thirdly, the lack of functional validation experiments directly testing the role of DNASE1L3 or the six-gene signature. Fourthly, our investigations primarily focused on the prognostic implications of this signature, so further mechanistic studies are essential to elucidate how stromal cell senescence-related gene alterations contribute to HCC progression.

Conclusions

In conclusion, this study developed a stromal senescence-associated gene signature capable of robustly predicting survival outcomes in HCC through the integration of scRNA-seq and transcriptome RNA-seq data. It holds the potential to steer future investigations into stromal senescence mechanisms in HCC and facilitate personalized selection of chemotherapy agents and immune checkpoint inhibitors (ICIs) for patients. STL167676, STK598132, and STK119225 hold promise as novel DNASE1L3-targeted candidates to overcome resistance to sorafenib plus anti-PD-1 mAb combination therapy.

Acknowledgments

This study analyzed publicly available datasets. The data that support the results of the current study are available on TCGA (https://portal.gdc.cancer.gov/), GEO websites (http://www.ncbi.nlm.nih.gov/geo), and CPTAC (https://cptac-data-portal.georgetown.edu/). We extend our gratitude to the curators and teams behind the public databases that made this work possible.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2894/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2894/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2894/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Craig AJ, von Felden J, Garcia-Lezana T, et al. Tumour evolution in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol 2020;17:139-52. [Crossref] [PubMed]
3. Llovet JM, Kelley RK, Villanueva A, et al. Hepatocellular carcinoma. Nat Rev Dis Primers 2021;7:6. [Crossref] [PubMed]
4. Zheng R, Qu C, Zhang S, et al. Liver cancer incidence and mortality in China: Temporal trends and projections to 2030. Chin J Cancer Res 2018;30:571-9. [Crossref] [PubMed]
5. Beaufrère A, Caruso S, Calderaro J, et al. Gene expression signature as a surrogate marker of microvascular invasion on routine hepatocellular carcinoma biopsies. J Hepatol 2022;76:343-52. [Crossref] [PubMed]
6. Tian Z, Song J, She J, et al. Constructing a disulfidptosis-related prognostic signature of hepatocellular carcinoma based on single-cell sequencing and weighted co-expression network analysis. Apoptosis 2024;29:1632-47. [Crossref] [PubMed]
7. Yang D, Zhou Y, Zhang Y, et al. Comprehensive analysis of scRNA-Seq and bulk RNA-Seq data reveals dynamic changes in tumor-associated neutrophils in the tumor microenvironment of hepatocellular carcinoma and leads to the establishment of a neutrophil-related prognostic model. Cancer Immunol Immunother 2023;72:4323-35. [Crossref] [PubMed]
8. Calcinotto A, Kohli J, Zagato E, et al. Cellular Senescence: Aging, Cancer, and Injury. Physiol Rev 2019;99:1047-78. [Crossref] [PubMed]
9. Campisi J. Aging, cellular senescence, and cancer. Annu Rev Physiol 2013;75:685-705. [Crossref] [PubMed]
10. Liu XL, Ding J, Meng LH. Oncogene-induced senescence: a double edged sword in cancer. Acta Pharmacol Sin 2018;39:1553-8. [Crossref] [PubMed]
11. Wang B, Kohli J, Demaria M. Senescent Cells in Cancer Therapy: Friends or Foes? Trends Cancer 2020;6:838-57. [Crossref] [PubMed]
12. Dou X, Fu Q, Long Q, et al. PDK4-dependent hypercatabolism and lactate production of senescent cells promotes cancer malignancy. Nat Metab 2023;5:1887-910. [Crossref] [PubMed]
13. Marin I, Boix O, Garcia-Garijo A, et al. Cellular Senescence Is Immunogenic and Promotes Antitumor Immunity. Cancer Discov 2023;13:410-31. [Crossref] [PubMed]
14. De Blander H, Morel AP, Senaratne AP, et al. Cellular Plasticity: A Route to Senescence Exit and Tumorigenesis. Cancers (Basel) 2021;13:4561. [Crossref] [PubMed]
15. Lee S, Schmitt CA. The dynamic nature of senescence in cancer. Nat Cell Biol 2019;21:94-101. [Crossref] [PubMed]
16. Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov 2022;12:31-46. [Crossref] [PubMed]
17. Anderson NM, Simon MC. The tumor microenvironment. Curr Biol 2020;30:R921-5. [Crossref] [PubMed]
18. Bussard KM, Mutkus L, Stumpf K, et al. Tumor-associated stromal cells as key contributors to the tumor microenvironment. Breast Cancer Res 2016;18:84. [Crossref] [PubMed]
19. Bardou P, Mariette J, Escudié F, et al. jvenn: an interactive Venn diagram viewer. BMC Bioinformatics 2014;15:293. [Crossref] [PubMed]
20. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 2012;16:284-287. [Crossref] [PubMed]
21. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021;22:bbab260. [Crossref] [PubMed]
22. Wang J, Chen Y, Xu Y, et al. DNASE1L3-mediated PANoptosis enhances the efficacy of combination therapy for advanced hepatocellular carcinoma. Theranostics 2024;14:6798-817. [Crossref] [PubMed]
23. Krtolica A, Parrinello S, Lockett S, et al. Senescent fibroblasts promote epithelial cell growth and tumorigenesis: a link between cancer and aging. Proc Natl Acad Sci U S A 2001;98:12072-7. [Crossref] [PubMed]
24. Takasugi M, Okada R, Takahashi A, et al. Small extracellular vesicles secreted from senescent cells promote cancer cell proliferation through EphA2. Nat Commun 2017;8:15729. [Crossref] [PubMed]
25. Park HS, Kim SY. Endothelial cell senescence: A machine learning-based meta-analysis of transcriptomic studies. Ageing Res Rev 2021;65:101213. [Crossref] [PubMed]
26. Yoshimoto S, Loo TM, Atarashi K, et al. Obesity-induced gut microbial metabolite promotes liver cancer through senescence secretome. Nature 2013;499:97-101. [Crossref] [PubMed]
27. Krouwer VJ, Hekking LH, Langelaar-Makkinje M, et al. Endothelial cell senescence is associated with disrupted cell-cell junctions and increased monolayer permeability. Vasc Cell 2012;4:12. [Crossref] [PubMed]
28. Kandhaya-Pillai R, Miro-Mur F, Alijotas-Reig J, et al. TNFα-senescence initiates a STAT-dependent positive feedback loop, leading to a sustained interferon signature, DNA damage, and cytokine secretion. Aging (Albany NY) 2017;9:2411-35. [Crossref] [PubMed]
29. Wieland E, Rodriguez-Vita J, Liebler SS, et al. Endothelial Notch1 Activity Facilitates Metastasis. Cancer Cell 2017;31:355-67. [Crossref] [PubMed]
30. Wang L, Lankhorst L, Bernards R. Exploiting senescence for the treatment of cancer. Nat Rev Cancer 2022;22:340-55. [Crossref] [PubMed]
31. Bloom SI, Islam MT, Lesniewski LA, et al. Mechanisms and consequences of endothelial cell senescence. Nat Rev Cardiol 2023;20:38-51. [Crossref] [PubMed]
32. Behmoaras J, Gil J. Similarities and interplay between senescent cells and macrophages. J Cell Biol 2021;220:e202010162. [Crossref] [PubMed]
33. Gabai Y, Assouline B, Ben-Porath I. Senescent stromal cells: roles in the tumor microenvironment. Trends Cancer 2023;9:28-41. [Crossref] [PubMed]
34. Fridman AL, Tainsky MA. Critical pathways in cellular senescence and immortalization revealed by gene expression profiling. Oncogene 2008;27:5975-87. [Crossref] [PubMed]
35. Saul D, Kosinsky RL, Atkinson EJ, et al. A new gene set identifies senescent cells and predicts senescence-associated pathways across tissues. Nat Commun 2022;13:4827. [Crossref] [PubMed]
36. Casella G, Munk R, Kim KM, et al. Transcriptome signature of cellular senescence. Nucleic Acids Res 2019;47:7294-305. [Crossref] [PubMed]
37. Xie JY, Wei JX, Lv LH, et al. Angiopoietin-2 induces angiogenesis via exosomes in human hepatocellular carcinoma. Cell Commun Signal 2020;18:46. [Crossref] [PubMed]
38. Guo D, Ma D, Liu P, et al. DNASE1L3 arrests tumor angiogenesis by impairing the senescence-associated secretory phenotype in response to stress. Aging (Albany NY) 2021;13:9874-99. [Crossref] [PubMed]
39. Lee DC, Kang YK, Kim WH, et al. Functional and clinical evidence for NDRG2 as a candidate suppressor of liver cancer metastasis. Cancer Res 2008;68:4210-20. [Crossref] [PubMed]
40. Kuefer R, Day KC, Kleer CG, et al. ADAM15 disintegrin is associated with aggressive prostate and breast cancer disease. Neoplasia 2006;8:319-29. [Crossref] [PubMed]
41. Li X, Zhang CC, Lin XT, et al. Elevated expression of WSB2 degrades p53 and activates the IGFBP3-AKT-mTOR-dependent pathway to drive hepatocellular carcinoma. Exp Mol Med 2024;56:177-191. [Crossref] [PubMed]
42. Shi G, Xiao Y, Li Z, et al. CD48 is a novel immune checkpoint on tumour-associated macrophages in hepatocellular carcinoma. Gut 2026;gutjnl-2025-336744.
43. Zhu J, Wu Y, Yu Y, et al. MYBL1 induces transcriptional activation of ANGPT2 to promote tumor angiogenesis and confer sorafenib resistance in human hepatocellular carcinoma. Cell Death Dis 2022;13:727. [Crossref] [PubMed]
44. Rimassa L, Finn RS, Sangro B. Combination immunotherapy for hepatocellular carcinoma. J Hepatol 2023;79:506-15. [Crossref] [PubMed]
45. Chiang HC, Lee YC, Chang TT, et al. Real-World Effectiveness of Sorafenib versus Lenvatinib Combined with PD-1 Inhibitors in Unresectable Hepatocellular Carcinoma. Cancers (Basel) 2023;15:854. [Crossref] [PubMed]