The development and validation of a novel senescence-related long-chain non-coding RNA (lncRNA) signature that predicts prognosis and the tumor microenvironment of patients with hepatocellular carcinoma

Introduction

With a 5-year survival rate of only 18%, liver cancer is the second most lethal tumor after pancreatic cancer (1). Hepatocellular carcinoma (HCC) accounts for between 85% and 90% of primary liver cancers (2). Surgery, surgical ablation, and orthotopic liver transplantation are the most popular locoregional treatment options for HCC (3). Despite the approval of several tyrosine kinase inhibitors (TKIs) for first- and second-line systemic treatment, the overall survival (OS) of advanced HCC is poor and has not improved in the last decade (4). More recently, the FDA approved immune checkpoint inhibitors (ICIs) for HCC treatment. Unfortunately, even with the combination of anti-CTLA-4 and anti-PD-L1, the objective response rate (ORR) was only 20.1% (5). Due to the lack of appropriate patient-selective biomarkers, most HCC patients do not respond to ICI therapies (6). Therefore, new effective biomarkers must be developed in order to improve patient outcomes.

The known traditional biomarkers of HCC, such as Alpha-fetoprotein (AFP), Glypican-3, Osteopontin and Des-γ-Carboxy Prothrombin, despite common use, are characterized by poor sensitivity and specificity (7,8). It is now believed that abnormal miRNA and long non-coding RNAs (lncRNAs) expression plays an important role in the oncogenic role of HCC (9,10). For example, miR-224 concentrations are higher in the serum of HCC patients, and the higher the concentration, the shorter the survival time (11). A high level of LINC00958 aggravated HCC lipogenesis and progression through sponging miR-3619-5p, further upregulating the hepatoma-derived growth factor expression (12). Although several lncRNA prognostic models for HCC have been developed (13-18), the potential immune microenvironment predictive ability and therapeutic guidance of the models, especially for immunotherapy, have yet to be evaluated.

Cellular senescence is a protective mechanism of tissue homeostasis that has long been thought to be a cancer-preventive factor (19). It results in irreversible inhibition of proliferation, which, as a supplement to programmed cell death, inactivates cells and eventually eliminates diseased or defective cells (20). However, there is growing evidence that most senescent cells in the tumor microenvironment (TME), regardless of cell origin, participate in the shaping of cancer phenotypes through the so-called senescence-associated secretory phenotype (SASP) (21). The first evidence that the SASP promotes the proliferation of cancer cells came from studying co-cultures of preneoplastic epithelial cell lines with senescent WI-38 lung fibroblasts. When cultured together in vitro, the proliferation rate of preneoplastic epithelial cells was quicker than co-cultures with non-senescent fibroblasts (22). Co-injection of breast cancer cells with senescent fibroblasts in nude mice also results in faster-growing tumors (23). Interleukin (IL)-6 is a known SASP factor that binds to the IL-6 receptor of cancer cells through senescence-related paracrine effects and then activates STAT3, which transcribes oncogenes including the complex of CyclinD1 and mammalian target of rapamycin (mTOR) protein (21,24). Thus, with the increasing understanding of cellular senescence, it has recently been included as a significant hallmark of cancer (20). To date, however, there is a general lack of understanding of the underlying roles and mechanisms of genes involved in the senescence regulatory network of the TME of HCC.

LncRNAs participate in a variety of biological processes, and accumulating evidence indicates that various HCC-related lncRNAs are abnormally expressed in tumor tissues and participate in cancer phenotypes by binding to DNA, RNA, or proteins or by encoding small peptides (25). Similarly, lncRNAs are also significantly involved in regulatory processes in senescent cells (26). A recent transcriptional sequence study showed that lncRNA PURPL is one of the most abundant transcripts of lncRNA in senescent cells induced by radiation or drug therapy. PURPL reduces the stability of p53 by competitively binding to MYBBP1A with p53, thus maintaining a low concentration of intracellular p53 (27). Since PURPL promotes the tumorigenesis of colon cancer cells by negatively regulating p53 as a survival factor, it may also contribute to the survival phenotype of senescent cells (26). In addition, short peptide PINT87aa, which is encoded by LINC-PINT, was demonstrated to mediate cell cycle arrest, cellular senescence, and mitophagy in HCC cells (28). In view of this, identifying lncRNA transcriptional changes helps to define senescence and cancer characteristics.

Herein, we constructed and validated an HCC senescence-related lncRNA signature (HCCSenLncSig) and demonstrated that it can predict the prognosis of HCC patients. In addition, we developed a nomogram that included HCCSenLncSig and clinical factors and further compared gene enrichment, mutations, immune cell infiltration, and potential response to targeted therapy and immunotherapy in the HCCSenLncSig high- and low-risk groups. The present study highlights the regulatory network of cellular senescence, which is expected to improve the efficacy of individualized treatment for HCC. We present the following article in accordance with the TRIPOD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-22-3348/rc).

Methods

Dataset and sample extraction

Figure 1 depicts the flow chart of the present study. RNA-sequencing data (RNA-seq), clinical characteristics, and mutation data of HCC were obtained from The Cancer Genome Atlas-Liver Hepatocellular Carcinoma (TCGA-LIHC) database (https://portal.gdc.cancer.gov/). Patients lacking complete follow-up information, surviving less than 30 days, or lacking complete clinicopathological information will be excluded. The 279 cellular senescence genes were obtained from the online database CellAge (http://genomics.senescence.info/cells), which was developed after scouring the literature for gene manipulation experiments that caused cells to induce or inhibit cellular senescence (29). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Figure 1 Flowchart of the present study. TCGA, The Cancer Genome Atlas; RNA-seq, RNA sequencing; HCC, hepatocellular carcinoma; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNAs, long non-coding RNAs; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; ROC, receiver operating characteristic; TMB, tumor mutation burden.

Cellular senescence genes with differential expression in HCC and normal tissues

The differentially expressed genes (DEGs) involved in cellular senescence between tumor and normal tissues were identified with a log₂ fold change absolute value greater than 1 and a false discovery rate (FDR) of <0.05. The R package “limma” was used to visualize the DEGs. The DEGs were then analyzed by Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis.

Identifying an HCCSenLncSig

The absolute values of Pearson correlation coefficient >0.3 and P<0.001 were used as thresholds in the establishment of a senescence-related messenger RNA (mRNA)-lncRNA co-expression network to identify the senescence-related lncRNAs. The lncRNA-mRNA co-expression network was visualized using a Sankey diagram generated by the R package “ggalluvial” and Cytoscape software (version 3.7.2). We first used univariate Cox regression analysis to identify lncRNAs associated with HCC prognosis, and then included them in multivariate Cox regression analysis to calculate risk scores. The risk score of each patient was calculated by the following formula: risk score = explncRNA1 × coef lncRNA1 + explncRNA2 × coef lncRNA2 + ... + explncRNAi × coef lncRNAi. Then, HCC patients were classified into high- and low-risk groups based on median risk scores. We termed this lncRNA predictive signature HCCSenLncSig. To determine whether the HCCSenLncSig risk score is an independent prognostic indicator, demographic data were included in the univariate and multivariate Cox regression analyses. The risk score level was used to investigate the distribution of survival status. The receiver operating characteristic (ROC) curve was used to assess the accuracy of the risk model. To visualize clinicopathological variables in high- and low-risk groups, the R package “pheatmap” was used. The distribution of patients with varying risk scores was assessed using principal component analysis (PCA), which was then visualized using the R package “scatterplot3d.” To test whether HCCSenLncSig is generally applicable to predict the prognosis of HCC, we randomly divided TCGA-LIHC patients into 2 cohorts and repeated the grouping 1,000 times, and randomly selected 1 of the cohorts for internal validation. The K-M curves and time-dependent ROC curves of the two internal validation sets grouped according to the HCCSenLncSig median risk score were tested.

Construction of the nomogram

By combining risk scores with clinicopathological features, we constructed a nomogram using the R package “rms” that predicts HCC patients’ survival at 1, 3, and 5 years. A calibration curve was used to determine whether the predicted survival rate matched the actual survival rate.

Functional enrichment analysis of the senescence-related lncRNA predictive signature

According to the prognostic model, gene set enrichment analysis (GSEA) was performed for the high- and low-risk groups. Significant biological processes and pathways were enriched when the nominal (NOM) P<0.05 and FDR q<0.25 were used. The reference files from the Molecular Signatures Database (MSigDB) were c2.cp.kegg.v7.4.symbols.gmt and c5.go.v7.4 symbols.gmt (http://software.broadinstitute.org/gsea/msigdb/index.jsp). The “ggplot2” R package was used to visualize the results.

Analysis of somatic mutation data and tumor mutation burden (TMB)

The TMB was calculated as the number of somatic, coding, base replacement, and insert-deletion mutations discovered per megabase of the genome using non-synonymous and code-shifting indels and a 5% detection limit. The R package “maftools” (30) was used to count the number of somatic non-synonymous point mutations in each sample. In addition, TMB was compared between high- and low-risk groups, and survival curves for TMB and risk score integration were plotted.

Estimation of immune cell infiltration

Single-sample gene set enrichment analysis (ssGSEA) was performed using the R package “GSVA” to assess the infiltration fraction of 16 immune cells and the activity of 13 immune-related pathways (31). In order to cross-reference to the ssGSEA results, the CIBERSORT algorithm (32) was used to estimate the infiltration proportions of 22 immune cell types in HCC samples. We retained samples with P<0.05 for further investigation. To determine whether there was a significant difference in the proportion of immune cells between the low- and high-risk groups, the Wilcoxon rank-sum test was used.

Potential relationship between the HCCSenLncSig and immunotherapy, chemotherapy, and targeted therapy

First, the differential expression of 40 immune checkpoints was compared between high- and low-risk HCCSenLncSig groups. The tumor immune dysfunction and exclusion module (TIDE, http://tide.dfci.harvard.edu/), which predicts treatment responses to anti-PD1 and anti-CTLA4 therapies based on pretreatment transcript expression profiles, was then used to distinguish between the potential treatment responses of high- and low-risk groups. In addition, to further assess the role of the HCCSenLncSig in predicting response to HCC treatment, we calculated the half-maximal inhibitory concentration (IC50) of chemotherapy and targeted therapy drugs commonly used in clinical HCC treatment. The Wilcoxon signed-rank test and the R package “pRRophetic” were used to compare the IC50 values between the high- and low-risk groups.

Statistical analysis

To compare the levels of expression of cell senescence-related DEGs in cancer and normal tissues, the Wilcoxon test was used. The Kaplan-Meier method and log-rank test were used to compare the OS between patients in the high- and low-risk groups. The “survivalROC” package was used to develop the 5 years’ time-dependent ROC curves, calculate the area under the curve (AUC), and determine the sensitivity and specificity of the model at the optimal cut-off value. The model is considered to have good performance when the AUC value is greater than 0.7. To compare differences between groups, the Kruskal-Wallis test was used. Clinical data was analyzed using the chi-squared test or Fisher’s exact test. The relationships between lncRNA expression, immune infiltration, and immune checkpoint gene expression were assessed using Spearman or Pearson correlation coefficients. All statistical analyses were performed using R software (version 4.1.2), and a corrected P value of less than 0.05 was considered statistically significant.

Results

Enrichment analysis of senescence-related DEGs in HCC

There were 424 HCC patients in total. After excluding those who lacked complete follow-up information, survived less than 30 days, or lacked complete clinicopathological data, 343 patients were eventually recruited for follow-up analysis. To determine whether the senescence-related genes are abnormally expressed in tumor tissues, we first compared the expression levels of 279 genes from the CellAge database between HCC tumor and normal tissues (available in https://cdn.amegroups.cn/static/public/atm-22-3348-1.xlsx). A total of 71 DEGs were obtained, of which 27 were down-regulated and 44 were up-regulated in HCC tumor tissues (Figure 2A,2B, available in https://cdn.amegroups.cn/static/public/atm-22-3348-2.xlsx). KEGG pathway analysis revealed that the 5 most enriched pathways were cellular senescence, human T-cell leukemia virus 1 infection, bladder cancer, endocrine resistance, and hepatitis B (Figure 2C). On the other hand, GO analysis revealed that the 5 most enriched categories were cell aging, aging, cellular response to chemical stress, epithelial cell proliferation, and regulation of DNA-binding transcription factor activity (Figure 2D). These findings suggested that the DEGs are primarily involved in cell senescence, tumorigenesis, cell proliferation, stress, and transcriptional activation.

Figure 2 Identification of senescence-related DEGs. (A) Heatmap of 71 senescence-related DEGs in normal and tumor HCC tissues. (B) Volcano plot of 279 senescence-related genes in normal and HCC tumor tissues. Pink dots represent genes that are up-regulated, while blue dots represent genes that are down-regulated. (C) KEGG analysis of senescence-related DEGs. (D) GO analysis of senescence-related DEGs. N, normal tissues; T, tumor tissues; FDR, false discovery rate; FC, fold change; PD-L1, programmed cell death ligand-1; PD-1, programmed cell death-1; BP, biological process; CC, cellular components; MF, molecular function; DEGs, differentially expressed genes; HCC, hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

Construction of the HCCSenLncSig

Pearson correlation analysis was used to identify 734 senescence-related lncRNAs (available in https://cdn.amegroups.cn/static/public/atm-22-3348-3.xlsx). Univariate Cox regression analysis revealed that 54 lncRNAs were associated with the prognosis of HCC (Figure 3A). Using multivariate Cox regression analysis, 8 senescence-related lncRNAs (AL117336.3, AC103760.1, FOXD2-AS1, AC009283.1, AC026401.3, AC021491.4, AC124067.4, and RHPN1-AS1) were identified to form the HCCSenLncSig predictive signature. A heatmap of the expression levels of these HCCSenLncSig lncRNAs in individual HCC patients is shown in Figure 3B. Cytoscape software was used to visualize the co-expression network of 36 pairs of senescence-related lncRNAs-mRNAs (Figure 3C, |R²| >0.3 and P<0.001). A Sankey diagram in Figure 3D shows that AC009283.1, AC021491.4, AC026401.3, and AC103760.1 were protective factors, whereas AC124067.4, AL117336.3, FOXD2-AS1, and RHPN1-AS1 were risk factors.

Figure 3 Identification of senescence-related lncRNAs related to HCC prognosis and lncRNA-mRNA co-expression network construction. (A) Forest plot showing 54 lncRNAs with hazard ratios, 95 percent confidence intervals, and P values for their associated HCC prognosis based on univariate Cox proportional hazards analysis. (B) Heatmap depicting the expression levels of 8 senescence-related lncRNAs identified by multivariate Cox regression analysis as being associated with HCC prognosis. (C) The lncRNA-mRNA co-expression network of the prognostic senescence-related lncRNA signature. Yellow squares represent prognostic lncRNAs, while green ellipses represent senescence-related mRNAs. Levels of expression of the 8 senescence-related lncRNAs were linked to the levels of 120 senescence mRNAs. (D) Sankey diagram showing the associations between prognostic senescence-related lncRNAs, mRNAs, and risk type. N, normal tissues; T, tumor tissues; lncRNAs, long-chain non-coding RNAs; HCC, hepatocellular carcinoma; mRNA, messenger RNA.

Correlation between HCCSenLncSig and the prognosis of HCC patients

The risk score of HCCSenLncSig was calculated as follows: risk score = (0.6488 × AL117336.3 expression) + (−0.3616 × AC103760.1 expression) + (0.5888 × FOXD2-AS1 expression) + (−0.7890 × AC009283.1 expression) + (−0.3235 × AC026401.3 expression) + (−1.1729 × AC021491.4 expression) + (0.3475 × AC124067.4 expression) + (0.9231 × RHPN1-AS1 expression). The formula was used to calculate the risk score for each patient, and patients were divided into 2 groups based on the median risk score, with 178 in the high-risk group and 165 in the low-risk group (Figure 4A). Kaplan-Meier analysis showed that the OS time of the high-risk group was significantly shorter than that of the low-risk group (Figure 4A). Figure 4B,4C show individual patient risk scores and survival statistics, with the number of deaths increasing as the risk score increases. Univariate and multivariate Cox regression analysis showed that the HCCSenLncSig risk score was an independent prognostic factor (Figure 4D,4E) and its AUC was 0.783, sensitivity of 0.600, and specificity of 0.896 at the cut-off value of 1.447, which showed that it was a better predictor of HCC prognosis than other clinicopathological variables (Figure 4F). The AUCs of the 1-, 3-, and 5-year ROC curves were 0.774, 0.713, and 0.793, respectively, indicating that the HCCSenLncSig had good prognostic performance (Figure 4G).

Figure 4 Prognostic value of the risk score determined by the HCCSenLncSig predictive model. (A) Kaplan-Meier curves for OS in the high- and low-risk groups stratified by the median of risk scores determined by the HCCSenLncSig. (B) Risk curve based on the risk score for each sample, where the yellow dot indicates a high-risk and blue dot indicates a low-risk. (C) Scatterplot based on the survival status of each sample. Yellow and blue dots indicate survival and death, respectively. (D) Forest plot for univariate Cox regression analysis. (E) Forest plot for multivariate Cox regression analysis. (F) ROC curve of the risk score and other clinicopathological variables. (G) Time-dependent ROC curves for 1-, 3-, and 5-year survival for the predictive signature. T, tumor; M, distant metastasis; N, lymph node metastasis; CI, confidence interval; AUC, area under the curve; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; OS, overall survival; ROC, receiver operating characteristic.

The expression levels of 8 lncRNAs from the HCCSenLncSig model and clinicopathological factors are shown in Figure 5A. We used PCA with whole-genome, senescence-related mRNAs, senescence-related lncRNAs, and the HCCSenLncSig to distinguish between high- and low-risk patients (Figure 5B-5E). The HCCSenLncSig (Figure 5E) model clearly distinguishes between low- and high-risk groups, showing the accuracy of the model.

Figure 5 Visualization of the expression levels of the 8 prognostic lncRNAs based on clinicopathological variable stratification, and PCA of the gene sets performed to classify patients into different risk groups. (A) Eight prognostic senescence-related lncRNAs and clinicopathological variables were distributed in a heatmap for high- and low-risk groups. PCA of low- and high-risk groups based on (B) whole-genome mRNA transcripts, (C) senescence-related mRNAs, (D) senescence-related lncRNAs, and (E) risk model including the 8 HCCSenLncSig senescence-related lncRNAs. Patients with high-risk scores are indicated in red, and those with low-risk scores are indicated in green. N, lymph node metastasis; M, distant metastasis; T, tumor; PC1, principal components 1; PC2, principal components 2; PC3, principal components 3; mRNA, messenger RNA; lncRNAs, long-chain non-coding RNAs; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; PCA, principal component analysis.

In addition, we assessed whether the HCCSenLncSig had prognostic value in subgroups with different clinicopathological parameters (Figure 6A-6P). There were significant correlations between the risk score and age (>65 and 65 years, Figure 6A,6B), male patients (Figure 6D), grade 1 to 3 (Figure 6E-6G), stage I and III (Figure 6I,6K), and T1, T3, and T4 (Figure 6M,6O,6P) when the correlations between the risk score and clinicopathological factors were assessed. In contrast, female patients (Figure 6C), grade 4 (Figure 6H), stage II and IV (Figure 6J,6L), and T2 (Figure 6N) were not associated with the risk score. The subgroup case numbers of M and N stage were too small to be evaluated. Taken together, the HCCSenLncSig risk score was proven to be an independent prognostic risk factor for patients with HCC.

Figure 6 Kaplan-Meier survival curves for high- and low-risk groups of patients sorted by clinicopathological variables. (A,B) Age; (C,D) sex; (E-H) grade; (I-L) overall stage; (M-P) T stage. T, tumor.

Construction of a predictive nomogram

We combined the HCCSenLncSig risk score with other clinicopathological factors to develop a nomogram to guide the clinical assessment of prognosis to estimate the 1-, 3-, and 5-year survival probability of HCC patients (Figure 7A). According to the 3 calibration diagrams, the nomogram estimated that the mortality rate was similar to the actual mortality rate (Figure 7B-7D).

Figure 7 Construction and verification of the nomogram. (A) Nomogram combining clinicopathological variables and risk score predicts HCC patient 1-, 3-, and 5-year survival probabilities. The calibration curves assess the consistency between the actual and predicted OS at (B) 1 year, (C) 3 years, and (D) 5 years. T, tumor; N, lymph node metastasis; M, distant metastasis; HCC, hepatocellular carcinoma; OS, overall survival.

Internal validation of the HCCSenLncSig

To test whether the HCCSenLncSig was universally applicable to HCC cancer OS predictive values, 343 TCGA-LIHC patients were randomly divided into 2 cohorts for internal validation (n=172 in the first internal cohort and n=171 in the second internal cohort). In both the first and second internal cohorts, the OS of patients in the high-risk group was lower than that of the low-risk group (Figure 8A,8B), which is consistent with the results generated from the entire TCGA-LIHC dataset. In addition, based on the time-dependent ROC curves, the AUCs for 1-, 3-, and 5-year survival of the first internal cohort were 0.809, 0.741, and 0.833, respectively (Figure 8C). In the second internal cohort, the AUCs for 1-, 3-, and 5-year survival were 0.743, 0.66, and 0.765, respectively (Figure 8D). These results show that the HCCSenLncSig has good predictive performance in both internal validation cohorts, indicating that the prediction model is robust.

Figure 8 HCCSenLncSig internal validation for overall survival using the entire TCGA dataset. (A) Kaplan-Meier survival curves for the first internal cohort. (B) Kaplan-Meier survival curves for the second internal cohort. (C) In the first internal cohort, the ROC curve and AUCs at 1-, 3-, and 5-year survival were calculated. (D) In the second internal cohort, the ROC curve and AUCs at 1-, 3-, and 5-year survival were calculated. AUC, area under the curve; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic.

Identification of biological pathways linked to the HCCSenLncSig

GSEA software was used for gene set functional annotation in HCCSenLncSig-identified high- and low-risk HCC patient groups. Specifically, the KEGG pathways of cell cycle [normalized enrichment score (NES) =2.18, NOM P<0.001, FDR q<0.002], pathogenic Escherichia coli infection (NES =2.12, NOM P<0.001, FDR q<0.001), homologous recombination (NES =2.11, NOM P<0.001, FDR q<0.001), Fc gamma R-mediated phagocytosis (NES =2.09, NOM P<0.001, FDR q<0.001), and oocyte meiosis (NES =2.02, NOM P<0.001, FDR q<0.004) were enriched in the high-risk group (Figure 9A). Drug metabolism-cytochrome p450 (NES =−2.2, NOM P<0.001, FDR q<0.001), tryptophan metabolism (NES =−2.09, NOM P<0.001, FDR q<0.001), fatty acid metabolism (NES =−2.09, NOM P<0.001, FDR q<0.001), tyrosine metabolism (NES =−2.04, NOM P<0.001, FDR q<0.002), and peroxisome (NES =−2.03, NOM P<0.001, FDR q<0.002) were enriched in the low-risk group (Figure 9A). On the other hand, the GO terms enriched in HCC patients with high-risk scores were positive regulation of spindle midzone (NES =2.25, NOM P<0.001, FDR q=0.005), chloride channel complex (NES =2.21, NOM P<0.001, FDR q=0.007), protein depolymerization (NES =2.19, NOM P<0.001, FDR q=0.009), regulation of microtubule cytoskeleton organization (NES =2.18, NOM P<0.001, FDR q=0.008), and microtubule (NES =2.18, NOM P<0.001, FDR q=0.007) (Figure 9B). In contrast, acylglycerol homeostasis (NES =−2.23, NOM P<0.001, FDR q<0.001), monocarboxylic acid catabolic process (NES =−2.2, NOM P<0.001, FDR q<0.001), organic acid catabolic process (NES =−2.19, NOM P<0.001, FDR q<0.001), lipid oxidation (NES =−2.18, NOM P<0.001, FDR q<0.001), and fatty acid catabolic process (NES =−2.18, NOM P<0.001, FDR q<0.001) were enriched in tumors of HCC patients with low-risk scores (Figure 9B).

Figure 9 HCCSenLncSig-based GSEA of different risk groups. (A) According to the GSEA results, KEGG genes were differentially enriched for senescence-related lncRNA expression. Five KEGG items, namely the cell cycle, pathogenic Escherichia coli infection, homologous recombination, Fc gamma R-mediated phagocytosis, and oocyte meiosis, were significantly differentially enriched in the high expression phenotype. Drug metabolism-cytochrome p450, tryptophan metabolism, fatty acid metabolism, tyrosine metabolism, and peroxisome were enriched in the low-risk group based on the NES, NOM P value, and FDR value. (B) Differential enrichment of genes in GO with senescence-related lncRNAs (5 GO items, namely positive regulation of spindle midzone, chloride channel complex, protein depolymerization, regulation of microtubule cytoskeleton organization, and microtubule) revealed significant differential enrichment in the high expression phenotype. The other 5 GO terms, namely acylglycerol homeostasis, monocarboxylic acid catabolic process, organic acid catabolic process, lipid oxidation, and fatty acid catabolic process, were found to be significantly enriched in the low expression phenotype based on the NES, NOM P value, and FDR value. KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; BP, biological process; CC, cellular components; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; GSEA, gene set enrichment analysis; lncRNA, long non-coding RNA; Fc, fold change; NES, normalized enrichment score; NOM, nominal; FDR, false discovery rate.

The relationship between HCCSenLncSig risk scores, driver mutations, and TMB

The somatic mutations in patients in low- and high-risk subgroups were assessed separately (Figure 10A,10B). TP53 (37% vs. 15%) had a higher rate of somatic mutation in the high-risk group, while CTNNB1 and TTN (33% vs. 17% and 27% vs. 19%, respectively) had a higher rate of somatic mutation in the low-risk group. In addition, for analysis of TMB levels, although there was no difference in TMB between the 2 groups (Figure 10C), the survival time of patients with higher TMB was significantly reduced (Figure 10D). In addition, the combination of high TMB in the high-risk group led to a worse prognosis (Figure 10E), demonstrating a significant synergistic effect between the 2 indicators.

Figure 10 The relationship between HCCSenLncSig risk scores, epigenetic mutation, and TMB. Waterfall plot showing the genetic mutations between high-risk (A) and low-risk (B) HCC patients. (C) Difference in TMB between patients from the low- and high-risk score subgroups. (D) Kaplan-Meier curves for the high and low TMB groups. (E) Kaplan-Meier curves for patients stratified by both TMB and risk scores. The P value represents the ANOVA test between the subgroups. TMB, tumor mutation burden; HCCSenLncSig, hepatocellular carcinoma senescence-related lncRNA predictive signature; HCC, hepatocellular carcinoma; ANOVA, analysis of variance.

The landscape of immune infiltration in various HCC risk subgroups

Immune cells in the TME have been proven to undergo senescence, and thereby modulate cancer hallmarks and consequent tumor phenotypes (33). Therefore, we investigated the differences in immune cells and immune functions between the 2 risk groups. The results showed that immune cell types such as activated dendritic cells (aDCs), dendritic cells (DCs), macrophages, T follicular helper cells (Tfh), T helper type 2 cells (Th2), and T regulatory cells (Treg) in the high-risk group were significantly different from those in the low-risk group (Figure 11A). Moreover, the scores of the immune functions, such as antigen-presenting cell (APC) co-inhibition, APC co-stimulation, C-C chemokine receptor (CCR), checkpoint, human leukocyte antigen (HLA), major histocompatibility complex (MHC) class I, parainflammation, T cell co-inhibition, T cell co-stimulation, and type II IFN response were also significantly different between the 2 groups (Figure 11B). To further investigate the relationship between risk score and immune cell infiltration, we used the CIBERSORT algorithm to estimate the proportions of 22 immune cell types in the HCC tumor immune microenvironment. As shown in Figure 11C, activated memory CD4 T cells (P=0.013), activated NK cells (P=0.009), monocytes (P<0.001), M0 macrophages (P<0.001), resting DCs (P=0.037), resting mast cells (P<0.001), and neutrophils (P=0.005) varied significantly between high- and low-risk score patients. The CIBERSORT findings partially confirmed the ssGSEA findings and, at the very least, revealed differences in immune infiltration between the high- and low-risk groups. Based on the theory that immunotherapy must rely on a pre-existing immune-hot microenvironment (34), these differences provide a potential therapeutic guide for immunotherapy.

Figure 11 Immune cell infiltration and immune-related functions in different risk groups. (A) The ssGSEA algorithm was used to compute the levels of infiltration of 16 immune cells in the high- and low-risk groups. (B) The relationship between risk score and 13 immune-related functions. (C) The Wilcoxon rank-sum test was used to explore whether there were differences between the 22 kinds of immune cells in different groups. *, P<0.05; **, P<0.01; ***, P<0.001; ns, non-significant. aDCs, activated dendritic cells; iDCs, immature dendritic cells; NK, natural killer; pDCs, plasmacytoid dendritic cells; Tfh, T follicular helper; Th1, T helper type 1; Th2, T helper type 2; TIL, tumor-infiltrating lymphocyte; Treg, T regulatory cell; APC, antigen-presenting cell; CCR, C-C chemokine receptor; HLA, human leukocyte antigen; MHC, major histocompatibility complex; IFN, interferon; ssGSEA, single-sample gene set enrichment analysis.

Potential relationship between the HCCSenLncSig and HCC immunotherapy, chemotherapy, and targeted therapy

There were 27 genes out of 47 evaluated immune checkpoints whose expression levels differed between the high- and low-risk groups (Figure 12A). Immunotherapy markers such as PDCD-1 and CTLA-4, which are now well established in clinical use, were found to be significantly higher in the high-risk group (Figure 12A), implying potential immunotherapeutic responses in high-risk patients. Moreover, when we used online software “TIDE” to predict the outcome of cancer patients treated with first-line anti-PD1 or anti-CTLA4 therapy, as shown in Figure 12B, we found no difference in TIDE scores between the 2 HCC subgroups. This finding could be explained by the fact that the TIDE database was created using melanoma data and only evaluated the efficacy of anti-CTLA4 and anti-PD-1 therapies without considering other ICIs, such as anti-PD-L1 (35).

Figure 12 Comparison of immune checkpoints, TIDE scores, and sensitivity of chemotherapy and targeted therapy drugs in high- and low-risk groups. (A) The expression of 27 immune checkpoint genes differed between the high- and low-risk groups. Red boxes represent high-risk patients, while blue boxes represent low-risk patients. (B) The online software “TIDE” predicted the TIDE score of the outcome of subgroups of HCC patients treated with anti-PD-1 or anti-CTLA4 therapies. The IC50 values for (C) 5-fluorouracil, (D) XL-184 (cabozantinib), (E) sunitinib, (F) gemcitabine, (G) paclitaxel, (H) imatinib, (I) bortezomib, and (J) erlotinib in the high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001; ns, non-significant; P<0.05 indicates statistical significance. TIDE, tumor immune dysfunction and exclusion module; HCC, hepatocellular carcinoma; PD-1, programmed cell death-1; IC50, half-maximal inhibitory concentration.

Finally, we investigated the relationship between the HCCSenLncSig risk score and the efficacy of chemotherapy and targeted therapy for HCC. Most of the drugs commonly used in preclinical and clinical systemic therapy for HCC were more sensitive in the low-risk group, such as 5-fluorouracil (Figure 12C), cabozantinib (XL-184, Figure 12D), sunitinib (Figure 12E), gemcitabine (Figure 12F), paclitaxel (Figure 12G), imatinib (Figure 12H), and bortezomib (Figure 12I), while erlotinib (Figure 12J) was more sensitive in the high-risk group. Taken together, our drug sensitivity study demonstrates that the HCCSenLncSig has the potential to play a role in the clinical development of personalized treatment strategies.

Discussion

Due to the control of hepatitis B and C, the incidence of liver cancer in China has gradually decreased, from 29.2 cases per 100,000 people in 1998 to 21.9 cases per 100,000 people in 2012 (36). However, the HCC cure rate remains unsatisfactory, which is due in part to a lack of biomarkers for treatment and prognosis. Traditionally used biomarkers, such as alpha-fetoprotein (AFP), can guide the diagnosis and monitoring of recurrence to some extent, but do not provide guidance for treatment (37). In addition, compared with a single clinical biomarker, combining multiple biomarkers into a single model can improve the prediction accuracy and help to develop an accurate personalized treatment plan (38). Senescent cells (of any cell type) have recently been identified as functionally significant cell types in the TME, including HCC (20). When HCC tissues expressed more senescence-related genes, patients had a lower OS and a shorter recurrence-free survival (39). When senescent cells (expressing p16INK4a+) were pharmacologically eliminated in aging mice, the incidence of spontaneous tumorigenesis and cancer-related death was reduced, and some age-related symptoms were delayed (40). The detection of senescent cells remains controversial at present. They are mostly identified by a combination of SASP expression, DNA damage, and β-galactosidase activity, which are neither specific nor universal (26,41-44). Recently, several lncRNAs have been found to be involved in the regulation of HCC cellular senescence (28,45-48), but the overall landscape of the regulatory network of senescent cells in HCC remains largely unknown. This situation prompted us to develop a HCCSenLncSig signature that could represent both senescence and HCC.

Based on the ROC curve, the HCCSenLncSig displayed moderate predictive performance for HCC OS. In addition, the nomogram incorporating the HCCSenLncSig and other clinical parameters is expected to improve clinical decision making and guide the development of treatment strategies. In the lncRNA model, both FOXD2-AS1 and RHPN1-AS1 have previously been identified as oncogenes in HCC, where FOXD2-AS1 aggravates HCC tumorigenesis by regulating the miR-206/MAP3K1 axis (49), and RHPN1-AS1 accelerates proliferation, migration, and invasion via regulating the miR-485-5p/BSG axis in HCC (50). Furthermore, AL117336.3 (51) and AC026401.3 (52) have been identified as HCC prognosis signatures. However, research on the prognostic value in cancer and contributions to senescence for AC103760.1, AC009283.1, AC021491.4, and AC124067.4 is lacking. As a result, more research is needed to further investigate the effects of these lncRNAs in HCC and senescence.

The highlight of this study is to emphasize the important relationship between cell senescence and the prognosis and microenvironment of HCC. At present, the role of the senescent microenvironment in tumors is often ignored in preclinical study, which are usually designed in young mice rather than old mice (53). This may help explain why many successful preclinical responses are not reproduced after entering real clinical trials. In addition, compared with other studies using the TCGA database to develop lncRNA prognostic models (13-18), our HCCSenLncSig model is comparable. This shows that the senescent cells in the TME are helpful in predicting the prognosis of HCC, and are at least as valuable as other biomarkers.

Clinical trials of immunotherapy-based treatment for HCC are ongoing. However, the ORR of immunotherapy is not satisfactory either as monotherapy or in combination. The ORR for anti-PD-1 was 14.3% (54,55), 17.6% for anti-CTLA-4 (6), and 20.1% for the combination of anti-CTLA-4 and anti-PD-L1 (5). These findings support the fact that the efficacy of immunotherapy depends on appropriate patient selection. Senescence, an important parameter, has recently been included in preclinical or clinical retrospective studies to screen individuals who may benefit from ICI therapy. In a melanoma mouse model, anti-CTLA-4, anti-PD-1, and anti-PD-L1 were all effective in young mice, but the latter was not effective in older mice (56). In clinical retrospective studies, patients over the age of 60 responded better to PD-1 than patients under the age of 60 (57). Similarly, the current study also demonstrated the relationship between the HCCSenLncSig model derived from senescent genes and the HCC immune microenvironment. According to HCCSenLncSig stratification, the expression of most immune checkpoints, activation of immune pathways, and infiltration of anti-tumor immune cells were higher in the high-risk group than in the low-risk group (Figures 11,12), suggesting that high-risk patients may benefit more from immunotherapy. Notably, existing studies have provided strong evidence that pre-existing anti-tumor immune responses in patients have been shown to be effective in the setting of immune checkpoint blockade therapy (34). On the other hand, high TMB has been linked to an increase in CD8+T cells, which recognize tumor neoantigens and thus produce a strong tumor killing effect (58). Although there was no difference in TMB in our model between the hierarchical TCGA-LIHC high-risk and the low-risk group (Figure 10C), individuals with a high-risk score and a high TMB had a poor prognosis (Figure 10E), indicating a significant synergistic effect between the 2 indicators. These findings emphasize the potential future direction of tumor immunotherapy by focusing on cellular senescence treatment.

Another aspect of the senescent impact on cancer phenotypes is the senescence of stromal cells. Despite the fact that most traditional drugs have been used to induce senescence and thus programmed cell death, stromal cells in the TME, despite having experienced a senescence phenotype, still mediate tumor survival via SASP (20,59). In a breast cancer model, senescent cells significantly increased cell proliferation, invasion, and metastasis by secreting IL-6, IL-8, IL-10 (60), CXCL11 (61), and VCAM1 (62). In such cases, senescent cellular targeting treatment strategies should be well considered during therapy (53). By using HCCSenLncSig-based risk scores, we can identify patients who are likely to benefit from current first-line HCC systemic drug therapy. Because the HCCSenLncSig model takes senescence into account, it may benefit clinical treatment options.

There are several limitations to the current study. First and foremost, additional external data should be considered to test whether the HCCSenLncSig model adequately fits the dataset. Second, some prognostic factors such as surgical data were not included in the nomograms owing to the lack of these data in TCGA. This may affect the accuracy of the model. Third, functional studies are required to better understand the molecular mechanism of senescence-related lncRNA effects.

In conclusion, we developed a HCCSenLncSig lncRNA signature that can be used to predict the prognosis of HCC. Importantly, the level of immune infiltration and the potential efficacy of tumor immunotherapy are related to the HCCSenLncSig.

Acknowledgments

The authors would like to thank The Cancer Genome Atlas (TCGA) for providing useful RNA-seq data with detailed clinical information (https://tcga-data.nci.nih.gov/tcga/). The authors also appreciate the academic support from the AME Hepatocellular Carcinoma Collaborative Group.

Funding: This work was supported by National Key Clinical Discipline, the Basic and Applied Basic Research Fund Project of Guangdong Province (No. 2021A1515410004 to TCZ).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://atm.amegroups.com/article/view/10.21037/atm-22-3348/rc

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://atm.amegroups.com/article/view/10.21037/atm-22-3348/coif). TBC received support from Astra Zeneca, Bayer, Eli Lilly, Genentech/Roche, Bristol-Myers Squibb, Merck KGaA, Merck Sharp & Dohme, and Astellas. TCZ reports funding support from the National Key Clinical Discipline, the Basic and Applied Basic Research Fund Project of Guangdong Province (No. 2021A1515410004). The other 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 (as revised in 2013).

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. Jemal A, Ward EM, Johnson CJ, et al. Annual Report to the Nation on the Status of Cancer, 1975-2014, Featuring Survival. J Natl Cancer Inst 2017;109:djx030. [Crossref] [PubMed]
2. El-Serag HB, Rudolph KL. Hepatocellular carcinoma: epidemiology and molecular carcinogenesis. Gastroenterology 2007;132:2557-76. [Crossref] [PubMed]
3. Zhang FJ, Yang JT, Tang LH, et al. Effect of X-ray irradiation on hepatocarcinoma cells and erythrocytes in salvaged blood. Sci Rep 2017;7:7995. [Crossref] [PubMed]
4. Giraud J, Chalopin D, Blanc JF, et al. Hepatocellular Carcinoma Immune Landscape and the Potential of Immunotherapies. Front Immunol 2021;12:655697. [Crossref] [PubMed]
5. Abou-Alfa GK, Chan SL, Kudo M, et al. Phase 3 randomized, open-label, multicenter study of tremelimumab (T) and durvalumab (D) as first-line therapy in patients (pts) with unresectable hepatocellular carcinoma (uHCC): HIMALAYA. J Clin Oncol 2022;40:abstr 379.
6. Uson PLS Junior, Liu AJ, Sonbol MB, et al. Immunotherapy and chimeric antigen receptor T-cell therapy in hepatocellular carcinoma. Chin Clin Oncol 2021;10:11. [Crossref] [PubMed]
7. Zhang G, Ha SA, Kim HK, et al. Combined analysis of AFP and HCCR-1 as an useful serological marker for small hepatocellular carcinoma: a prospective cohort study. Dis Markers 2012;32:265-71. [Crossref] [PubMed]
8. Baj J, Bryliński Ł, Woliński F, et al. Biomarkers and Genetic Markers of Hepatocellular Carcinoma and Cholangiocarcinoma-What Do We Already Know. Cancers (Basel) 2022;14:1493. [Crossref] [PubMed]
9. Shen S, Lin Y, Yuan X, et al. Biomarker MicroRNAs for Diagnosis, Prognosis and Treatment of Hepatocellular Carcinoma: A Functional Survey and Comparison. Sci Rep 2016;6:38311. [Crossref] [PubMed]
10. Zhu Y, Wang S, Xi X, et al. Integrative analysis of long extracellular RNAs reveals a detection panel of noncoding RNAs for liver cancer. Theranostics 2021;11:181-93. [Crossref] [PubMed]
11. Cui Y, Xu HF, Liu MY, et al. Mechanism of exosomal microRNA-224 in development of hepatocellular carcinoma and its diagnostic and prognostic value. World J Gastroenterol 2019;25:1890-8. [Crossref] [PubMed]
12. Zuo X, Chen Z, Gao W, et al. M6A-mediated upregulation of LINC00958 increases lipogenesis and acts as a nanotherapeutic target in hepatocellular carcinoma. J Hematol Oncol 2020;13:5. [Crossref] [PubMed]
13. Chen W, Hu MJ, Zhong XL, et al. Screening of a novel autophagy-related prognostic signature and therapeutic targets in hepatocellular carcinoma. J Gastrointest Oncol 2021;12:2985-98. [Crossref] [PubMed]
14. Zhao X, Bai Z, Li C, et al. Identification of a Novel Eight-lncRNA Prognostic Signature for HBV-HCC and Analysis of Their Functions Based on Coexpression and ceRNA Networks. Biomed Res Int 2020;2020:8765461. [Crossref] [PubMed]
15. Bai Z, Li H, Li C, et al. Integrated analysis identifies a long non-coding RNAs-messenger RNAs signature for prediction of prognosis in hepatitis B virus-hepatocellular carcinoma patients. Medicine (Baltimore) 2020;99:e21503. [Crossref] [PubMed]
16. Tang P, Qu W, Wang T, et al. Identifying a Hypoxia-Related Long Non-Coding RNAs Signature to Improve the Prediction of Prognosis and Immunotherapy Response in Hepatocellular Carcinoma. Front Genet 2021;12:785185. [Crossref] [PubMed]
17. Feng Y, Hu X, Ma K, et al. Genome-Wide Screening Identifies Prognostic Long Noncoding RNAs in Hepatocellular Carcinoma. Biomed Res Int 2021;2021:6640652. [Crossref] [PubMed]
18. Deng B, Yang M, Wang M, et al. Development and validation of 9-long Non-coding RNA signature to predicting survival in hepatocellular carcinoma. Medicine (Baltimore) 2020;99:e20422. [Crossref] [PubMed]
19. Kowald A, Passos JF, Kirkwood TBL. On the evolution of cellular senescence. Aging Cell 2020;19:e13270. [Crossref] [PubMed]
20. Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov 2022;12:31-46. [Crossref] [PubMed]
21. Wang B, Kohli J, Demaria M. Senescent Cells in Cancer Therapy: Friends or Foes? Trends Cancer 2020;6:838-57. [Crossref] [PubMed]
22. 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]
23. Kessenbrock K, Plaks V, Werb Z. Matrix metalloproteinases: regulators of the tumor microenvironment. Cell 2010;141:52-67. [Crossref] [PubMed]
24. Fisher DT, Appenheimer MM, Evans SS. The two faces of IL-6 in the tumor microenvironment. Semin Immunol 2014;26:38-47. [Crossref] [PubMed]
25. Huang Z, Zhou JK, Peng Y, et al. The role of long noncoding RNAs in hepatocellular carcinoma. Mol Cancer 2020;19:77. [Crossref] [PubMed]
26. Casella G, Munk R, Kim KM, et al. Transcriptome signature of cellular senescence. Nucleic Acids Res 2019;47:7294-305. [Crossref] [PubMed]
27. Li XL, Subramanian M, Jones MF, et al. Long Noncoding RNA PURPL Suppresses Basal p53 Levels and Promotes Tumorigenicity in Colorectal Cancer. Cell Rep 2017;20:2408-23. [Crossref] [PubMed]
28. Xiang X, Fu Y, Zhao K, et al. Cellular senescence in hepatocellular carcinoma induced by a long non-coding RNA-encoded peptide PINT87aa by blocking FOXM1-mediated PHB2. Theranostics 2021;11:4929-44. [Crossref] [PubMed]
29. Avelar RA, Ortega JG, Tacutu R, et al. A multidimensional systems biology analysis of cellular senescence in aging and disease. Genome Biol 2020;21:91. [Crossref] [PubMed]
30. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
31. Rooney MS, Shukla SA, Wu CJ, et al. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 2015;160:48-61. [Crossref] [PubMed]
32. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
33. Lian J, Yue Y, Yu W, et al. Immunosenescence: a key player in cancer development. J Hematol Oncol 2020;13:151. [Crossref] [PubMed]
34. Mlecnik B, Bindea G, Angell HK, et al. Integrative Analyses of Colorectal Cancer Show Immunoscore Is a Stronger Predictor of Patient Survival Than Microsatellite Instability. Immunity 2016;44:698-711. [Crossref] [PubMed]
35. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
36. Ferlay J, Colombet M, Bray F. Cancer incidence in five continents, CI5plus: IARC CancerBase No. 9. Lyon, France: International Agency for Research on Cancer, 2018.
37. Muscari F, Maulat C. Preoperative alpha-fetoprotein (AFP) in hepatocellular carcinoma (HCC): is this 50-year biomarker still up-to-date? Transl Gastroenterol Hepatol 2020;5:46. [Crossref] [PubMed]
38. Guo Y, Qu Z, Li D, et al. Identification of a prognostic ferroptosis-related lncRNA signature in the tumor microenvironment of lung adenocarcinoma. Cell Death Discov 2021;7:190. [Crossref] [PubMed]
39. Eggert T, Wolter K, Ji J, et al. Distinct Functions of Senescence-Associated Immune Responses in Liver Tumor Surveillance and Tumor Progression. Cancer Cell 2016;30:533-47. [Crossref] [PubMed]
40. Baker DJ, Childs BG, Durik M, et al. Naturally occurring p16(Ink4a)-positive cells shorten healthy lifespan. Nature 2016;530:184-9. [Crossref] [PubMed]
41. van Deursen JM. The role of senescent cells in ageing. Nature 2014;509:439-46. [Crossref] [PubMed]
42. Kastenhuber ER, Lowe SW. Putting p53 in Context. Cell 2017;170:1062-78. [Crossref] [PubMed]
43. Georgakilas AG, Martin OA, Bonner WM. p21: A Two-Faced Genome Guardian. Trends Mol Med 2017;23:310-9. [Crossref] [PubMed]
44. Lee BY, Han JA, Im JS, et al. Senescence-associated beta-galactosidase is lysosomal beta-galactosidase. Aging Cell 2006;5:187-95. [Crossref] [PubMed]
45. Jafari-Oliayi A, Asadi MH. SNHG6 is upregulated in primary breast cancers and promotes cell cycle progression in breast cancer-derived cell lines. Cell Oncol (Dordr) 2019;42:211-21. [Crossref] [PubMed]
46. Zhao L, Hu K, Cao J, et al. lncRNA miat functions as a ceRNA to upregulate sirt1 by sponging miR-22-3p in HCC cellular senescence. Aging (Albany NY) 2019;11:7098-122. [Crossref] [PubMed]
47. Jia Y, Jin H, Gao L, et al. A novel lncRNA PLK4 up-regulated by talazoparib represses hepatocellular carcinoma progression by promoting YAP-mediated cell senescence. J Cell Mol Med 2020;24:5304-16. [Crossref] [PubMed]
48. Si L, Chen J, Yang S, et al. lncRNA HEIH accelerates cell proliferation and inhibits cell senescence by targeting miR-3619-5p/CTTNBP2 axis in ovarian cancer. Menopause 2020;27:1302-14. [Crossref] [PubMed]
49. Hu W, Feng H, Xu X, et al. Long noncoding RNA FOXD2-AS1 aggravates hepatocellular carcinoma tumorigenesis by regulating the miR-206/MAP3K1 axis. Cancer Med 2020;9:5620-31. [Crossref] [PubMed]
50. Zhang W, Han L, Xing P, et al. LncRNA RHPN1-AS1 accelerates proliferation, migration, and invasion via regulating miR-485-5p/BSG axis in hepatocellular carcinoma. Naunyn Schmiedebergs Arch Pharmacol 2020;393:2543-51. [Crossref] [PubMed]
51. Jia Y, Chen Y, Liu J. Prognosis-Predictive Signature and Nomogram Based on Autophagy-Related Long Non-coding RNAs for Hepatocellular Carcinoma. Front Genet 2020;11:608668. [Crossref] [PubMed]
52. Cheng M, Zhang J, Cao PB, et al. Prognostic and predictive value of the hypoxia-associated long non-coding RNA signature in hepatocellular carcinoma. Yi Chuan 2022;44:153-67. [PubMed]
53. Fane M, Weeraratna AT. How the ageing microenvironment influences tumour progression. Nat Rev Cancer 2020;20:89-106. [Crossref] [PubMed]
54. El-Khoueiry AB, Sangro B, Yau T, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet 2017;389:2492-502. [Crossref] [PubMed]
55. El-Khoueiry AB, Melero I, Yau TC, et al. Impact of antitumor activity on survival outcomes, and nonconventional benefit, with nivolumab (NIVO) in patients with advanced hepatocellular carcinoma (aHCC): Subanalyses of CheckMate-040. J Clin Oncol 2018;36:abstr 475.
56. Padrón Á, Hurez V, Gupta HB, et al. Age effects of distinct immune checkpoint blockade treatments in a mouse melanoma model. Exp Gerontol 2018;105:146-54. [Crossref] [PubMed]
57. Kugel CH 3rd, Douglass SM, Webster MR, et al. Age Correlates with Response to Anti-PD1, Reflecting Age-Related Differences in Intratumoral Effector and Regulatory T-Cell Populations. Clin Cancer Res 2018;24:5347-56. [Crossref] [PubMed]
58. Chan TA, Yarchoan M, Jaffee E, et al. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol 2019;30:44-56. [Crossref] [PubMed]
59. 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]
60. Birch J, Gil J. Senescence and the SASP: many therapeutic avenues. Genes Dev 2020;34:1565-76. [Crossref] [PubMed]
61. Hwang HJ, Lee YR, Kang D, et al. Endothelial cells under therapy-induced senescence secrete CXCL11, which increases aggressiveness of breast cancer cells. Cancer Lett 2020;490:100-10. [Crossref] [PubMed]
62. Wang D, Xiao F, Feng Z, et al. Sunitinib facilitates metastatic breast cancer spreading by inducing endothelial cell senescence. Breast Cancer Res 2020;22:103. [Crossref] [PubMed]

(English Language Editor: C. Betlazar-Maseh)