The gene expression profiles and clinical data of patients with ECa were extracted from the TCGA public database (https://portal.gdc.cancer.gov) and used as the training cohort. A total of 194 patients (12 normal esophageal tissue samples and 182 ECa samples) were included in the TCGA database. Independent probe matrix file (GSE53624) and platform file (GPL18109) containing 179 samples were derived from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) and served as validation sets. According to the probe sequence of the platform file, the gene names were obtained by chip re-annotation technology, to obtain the corresponding relationship between the probe matrix and the gene names. Finally, we selected a list of 279 CSRGs ( Supplementary Table S1 ) from the CellAge database (https://genomics.senescence.info/cells/). The flow diagram of this study is depicted in Figure 1 .

The expression differences of CSRGs between 12 normal esophageal tissues and 182 ECa tissues were analyzed by the “limma” R package. The data were analyzed strictly according to the screening criteria of false discovery rate (FDR)<0.05 and |log₂ (fold change, FC)|>0.585. Next, differentially expressed CSRGs were visualized by plotting heat map and volcano map with the “pheatmap” package of the R software system.

We first combined the expression data in TCGA-ECa with the survival data and then analyzed and obtained the expression levels of differentially expressed CSRGs in ECa samples. Then, the CSRGs related to the prognosis of ECa were obtained by univariate Cox regression analysis, and the screening criteria was P<0.05. Finally, the Least Absolute Shrinkage and Selection Operator (lasso) regression analysis with the “glmnet” package in R was executed to pick up the preliminary hub CSRGs. The calculation formula of CSRGs-related risk scores was as follows: risk scores =Σ_i(Lasso_Coef_i*GeneExp_i). “Lasso_Coef”, Lasso regression coefficient; “GeneExp”, amount of gene expression. According to the median risk score derived from CSRGs, patients with ECa were divided into high- and low-risk groups.

The effectiveness of the risk signature was demonstrated by principal component analysis (PCA), Kaplan–Meier survival curves, and receiver operating characteristic (ROC) curves. In order to verify the independence of the risk signature in prediction, univariate and multivariate Cox analyses on clinical variables and risk scores were generated using the “survival” package of R. The selected clinical variables are mainly indicators that have important prognostic significance for patients with ECa, such as age, gender, clinicopathological grade, and TNM stage. The “ggpubr” package of R was used to investigate whether there were differences in clinical variables and immune subtypes among different risk groups.

Total RNA was extracted from 12 paired human ECa tissues and adjacent non-tumorous esophageal tissues using TRIzol Reagent (Abcam, China). The reverse transcription was conducted with TransScrip All-in-One SuperMix for qPCR reagent Kit (TaKaRa, Japan). Real-time fluorescent PCR was performed by SYBR Green assay. The experiment adopted a 20μL reaction system, Tip Green qPCR SuperMix (TaKaRa, Japan), cDNA template, upstream and downstream primers were successively added into 8 reaction tubes, and 3 repeated experiments were performed for each sample. We used glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as the internal reference, and the data were analyzed using the 2^–ΔΔCt approach. The sequences of primers were listed in Supplementary Table S2 .

In order to predict the survival rate of patients with ECa at 1, 2, and 3 years, clinicopathological factors and risk score integrated nomogram was generated using the”regplot” and “rms” packages in R. Subsequently, by generating a calibration curve, we evaluated the preliminary consistency between the survival rate of ECa predicted by the nomogram and the actual survival rate. The ROC curve was utilized to explore the accuracy of the nomogram and clinical features in predicting the survival rate of patients with ECa. At last, univariate and multivariate Cox analyses were utilized to explore the potential of the nomogram to independently predict prognosis in ECa.

R software combined with the CIBERSORT algorithm was used to explore the differences of 22 kinds of human immune cell subpopulations between the two risk groups in TCGA-ECa. Immediately after, we explored the differences in pathway enrichment between the two risk groups by gene set variation analysis (GSVA). The TIDE (Tumor Immune Dysfunction and Exclusion) algorithm was used to evaluate the weight of immunological rejection in two risk groups via an online website (http://tide.dfci.harvard.edu/).

Firstly, the “limma” package of R software was used to screen the differentially expressed genes (DEGs) between different groups, and the screening criteria were FDR<0.05 and |log₂FC|>1. Subsequently, we performed GO and KEGG analysis on the two risk groups to explore the differences in their potential biological functions and pathways.

Using the curated gene set (kegg.v7.4.symbols.gmt), broad GSEA v.4.2.3 was applied to detect high- and low-risk group correlation pathways with the criteria: NOM P<0.05 and |NES|>1 (22).

The STRING online database (https://string-db.org/) was first applied to obtain the PPI (protein-protein interaction) information (interaction score >0.70) of the DEGs between different groups. Next, we visualized the PPI network using the Cytoscape software (version 3.9.1). In addition, the plugin of cytoHubba in Cytoscape was utilized to screen the hub genes with the most complex connections in the PPI network. Finally, the clinical significance of hub genes in ECa was further explored.

Data processing, statistical analysis, and diagram formation were all conducted in R software (version R-4.2.2), GraphPad Prism (version 9.0), and Cytoscape software (version 3.9.1). The Kaplan–Meier curve plotted by the “survminer” package of R was used to compare differential survival probability. Univariate and multivariate Cox regression analyses of independent prognostic factors were performed using the “forestplot” package of R. ROC curve plotted by the “timeROC” package of R was used to assess the predictive efficacy of the CSRGs prognostic signature and nomogram. Results with two-sided P<0.05 were deemed statistically significant.

In TCGA data, among 279 CSRGs, 119 were differentially expressed between ECa samples and normal esophageal samples (FDR <0.05, |log₂FC| > 0.585), including 15 downregulated genes and 104 upregulated genes in the ECa tissue samples. These above differences are displayed in the heat map and volcano map ( Figures 2A, B ).

Based on univariate Cox regression analysis of the TCGA-ECa cohort, 7 CSRGs that significantly affect the prognosis of ECa were disclosed, such as TERT ( Figure 3A ). Subsequently, in order to develop a CSRGs signature for survival prediction of ECa, the 6 OS (overall survival)-associated CSRGs were analyzed using the LASSO analysis. Finally, a total of 6 genes (DEK, RUNX1, SMARCA4, SREBF1, TERT, TOP1) were built ( Figures 3B, C ). The risk scores for all patients were calculated as follows: risk score = (0.486102936551265×DEK level) + (-0.466577382724787×RUNX1 level) + (-0.379803436002195×SMARCA4 level) + (-0.216260668499574×SREBF1 level) + (0.326358456240447×TERT level) + (0.126266167208409×TOP1 level) ( Supplementary Table S3 ). The results of our further analysis of the expression levels of signature genes in ECa and normal esophageal tissues are as follows. Compared with the normal esophageal tissues group, DEK, RUNX1, SMARCA4, SREBF1, TERT, and TOP1 all showed elevated expression abundance in the ECa tissues group ( Figures 3D-I ).

Real-time fluorescence quantitative PCR (RT-qPCR) was performed to verify the mRNA level of the six signature genes in 12 paired ECa tissues and adjacent normal esophageal tissues. The results indicated that the expression of DEK, RUNX1, SMARCA4, SREBF1, TERT, and TOP1 were all up-regulated in ECa tissues than that in adjacent normal tissues ( Figures 4A-F ). These results were consistent with the expression tendencies of the previous signature genes in esophageal cancer and normal esophageal tissues.

It was found by PCA that our constructed CSRGs signature can accurately divide ECa samples into high- and low-risk groups ( Figures 5A, B ). Similar to the results obtained from the training cohort, patients in the high-risk group of the testing cohort were more likely to encounter a worse prognosis ( Figure 5C ). The results of univariate Cox regression analysis combined with multivariate Cox regression analysis in the TCGA-ECa cohort and GEO-ECa cohort suggested that the risk score based on CSRGs remained an independent risk factor affecting the prognosis of ECa ( Figures 5D, E ). In predicting survival, the risk score had a larger area under the ROC curve than other clinical features in the TCGA-ECa cohort and GEO-ECa cohort, suggesting that the risk score could serve as a more accurate prognostic factor ( Figure 5F ). Afterward, ROC analysis was employed to assess the risk signature in OS, with AUC values of 0.779, 0.720, and 0.761 at the 1-, 2-, and 3-year in TCGA-ECa cohort and 0.705, 0.721, and 0.745 at the 1-, 2-, and 3-year in the GEO-ECa cohort, respectively ( Figure 5G ).

Based on the results of risk grouping, we are interested in whether there are differences in clinical features among different risk groups and conducting further analysis. We found a significantly increased risk in patients older than 65 years and in patients with grade G3 ( Figures 6A, B ). However, other clinical features did not show statistical significance between the different risk score groups ( Supplementary Figure S1 ). Thorsson et al. performed an extensive immunogenomic analysis of over 10,000 tumors comprising 33 diverse cancer types utilizing data compiled by TCGA (23). Across cancer types, they identified six immune subtypes by RNA sequencing: C1-Wound Healing, C2-IFN-γ Dominant, C3-Inflammatory, C4-Lymphocyte Depleted, C5-Immunologically Quiet, and C6-TGF-β Dominant (23). Next, we also found that the risk score between the different immune subtypes was also not statistically different ( Figure 6C ). Finally, we investigated the prognostic value of CSRGs risk score signature in different subgroups of patients with ECa ( Figures 6D–H ). CSRGs risk score can accurately determine prognosis in ECa patients with either aged less than 65 years (P<0.001) or M0 stage (P<0.001), as well as in patients with ECa in T1&T2 stages (P=0.007) or T3&T4 stages (P=0.003) or N0 stage (P=0.015) or N1&N2&N3 stages (P=0.006) or pathological stage I & II (P=0.003) or pathological stage III & IV (P=0.004). However, CSRGs risk score signature was not a good predictor of prognostic outcome for ECa patients aged more than 65 years (P=0.170).

In order to broaden the application of the CSRGs in ECa patients, we developed a quantitative nomogram to compute OS at 1, 2, and 3 years in the TCGA-ECa cohort and GEO-ECa cohort ( Figure 7A ). In the TCGA-ECa cohort, when the prognostic parameters’ point is 255, the predicted OS of patients with ECa is 0.544 at 1 year, 0.186 at 2 years, and 0.054 at 3 years. The OS prediction lines of the nomogram are close to the 45° standard curve of the calibration analysis in the TCGA-ECa cohort and GEO-ECa cohort, demonstrating that the established clinical nomogram performs excellent ( Figure 7B ). It confirmed the high predictive efficiency of the nomogram in OS of patients with ECa. We also found that the AUC value of the nomogram in the ROC curve reached 0.808, suggesting that the nomogram was a better predictor of survival than other prognostic indicators of ECa patients in the TCGA cohort ( Figure 7C ). In addition, the results of univariate Cox regression analysis combined with multivariate Cox regression analysis ( Figures 7D, E ) suggested that the nomogram (HR=1.164, 1.097-1.235, P<0.001) based on risk score remained an independent risk factor affecting the prognosis of ECa in TCGA cohort.

By performing GSVA (Gene Set Variation Analysis), we evaluate the relative expression difference of the pathways between two risk groups. Many differentially expressed pathways were enriched by GSVA analysis and finally visualized by heatmap ( Figure 8A ). Compared with the low-risk group, the expression of pathways associated with metabolism and transport of cellular processes were remarkably activated in the high-risk group, whereas the expression of tumor and genetic information processing associated pathways were significantly lower. Analysis of immune cell infiltration revealed higher naïve B cells and regulatory T cells (Tregs) in the high-risk group, however, activation of M0 macrophages and activated memory CD4+ T cells was higher in the low-risk group ( Figure 8B ). In addition, compared with the high-risk group, immune function analysis displayed that APC_co_inhibition and type_II_IFN_response were more activated in the low-risk group ( Figure 8C ). Finally, the TIDE algorithm showed that the high-risk group has a greater risk of immune escape during immunotherapy, which also implied that the low-risk group of ECa patients may benefit from immunotherapy ( Figure 8D ).

We screened 747 DEGs between the two risk groups. GO and KEGG analysis further elucidated differences in biological functions and pathways between the two risk groups. The biological process (BP) modules of GO analysis are mainly focused on epidermis development, digestion, epidermal cell differentiation, digestive system process, and so on ( Figures 9A, B ). However, the human papillomavirus infection, protein digestion and absorption, metabolism of xenobiotics by cytochrome P450, and tight junction pathways were observably enriched in the KEGG analysis ( Figures 9C, D ). To explore the different biological functions of 6 cellular senescence-related signature genes in two risk groups, the GSEA (gene set enrichment analysis) analysis was used to identify the top five pathways. GSEA analysis showed that complement and coagulation cascades, glycerolipid metabolism, maturity-onset diabetes of the young, PPAR signaling pathway, and tryptophan metabolism were enriched in high-risk groups ( Figure 9E ). However, dilated cardiomyopathy, ECM receptor interaction, hedgehog signaling pathway, pathways in cancer, and regulation of actin cytoskeleton were enriched in low-risk groups ( Figure 9F ).

In order to screen the differential hub genes that are involved in two risk score groups from the interaction level, the expression profiles of DEGs were visualized by the STRING database and Cytoscape software. Proteins encoded by 293 differential genes were used to construct the PPI network, which included the interaction relationship of 177 up-regulated genes and 116 down-regulated genes ( Figure 10A ). Finally, the 10 hub genes (APOA1, MUC5AC, GC, APOA4, AMBP, FABP1, APOA2, SOX2, MUC8, MUC17) with the highest interaction degrees were identified by Cytoscape ( Figure 10B ). By further analysis, 4 hub genes (APOA4, AMBP, FABP1, and APOA2) with survival differences were identified ( Figure 10C ; Supplementary Figure S2 ). The results demonstrated that patients with ECa with high expression of the above 4 genes had a lower probability of survival.

We further analyzed the relationship between the expression level of these 4 genes and the clinical features of ECa. In ECa patients with grade G2 and G3, the expression level of APOA4 was significantly higher than that in grade G1 patients ( Supplementary Figure S3 ). With the increase of G grade level in patients with ECa, the expression of AMBP gradually increased with statistical significance ( Figure 10D ). In ECa patients over 65 years old, the FABP1 expression level was significantly higher than that below 65 years old ( Supplementary Figure S4 ). Finally, the relationship between the expression levels of these 4 genes and different immune cells was further explored. The AMBP and FABP1 high expression groups had higher resting CD4 memory T cells infiltration, suggesting that those groups of patients are more suitable for immunotherapy ( Figure 10E ; Supplementary Figure S5 ).