According to the 6th edition of the American Joint Committee on Cancer TNM staging system, the Guangzhou cohort included patients with ESCC staged IIb-III who underwent NCRT prior to surgery from September 2007 to March 2012. The RNA-seq data for these patients are accessible in the Gene Expression Omnibus (GEO) database under accession number GSE45670. We monitored these patients until July 2023, selecting those who survived over 3 months after treatment to assess their survival outcomes. Another cohort, the Beijing cohort (GSE53624), comprised tumor and adjacent normal tissues from 119 ESCC patients. The Series Matrix Files for GSE45670 and GSE53624 were based on annotation platforms GPL570 and GPL18109, respectively.

We retrospectively included the real-world patients with stage II-III ESCC from Sun Yat-sen University Cancer Center. Tissue samples and pathological results of these patients were obtained by endoscopic biopsy prior to treatment. After pathological diagnosis, they received NCRT combined with esophagectomy. This study was approved by the Ethics Committee of the Sun Yat-sen University Cancer Center and conducted in accordance with the local legislation and institutional requirements.

The “limma” package (Ritchie et al., 2015) was utilized to identify DEGs between pCR and npCR patients, while univariate Cox regression analysis determined genes associated with survival. We then identified genes that were highly expressed in pCR patients and suggested a favorable prognosis, and those underexpressed in pCR patients and indicated a poorer prognosis. The random survival forest algorithm executed using the “randomForestSRC” package, selected genes based on their prognostic importance. Genes with a relative importance exceeding 0.7 were classified as key genes related to chemoradiotherapy sensitivity.

Using the gene expression of each patient, the CIBERSORT deconvolution algorithm estimated the relative proportions of 22 immune cell types, including B cell subsets, T cell subsets, NK cells, and macrophages (Newman et al., 2015). We analyzed and compared the immune cell fractions in tumor and adjacent non-tumor tissues of ESCC patients employing the “CIBERSORT” and “ggpubr” R packages. Correlations between chemoradiotherapy sensitivity-related genes and immune cell fractions were examined using the “corrplot” R package.

To explore the associations between chemoradiotherapy sensitivity-related genes and immunoregulatory genes, we sourced gene encodings for immunomodulators from the TISIDB database (Ru et al., 2019), including 24 immunoinhibitors and 46 immunostimulators. After intersecting these with the ESCC gene expression profile, 60 genes encoding immunomodulators were retained. Pearson correlation analysis was conducted on the expression of immunoregulatory genes and chemoradiotherapy sensitivity-related genes. P < 0.05 was deemed statistically significant.

The Genomics of Drug Sensitivity in Cancer (GDSC) database records the responsiveness of cancer cells to drugs and molecular markers associated with drug response (Yang et al., 2013). Utilizing this pharmacogenomics database, the “oncoPredict” R package assessed the half-maximal inhibitory concentration (IC50) of chemotherapeutics for ESCC patients. We evaluated the impact of chemoradiotherapy sensitivity-related gene expression on chemotherapy sensitivity in these patients.

Employing the GEPIA database (http://gepia.cancer-pku.cn), we analyzed the differential expression of key chemoradiotherapy sensitivity-related genes in tumor samples and paired normal tissues across 32 cancer types, including ESCC.

Cox regression analyses were performed on clinical factors to identify independent prognostic indicators. We constructed a nomogram incorporating these prognostic factors and chemoradiotherapy sensitivity-related genes using the “rms” R package, estimating the 1-, 3-, and 5-year survival probabilities for ESCC patients. The nomogram’s consistency and accuracy were validated with a calibration curve.

Using the “fgsea” and “enrichplot” R packages and the annotated gene set ‘h.all.v2023.2.Hs.symbols.gmt’ from the Molecular Signatures Database, we performed fast gene set enrichment analysis (fGSEA) to identify potential pathways and biological functions differing among chemoradiotherapy sensitivity-related gene expression groups. Based on the Net enrichment score (NES) and P-value, we identified significantly enriched hallmark pathways and explored the mechanisms through which key genes influence these pathways.

Utilizing data from 452,264 individuals in the United Kingdom Biobank and documented in the Gene Atlas database (http://geneatlas.roslin.ed.ac.uk/), which links 778 traits to 30 million variants, we pinpointed pathogenic regions associated with key genes by analyzing GWAS data.

Statistical analyses were conducted using R software (version 4.3.1) (https://www.r-project.org). Survival analysis was performed using the Kaplan-Meier method and comparisons were made using the log-rank test. Pearson correlation tests evaluated relationships between variables. P < 0.05 was considered statistically significant.

The Guangzhou cohort comprised 28 ESCC patients, including 39.3% (11/28) achieving pCR and 60.7% (17/28) exhibiting npCR. The clinicopathological features of these patients are summarized in Supplementary Table 1. Survival data for patients who underwent NRCT and surgical resection and survived at least 3 months after treatment was collected and performed Kaplan-Meier analysis to evaluate their survival probabilities in both groups. Patients who achieved pCR demonstrated a higher survival probability than those with npCR, particularly after 10 years (Figure 1A, n = 26). However, due to the limited sample size, the differences between the groups were not statistically significant. To identify key chemoradiotherapy sensitivity-related genes, we analyzed DEGs between pCR and npCR patients. We identified 1,726 DEGs, with 870 genes upregulated and 856 genes downregulated in npCR patients. The heatmap and volcano plot illustrating these DEGs are presented in Figures 1B, C.

To further identify key chemoradiotherapy sensitivity-related genes in ESCC, we conducted univariate Cox regression analysis using RNA-seq and prognostic data from the Beijing Cohort (Supplementary Table 1, n = 119). This analysis identified 126 prognosis-associated protein-encoding genes, comprising 52 associated with a good prognosis and 74 indicating a poor prognosis. After selecting genes highly expressed in pCR patients that indicate a good prognosis and genes lowly expressed in pCR patients indicating a poor prognosis, we identified 35 favorable and 40 unfavorable genes. Following random survival forest analysis, genes with a relative importance >0.7 were deemed key chemoradiotherapy sensitivity-related genes in ESCC. Ultimately, three genes, ATF2, SLC27A5, and ALOXE3, met the screening criteria and are illustrated in Figures 2A, B. Among these, SLC27A5 and ALOXE3 were not only highly expressed in pCR patients (P < 0.01 and P = 0.01, Figure 2C), but also significantly correlated with a favorable prognosis of ESCC (P < 0.001, Figures 2E, F). However, ATF2 was highly expressed in npCR patients (P = 0.001, Figure 2C) and was significantly associated with poor prognosis (P = 0.001, Figure 2D). By employing the X-tile software, the optimal expression level cut-off values of three genes were determined to be 12.5, 10.63 and 7.89, respectively.

By expanding the real-world sample (n = 50), we reanalyzed RNA-seq data from ESCC patients who underwent NCRT followed by surgery. Overall, their mean age was 58.2 years, 43 (86%) were male, 7 (14%) were female, 14 (28%) had stage II disease, and 36 (72%) had stage III disease. The expression of ATF in npCR patients was significantly upregulated compared with that in pCR patients, while the expression of SLC27A5 and ALOXE3 was significantly downregulated (Supplementary Figure 1A, B), which was consistent with our previous results.

We analyzed immune cell infiltration between cancerous and normal tissues to assess immunological variations in ESCC. The proportions of infiltrating immune cells in each group were displayed in Figures 3A, B. Compared to normal tissue, cancerous tissue in ESCC showed increased infiltration of memory B cells, activated memory CD4⁺ T cells, resting NK cells, M0 macrophages, M1 macrophages, resting dendritic cells, and activated dendritic cells and decreased infiltration of naive B cells, CD8⁺ T cells, follicular helper T cells, regulatory T cells, activated NK cells, monocytes, activated mast cells, and eosinophils (Figure 3C). Correlations between infiltrating immune cells in ESCC were illustrated in Figure 3D. Further analysis of the associations between the expression of the three key genes and tumor immune cell infiltration revealed that ATF2 negatively correlated with the infiltration of memory B cells, monocytes, and resting dendritic cells, among others. SLC27A5 negatively correlated with the infiltration of memory B cells, resting memory CD4⁺ T cells, and M0 macrophages, among others. ALOXE3 negatively correlated with the infiltration of regulatory T cells, activated NK cells, and resting mast cells, among others (Figure 4A).

Considering that cancer patients often exhibit immune abnormalities related to tumor immune escape mechanisms, and the necessity to tailor immunotherapy targets and strategies based on individual immune characteristics (Sadun et al., 2007), we explored the associations between three key genes and immunoregulators, including 24 genes encoding immunoinhibitors and 46 genes encoding immunostimulators. Results indicated that ATF2 was significantly positively correlated with TNFRSF25, TNFRSF14, and ADORA2A, and negatively correlated with CD40, CD70, ENTPD1, and HAVCR2. SLC27A5 was significantly positively correlated with CD274, and negatively correlated with TNFRSF9 and TGFBR1. ALOXE3 was significantly positively correlated with CD40 and TNFSF14, and negatively correlated with CXCL12 (Figures 4B, C, all p < 0.01).

The “pRRophetic” R package assessed the potential of key genes to predict drug sensitivity in ESCC patients. Compared with patients exhibiting higher expression of ATF2 and lower expression of SLC27A5 and ALOXE3, IC50 values for drugs including vinorelbine, paclitaxel, docetaxel, fluorouracil, cisplatin, oxaliplatin, erlotinib, gefitinib, and lapatinib, were lower in patients with decreased expression of ATF2 and increased expression of SLC27A5 and ALOXE3. Specifically, patients with higher ATF2 expression were significantly less responsive to erlotinib and gefitinib, whereas those with increased ALOXE3 expression were significantly more sensitive to vinorelbine, paclitaxel, docetaxel, fluorouracil, erlotinib, gefitinib, and lapatinib. Patients with higher SLC27A5 expression were more sensitive to oxaliplatin (Figures 5A–I, P < 0.05). This further confirmed the clinical utility of these key genes in ESCC patients.

Using the GEPIA database (http://gepia.cancer-pku.cn), we analyzed the differential expression of key chemoradiotherapy sensitivity-related genes in tumor samples and paired normal tissues across 32 cancers, including ESCC. The results showed that ATF2 is highly expressed in the tumor tissues of diffuse large B-cell lymphoma, esophageal carcinoma, pancreatic adenocarcinoma, stomach adenocarcinoma, and thymoma compared to normal tissues (Supplementary Figure 2A). However, the expressions of SLC27A5 and ALOXE3 were not significantly upregulated in esophageal cancer tissues compared to control tissues (Supplementary Figure 2B, C).

We conducted Cox regression analyses on clinical factors to identify independent prognostic factors. Due to the limited number of TNM stage I cases (n = 6) in the Beijing cohort, we combined stages I and II for analysis. The nomogram was then constructed based on TNM staging and the expression levels of ATF2, SLC27A5, and ALOXE3 to quantitatively predict 1-, 3-, and 5-year survival probabilities, providing a reference for clinical decision-making in ESCC patients. The results showed that all three chemoradiotherapy sensitivity-related genes had a greater impact on prognosis prediction than TNM staging, with SLC27A5 contributing the most (Figure 6A). Calibration curves for 1-, 3- and 5-year OS demonstrated high consistency between the predictions and actual observations, underscoring the prognostic predictive power of these genes in ESCC (Figures 6B–D).

Utilizing hallmark gene sets from the Molecular Signatures Database, we conducted fast gene set enrichment analysis (fGSEA) to discern the differences in biological processes between high and low expressions of key genes. For the high-expression ATF2 group, the top three upregulated pathways were epithelial-mesenchymal transition (EMT), mitotic spindle, and myogenesis, while the top three downregulated pathways were oxidative phosphorylation, fatty acid metabolism, and KRAS signaling DN. Additionally, the p53 pathway was also downregulated (Supplementary Figure 3A, B). Conversely, the high-expression SLC27A5 group showed downregulation in the EMT pathway (Supplementary Figure 3C, D). Moreover, compared to the high-expression ATF2 group, the high-expression ALOXE3 group exhibited increased activity in the p53 pathway and the EMT pathway (Supplementary Figure 3E, F). These findings suggest that chemoradiotherapy sensitivity-related genes may influence ESCC progression by regulating the epithelial-mesenchymal transition and P53 pathways.

Using GWAS data, we identified the pathogenic regions of key genes in esophageal cancer (Supplementary Figure 4A, B). ATF2, SLC27A5, and ALOXE3 were found in the pathogenic regions of chromosomes 2, 19, and 17, respectively (Supplementary Figure 4C–E).