Transcriptomic data and clinical information of esophageal cancer (EC) derived from the TCGA-ESCA cohort as a training set (https://portal.gdc.cancer.gov/), involving 162 EC samples and 11 normal samples. Clinical information not available or ambiguous was removed. Independent cohort GSE53625 as validation set available from GEO database. Cellular senescence-related genes (CSRGs) were selected by the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/) and Genecards (https://www.genecards.org/) tools (Supplementary Table S1). The procedure detailed in this study is shown in Figure 1.

Statistical analyses based on the TCGA database were performed with R. The differentially expressed genes (DEGs) in tumor and normal tissues of TCGA-ESCA cohort were screened by differential analysis. Combined with CS-related genes, CS-related DEGs in EC were initially screened by Venn analysis. The WGCNA weighting analysis of the distribution of correlation modules of these genes was performed, and CS-related prognostic genes were further obtained by univariate COX regression analysis. Finally, CS-related prognostic hub genes were identified by LASSO regression.

Based on the coefficients of CS-related prognostic hub genes gained from Lasso regression analysis, the CS-related risk scores (CSRS) were constructed as follows.CS_relatedriskscores(CSRS)=∑i=1nexpressiongene_i×lasso_coeffieicentgene_i

Independent prognostic factors were screened by univariate and multivariate COX regression analysis. These factors and CSRS were combined to construct a nomogram model for predicting survival in patients with EC. A preliminary assessment was performed with a calibration correction curve.

Data from the GSE53625 dataset was taken to validate the reliability of the model. The effectiveness of the nomogram model was demonstrated by the decision curve analysis (DCA) curve, Kaplan-Meier (KM) curve and receiver operating characteristic (ROC) curve.

The Wilcoxon signed-rank sum test was used to compare the differences in clinical characteristics of patients in high- and low- CSRS groups. The prognostic value of CSRS for patients of different age groups, pathological staging, and pathological stages was performed by Kaplan-Meier.

In the present study, three scoring systems for immuno-infiltration analysis were performed, namely ssGSEA analysis (14), ESTIMATE scores (15) and TIDE scores (16). Levels of infiltration of different immune cells in tumors were quantified by the ssGSEA algorithm through the GSVA package (17). The purity of tumor immune infiltration and abundance of stromal cells were calculated by ESTIMATE algorithm through the estimate package. The dysfunction score and exclusion scores from the TIDE scoring system were applied to predict the efficacy of immunotherapy in different CTL-related subgroups of patients.

GO analysis and KEGG analysis for probing the potential biological functions of gene networks in different modules of the WGCNA with the clusterProfiler package and org.HS.eg.db package (18). The biological mechanisms leading to differences in high and low CSRS groups were explored via gene set enrichment analysis (GSEA) by the clusterProfiler package (17, 18).

A total of 1,153 CS-related genes were derived by MSigDB and Genecards tools (Supplementary Table S1), of which 241 genes (Figure 2A) were differentially expressed between EC and normal tissues (Figure 2B, |log₂FC|>1, p < 0.05).

WGCNA analysis of TCGA-ESCA transcriptome data was performed to search for highly related gene modules. Based on the relationship between the soft threshold with the scale-free fit and the mean connectivity, a suitable soft threshold β was finally determined as 12 (Figure 2C). The network classified the CS-related DEGs into three different modules, blue, brown and turquoise (Figures 2D,E), by using a dynamic tree cutting and clustering algorithm. The correlation between modules was presented by a heat map, which showed that the turquoise module was highly genetically correlated with the brown module.

GO/KEGG analysis was performed to probe the biological functions associated with each module gene. The genes of the blue module were mainly enriched in cellular senescence and aging (Figure 2F). The genes of the brown and turquoise modules (Figure 2G) might play a role in biological processes such as cellular senescence, as well as, apoptosis-related signaling pathways. The detailed GO/KEGG annotations are presented in Table 1.

Univariate COX regression analysis of the modular genes identified seven genes that were strongly associated with overall survival (OS), namely SLC30A10, IGFBP1, H3C1, FBXO5, SOX5, CDHR4 and MT1E (Figure 2H). The above genes were identified as CS-related prognostic genes.

The regression coefficients (Table 2) of the above 7 CS-related prognostic genes were calculated by the Lasso algorithm (Figures 3A,B) using OS as an outcome indicator, with the CSRS=0.2901×H3C1+0.2158∗IGFBP1−0.7121∗CDHR4−0.1390∗MT1E−0.1184∗SOX5The prognostic DCA chart (Figure 3C) confirmed the utility of the CSRS scoring system in predicting survival outcomes in patients with EC.

We performed a COX regression analysis of the TCGA-ESCA cohort to uncover factors affecting the prognosis of esophageal patients. In the independent cohort GSE53625, EC patients were divided into high-risk and low-risk groups based on the median CSRS in TCGA-ESCA as the cutoff value for further analysis to verify the generalizability of the CSRS score. The results of the univariate COX analysis in the TCGA cohort (Figure 3D) suggested that N stage, M stage, pathological stage and CSRS (HR = 2.903, 95%CI = 1.497–5.629, p = 0.002) were risk factors affecting the prognosis of esophageal cancer, which was similarly validated in the GSE53625 cohort (Figure 3E, risks score group: HR = 1.742, 95%CI = 1.129–2.686, p = 0.012). Further multivariate COX analysis at TCGA-ESCA (Figure 3F) and GSE53625 (Figure 3G) indicated the reliability of the prediction of prognosis in patients with EC by CSRS. CSRS can accurately distinguish esophageal cancer patients with different survival times, which means that a higher CSRS represents a worse prognosis as reflected by the results of the KM analysis (Figures 3H,I).

Integrating the above analysis, we constructed a nomogram model to predict the 1-,2- and 3-year survival of EC patients based on N stage, M stage, pathological stage and CSRS (Figure 3J). The fit is around the diagonal and the C-index value is 0.744, indicating good consistency of the model (Figure 3K). In addition, we evaluated the efficacy of the nomogram model. The DCA curve (Supplementary Figure S1A) results showed that the prediction of survival outcome in patients with EC using the CSRS was superior to that using American Joint Committee on Cancer (AJCC) staging. The benefit of prediction using our constructed nomogram model was greater than that of CSRS and AJCC. The KM curve (Supplementary Figure S1B) results showed that patients with high nomogram scores had a worse prognosis (HR = 5.35, 95% CI = 2.61–10.96, p < 0.001). The accuracy of the nomogram model in predicting the 1-(AUC = 0.781),3-(AUC = 0.754) and 5 (AUC = 0.946) years’ prognosis of patients with EC was also assessed by time-dependent ROC analysis (Supplementary Figure S1C).

We observed no significant difference in the distribution of CSRS among EC groups by gender (Figure 4A), age (Figure 4B), and BMI (Figure 4C). However, in terms of pathological type (Figure 4D), CSRS was higher in patients with esophageal adenocarcinoma (EAC) than those with esophageal squamous cell cancer (ESCC).

For this reason, we investigated the prognostic value of CSRS in different subgroups of patients with EC (Figure 4E–J). CSRS accurately determined prognosis in patients with either EAC (HR = 3.12, 95%CI = 1.59–6.13, p = 0.001) or ESCC (HR = 5.68, 95%CI = 2.10–15.39, p = 0.001), as well as in patients with EC aged more than 65 years (HR = 4.35, 95%CI = 1.67–11.31, p = 0.003) or T3 stage (HR = 4.24, 95%CI = 1.77–10.14, p = 0.001) or N1&N2&N3 stage (HR = 3.34, 95%CI = 1.69–6.98, p = 0.001) or M0 stage (HR = 2.20, 95%CI = 1.16–4.16, p = 0.016) or pathological stage III & IV (HR = 2.50, 95%CI = 1.13–5.52, p = 0.023). However, for patients aged less than 65 years, T1 & T2 stages, N0 stages, and pathological stages I & II, CSRS scores were not good predictors of prognostic outcome.

We adopted three scoring systems to analyze tumor immune infiltration in EC patients with different CSRS groups, namely ssGSEA analysis (Figures 5A–C), ESTIMATE score (Figures 5D–F) and TIDE score (Figures 5G,H). EC patients in the high CSRS group were infiltrated by fewer Tc and Tgd cells, while there was a positive correlation with the infiltration of neutrophil cells. Stromal scores (r = −0.178, p = 0.024) and ESTIMATE scores (r = −0.189, p = 0.016) were observed to be negatively correlated with CSRS, whereas not immune scores. There were differences in all three scores between the high- and low-CSRS groups of esophageal cancer. The TIDE scoring system is commonly used to evaluate the efficacy of immunotherapy in oncology patients, including the exclusion score and dysfunction score of T cells. CSRS was negatively correlated with dysfunction scores (r = −0.214, p = 0.011), and no significant correlation was observed with exclusion scores. We subsequently compared the expression of immune checkpoint-related genes in different CSRS groups (Figure 5I). High expression of PDL1, LAG3 and TIGIT were observed in low-CSRS group (p < 0.05) than high-CSRS group.

In order to explore the biological mechanisms leading to differences between high- and low- CSRS groups, GSEA analysis was performed. The results showed that the high-CSRS group was positively enriched in acetylation- (Figure 6A) and methylation-related (Figure 6B) pathways, and negatively enriched in immunomodulatory (Figure 6C) and GPCR-related pathways (Figure 6D).