Association between cancer-associated fibroblasts and prognosis of neoadjuvant chemoradiotherapy in esophageal squamous cell carcinoma: a bioinformatics analysis based on single-cell RNA sequencing

Abstract

Background

Esophageal squamous cell carcinoma (ESCC) is a prevalent and aggressive subtype of esophageal cancer, posing a significant mortality and economic burden, especially in East and Southeast Asia. Current therapeutic strategies have limitations in improving patient survival, particularly regarding disease progression and resistance. This study aimed to investigate the impact of neoadjuvant chemoradiotherapy (NCRT) on the ESCC microenvironment.

Methods

We utilized single-cell RNA sequencing to systematically characterize the tumor and cancer-associated fibroblasts (CAFs) subtypes. Marker genes of myofibroblastic CAFs (myCAFs) were employed to establish a prognostic model and verify its application in other datasets. Other experiments were conducted on clinical samples to explore potential ESCC risk-related genes.

Results

Our bioinformatics and statistical analyses revealed an increased proportion of fibroblasts and epithelial cells in NCRT and identified the Ep_c1 subtype associated with a better prognosis. Further results indicated a complex communication network between Ep_c1 and myCAFs. The top 30 marker genes of myCAFs were used to construct a prognostic signature with a significant response to immunotherapy. Finally, experiments identified Complement C1s subcomponent (C1S), Decorin (DCN), and Neuroblastoma suppression of tumorigenicity 1 (NBL1) as potential ESCC risk-related genes.

Conclusion

Our findings highlight the dynamic alterations in the post-NCRT ESCC microenvironment and provide a foundation for the development of personalized treatment and immunotherapeutic approaches. Future studies are warranted to further validate these findings and explore their clinical implications.

Background

Esophageal cancer (EC) is one of the most diagnosed and deadly cancer types and comprises two epidemiologically and pathologically distinct diseases, esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC), of which ESCC accounts for more than 90% cases in China [1, 2]. More than half of the ESCC patients are already in advanced stages at the time of first diagnosis. Rapid local invasion and extensive metastases prevent patients from undergoing radical surgery, which is currently the only clinical approach for curing ESCC [3]. Currently, neoadjuvant or definitive chemoradiotherapy (NCRT or DCRT) is used for patients with locally advanced (including resectable and unresectable) ESCC who refuse or are unfit for surgery. However, the clinical outcomes of these patients remain unsatisfactory, with five-year overall survival (OS) rates ranging from 2.0 to 37.6% according to CONCORD-3 [4]. The recurrence rate after NCRT can reach 30–40% [5, 6]. Preoperative neoadjuvant therapy in conjunction with surgical resection represents the standard treatment for locally advanced ESCC. However, there is a lack of established strategies to identify ESCC patients at high risk of recurrence following standard NCRT and to improve their prognostic outcomes.

Recently, targeted therapies have achieved significant advancements across various malignancies, with particular emphasis on the noteworthy application of immune checkpoint inhibitors (ICIs) in ESCC [7,8,9,10]. In the second-line treatment of advanced esophageal cancer, the administration of programmed death receptor 1 (PD-1) inhibitor monotherapy has shown better antineoplastic activity than chemotherapy with controllable safety [11]. Notably, several recent studies have also explored first-line immunotherapy in advanced patients and demonstrated its effectiveness [9, 10]. Furthermore, immunotherapy combined with chemoradiotherapy (CRT) is currently being investigated as a novel treatment option for locally advanced ESCC with a latent synergistic effect and better outcome. Increasing evidence has demonstrated the pivotal role of the tumor microenvironment (TME) in the invasion and metastasis of ESCC. The TME plays a crucial role in promoting various processes associated with cancer progression, including cell proliferation, epithelial‒mesenchymal transition, angiogenesis/lymphangiogenesis, immunosuppression, inhibition of apoptosis, invasion and metastasis [12, 13]. The dynamic interplay between tumor cells and TME components drives tumor growth. Consequently, targeting the TME has emerged as a promising strategy for developing novel anticancer drugs against ESCC. Components of the TME include immune cells, fibroblasts, endothelial cells, the extracellular matrix and other constituents. Among them, cancer-associated fibroblasts (CAFs) represent the most prominent component of the TME. The phenotypes, origins, and functions of CAFs vary and include diverse subtypes. The vital role of CAFs has been identified in various malignancies, including lung cancers and ESCC [14,15,16]. Mounting evidence supports the utilization of distinctive features exhibited by CAFs as prognostic indicators for lung cancer and ESCC because of their different subtypes and the expression of several prognostic-related markers [17,18,19]. Moreover, activation of CAFs through multiple signaling pathways can promote tumor growth, angiogenesis, invasion and metastasis, as well as extracellular matrix (ECM) remodeling [20]. The secretion of a plethora of cytokines, growth factors and chemokines by CAFs can confer resistance to chemotherapeutic agents, thereby promoting tumor progression [21]. However, the mechanisms by which CAFs regulate antitumor immune responses in solid tumors are currently not fully understood.

Compared with conventional bulk RNA-seq analysis conducted on mixed cell populations, the latest advances in single-cell RNA sequencing (scRNA-seq) have facilitated a comprehensive characterization of intertumor (tumor-by-tumor) and intratumor (within a tumor) heterogeneity across various malignant tumors [22]. Notably, heterogeneity has been revealed in most ESCC interstitial cell types, particularly between fibroblasts and immune cells [23]. In addition, various tumor lineages and TME cells potentially associated with therapy resistance have been identified in ESCC [24]. The application of single-cell transcriptome sequencing has established a fundamental framework elucidating the key associations between cancer cells and diverse noncancer cells in the TME, which facilitates further exploration of ESCC progression and prognosis [25]. Consequently, additional research is warranted to investigate the precise role of interactions between tumor cells and CAFs in promoting the invasion and metastasis of ESCC.

Despite many scRNA-seq studies on ESCC, the systematic characteristics of tumor cells and CAFs in NCRT patients remain poorly understood, as do their potential interactions and mechanisms. In this study, we combined scRNA-seq data with bulk RNA sequencing data to characterize the subtypes of tumor cells and CAFs in NCRT patients with ESCC. Furthermore, we selected the subtypes associated with prognosis to analyze the potential interactions and molecular mechanisms between tumor cells and CAFs that may contribute to CRT resistance. Additionally, we identified a CAF-based prognostic signature for ESCC, validated its prognostic value, and analyzed the associated immune landscape and response to immunotherapy. Our findings offer new insights into the changes in the TME of ESCC patients treated with NCRT, supporting the combination of NCRT and immunotherapy for locally advanced ESCC.

Methods

Data collection and processing

The scRNA-seq data of ESCC patients in the GSE221561 [24] dataset were downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database, and the data of patients treated with CRT or surgery alone (SA) were downloaded from the clinical information of the samples. The scRNA-seq dataset GSE221561 contains scRNA-seq data of tumor tissues from 9 ESCC patients, which included 4 NCRT tumor samples (S1265T, S1535T, S2423T, and S6829T) and 2 SA tumor samples (S6417T and S9478T). The samples were integrated via the anchor method in the R package “Seurat”, and core cells were obtained by filtering via scRNA-seq. Ineligible cells include genes that can be detected in only 3 or fewer cells, and low-quality cells with fewer than 200 genes detected were excluded from subsequent analysis (min.cells = 3, min.features = 200). Quality control was then performed according to the data source of the original literature with the following requirements: nFeature > 600, nCount_RNA > 500, complex > 0.8, mitochondrial content < 15%, and Red BC content < 0.1%. After quality control, 6 samples were fused to generate 27,758 genes derived from 34,850 cells. The fused data were normalized by LogNormalize (scaling factor = 10000) followed by differential analysis to identify 2000 highly variable genes screened by the variance-stabilizing transformation (VST) method and then scaled to all genes via ScaleData. Principal component analysis (PCA) was performed on single-cell samples, and the top 17 principal components (PCs) were selected for subsequent analysis. Harmony was subsequently used to correct for batch effects. The uniform manifold approximation and projection (UMAP) algorithm was used to perform an overall dimensionality reduction analysis on the 17 PCs. With the R package “singleR”, HumanPrimaryCellAtlasData and BlueprintEncodeData were used as reference data for auxiliary annotation, and some known marker genes (Table S1) were used for manual annotation of different cell clusters. The expression of the marker genes was evaluated in all the ESCC samples.

The bulk expression profile data of GSE53625 [26] and GSE165252 [27] and the corresponding clinical information were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) to establish a prognostic signature. The GEO bulk data processing standard of the dataset is as follows: if the probe is converted into a symbol according to the probe correspondence of each platform, one probe corresponds to multiple genes, one probe is removed, and the average of multiple probes corresponding to the same symbol is taken. The FPKM expression profiles, OS information, and clinical information of the TCGA-Esophageal cancer (ESCA) cohort were downloaded from the Xena database (https://xena.ucsc.edu/). Samples with incomplete survival information or clinical information were excluded to obtain a training set of 164 ESCA patients for this study.

Identification and characterization of ESCC cell and CAF subtypes

In 4 NCRT samples (S1265T, S1535T, S2423T, and S6829T), epithelial cells were extracted and reprocessed via the “Seurat” package (resolution = 0.1 for clustering). After normalization, uniformization, identification of highly variable genes, removal of batch effects, PCA, clustering and identification were performed. Tumor cell subtypes were identified, and the marker gene expression of each subtype was evaluated. The malignant epithelial cells were obtained by deleting normal epithelial cells including endothelial cells and T/Natural killer (NK) cells as references via inferCNV (copy number variation) (R package: infercnv_1.14.2). The human genome (hg38.fa) and genome annotation files (gencode.v35.annotation.gtf) were downloaded from the UCSC Genome Browser Home (https://genome.ucsc.edu/). FindAllMarkers was used to identify the marker genes (avg_log2fc > 0.25, p_val_adj < 0.05) of each cluster.

Fibroblasts were extracted from NCRT samples in the same way as epithelial cells, and the 12 main components were selected for subsequent cluster analysis (resolution = 0.1 for clustering). We evaluated the results of using a resolution of 0.02-1.00 for clustering and found that the results of using 0.1/0.2 were reasonable and consistent. We finally also chose a resolution of 0.1 for clustering fibroblasts. The marker genes of fibroblasts, including myofibroblastic CAFs (myCAFs), inflammatory CAFs (iCAFs), antigen-presenting CAF (apCAFs) and extracellular matrix CAFs (eCAFs), used in previous studies [28,29,30,31] (Table S2) were used to identify the CAF subtypes. After identification, the marker gene expression for each subtype (avg_log2fc > 0.25, p_val_adj < 0.05) was evaluated, and the proportion distribution of different subtypes in each sample was visualized.

Gene set enrichment analysis and functional enrichment analysis of the marker genes of ESCC cell subtypes

Gene set enrichment analysis (GSEA) was performed on the marker genes of each ESCC cell subtype via the “clusterProfiler” package to explore their functions and associated pathways. Gene sets from the Kyoto encyclopedia of genes and genomes (KEGG), Gene ontology biological process (GOBP) and HALLMARK databases were downloaded from the Molecular Signature Database (MSigDB) (http://www.gsea-msigdb.org/gsea/index.jsp) and used to perform GSEA. Gene set variation analysis (GSVA) and single-sample gene set enrichment analysis (ssGSEA) were performed on the marker genes from each ESCC cell subtype, and the differences in ssGSEA scores among different subtypes were analyzed via the “limma” package. We selected 8 key pathways related to tumorigenesis and tumor progression for analysis and compared the ssGSEA scores of different subtypes.

Systematic characterization of the communication between tumor cells and CAFs based on ligand‒receptor interactions

CellPhoneDB version 5.0.0 was used to identify the ligand‒receptor interactions between tumor cells and CAFs in NCRT samples. On the basis of the chemokine gene family built into the R package “ktplots”, the marker genes utilized for the identification of CAF subtypes, as well as the specific marker genes for different CAF subtypes, the communication between tumor cell subtypes and CAF subtypes was comprehensively characterized and visualized via the R package “ktplots”.

Construction and validation of a prognostic signature

Based on the marker genes of each tumor cell and CAF subtype (logFC > 1.5, p_adj < 0.05), the ssGSEA algorithm of GSVA in the R package was used to calculate the enrichment score of each TCGA-ESCA sample. Univariate Cox regression analysis was used to determine the hazard ratio (HR) and prognostic significance of each tumor cell subtype.

According to the analysis of cell interactions, the CAF subtype that strongly interacts with the tumor cell subtype significantly related to prognosis was selected for the construction of a prognostic signature. The marker genes of the ESCC cell subtype (logFC > 1.5, p_adj < 0.05) and myCAFs (AUC > 0.7) were intersected to obtain the top 30 signature genes. The ssGSEA algorithm was subsequently used to calculate the signature scores of all the samples. The samples were divided into high- and low-score groups according to the median score. Survival curves were generated via the Kaplan‒Meier method, and the significance of the difference was determined via the log-rank test. In addition, univariate Cox analysis was performed to explore the independent prognostic value of the signature score.

Immune microenvironment analysis

To analyze the immune characteristics of the different signature score groups, we used ssGSEA with the R package ‘gsva’ to obtain 28 immune cell infiltration statuses for each sample in the TCGA-ESCA cohort. The Wilcoxon test was used to assess the differences in expression between the high- and low-score groups. The Pearson correlation coefficient was used to assess the relationships between the signature score and immune infiltrating cells, as well as between the signature genes and immune infiltrating cells.

Patients and specimens

We collected 20 frozen ESCC and adjacent normal tissues from Jinling Hospital, Medical School of Nanjing University, for qRT‒PCR. Prior to surgery, all patients underwent gastroscopy, biopsy, and upper abdomen contrast-enhanced computed tomography (CT) to confirm the diagnosis of ESCC, which included 10 patients who underwent NCRT followed by surgery and 10 patients who received surgery alone. All the samples were subsequently stored at -80 °C until used. In addition, 9 ESCC tissues for multiplex immunofluorescence staining were obtained from Jinling Hospital, Medical School of Nanjing University. Among these patients, 3 underwent surgery following NCRT combined with immunotherapy, 3 underwent surgery after NCRT, and 3 underwent surgery. All the samples were fixed in 4% paraformaldehyde and subsequently embedded in paraffin. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Jinling hospital ethics committee (Medical School of Nanjing University).

RNA isolation and qRT‒PCR assay

Total RNA from 20 ESCC tissues and paired paracancerous tissues was isolated using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The complementary DNA (cDNA) was synthesized according to the manufacturer’s instructions, utilizing the SuperScript III RT Reverse Transcription Kit (EXONGEN, Chengdu, China). qRT‒PCR was performed with a Sybr qPCR mix kit (Applied Biosystems, Carlsbad, CA, USA) on a StepOne Real-Time PCR system (Longgene, Hangzhou, China). The relative gene expression levels were quantified by employing the 2-△△CT method. The expression value of the target gene was normalized to that of the internal control gene Actin (CT Actin). All primers used in this research are listed in Table S3.

Multiplex immunofluorescence staining

Multiplex immunofluorescence staining was employed to investigate the expression of 3 selected genes in tissues from 9 ESCC patients. Two consecutive rounds of staining were carried out on each tumor section. The primary antibodies used in this study include Decorin (DCN) rabbit polyclonal antibody (1:200, Proteintech, #14667-1-AP), Neuroblastoma Suppression of Tumorigenicity 1 ( NBL1) rabbit polyclonal antibody (1:100, Affinity, #DF3177) and Complement C1s subcomponent (C1S) rabbit polyclonal antibody (1:200, Proteintech, #14554-1-AP). The signal was detected according to the manufacturer’s protocol using the PDOne TSA-RM-8275 kit (PANOVUE, Guangzhou, China), and ethylenediaminetetraacetic acid (EDTA) (pH = 9) repair mixture was also used. After two rounds of staining, the sections were counterstained with 4’,6-diamidino-2-phenylindole (DAPI) (PANOVUE, Guangzhou, China). Multiple stained sections were scanned with a Pannoramic SCAN II (3DHISTECH Ltd., Budapest, Hungary). Caseviewer 2.4.0 (3DHISTECH Ltd.) was used to analyze and acquire all the images.

Statistical analysis

All the statistical analyses were performed using R software (version 4.1.0) or GraphPad Prism 5.0. The Wilcoxon test was used for comparing two groups, whereas Spearman or Pearson correlation was used for correlation matrices. Survival differences through K‒M curves were assessed using the log-rank test. The differences in RNA expression between two groups were tested using Student’s t test. P value < 0.05 was considered statistically significant.

Results

Microenvironment cell landscape of post-NCRT and non-NCRT ESCC samples

To elucidate the changes in the TME of ESCC patients with NCRT, we used single-cell transcriptome data from GSE221561, which contains 4 NCRT samples and 2 SA samples, for analysis. Following the restriction of quality control, all samples were pooled for analysis (Figure S1). After quality control screening, 34,850 cells containing 27,758 genes were obtained. PCA was performed on all single-cell samples, followed by harmony correction to reduce batch effects (Fig. 1A, B). The UMAP plot was also followed by harmony correction to reduce batch effects (Fig. 1C, D). The core cells were classified into 12 independent cell clusters using the UMAP algorithm (Fig. 1E, F). The UMAP results were color-coded according to different cell clusters, which revealed that there was a significant difference between the NCRT samples and the SA samples (Fig. 1G). The different clusters were annotated via the “singleR " package, authoritative markers, and references [32], resulting in 7 cell clusters, namely, epithelial cells (n = 18,874), T/NK cells (n = 8743), endothelial cells (n = 1317), myeloid cells (n = 3060), fibroblasts (n = 1326), B cells (n = 1252) and mast cells (n = 278). After the cells were subgrouped, differential analysis was conducted to identify cluster markers. All the markers were screened using FindAllMarkers (min.pct = 0.25, logfc.threshold = 0.25). The 5 genes with the highest expression levels in each cluster were subsequently selected as the marker genes defining each cell type (Fig. 1H). Moreover, the bar charts of the cell cluster proportions revealed a discrepancy between the NCRT and SA samples in that the SA samples presented a greater abundance of T/NK cells, whereas the NCRT samples presented an increased presence of fibroblasts and epithelial cells with potentially more tumor cells (Fig. 1I).

Fig. 1

Identification of cell clusters according to scRNA-seq of all ESCC samples. A, B PCA plots before/after batch effects corrected using the harmony. C, D The UMAP plots before/after batch effects corrected using the harmony with clustering according to different samples. E, F The UMAP algorithm was applied to the top 12 PCs for dimensionality reduction, and 7 cell clusters were successfully classified in all ESCC samples. G The UMAP algorithm was used to analyze the cell composition of all ESCC samples, and it was found that there were significant differences between NCRT samples and SA samples. H Expression levels of top5 marker genes for each cell cluster. I The bar charts of the cell cluster proportion showed the discrepancy between the NCRT and SA samples

Identification of tumor cell subtypes and their marker gene functional enrichment analysis

Considering that malignant cells that survive after neoadjuvant therapy may acquire therapeutic resistance, NCRT samples were selected for identifying the cell subtypes associated with prognosis and resistance. A total of 18,036 epithelial cells were extracted from the NCRT samples. By using the Louvain algorithm, five epithelial cell subtypes were identified, namely, Ep_c0 to Ep_c4 (Figure S2A-C). To provide a concise overview of each subtype, marker genes were identified using the FindAllMarkers. Notably, Ep_c3 exhibited high expression levels of cell cycle-related genes and was thus designated the cycle-cell subtype (Ep_cycle) (Figure S2D, F). The expression of the cell cycle genes provided by the “Seraut” package is shown in Figure S2F. T-cell-related markers, such as CD3D, CD3E, CD3G, and PTPRC, were highly expressed in Ep-c4, which was subsequently excluded from further analysis (Figure S2D, E). Furthermore, infer CNV was employed to identify malignant epithelial cells into 5 groups (k_obs_groups = 5), with endothelial cells and T/NK cells serving as references. The inferCRV analysis yielded Figure S2G, which demonstrated that the result is consistent with the previous Louvain algorithm. Specifically, the cluster comprising the least measured cells with minimal variation (503) was classified as nonmalignant cells. After excluding Ep_c4 and nonmalignant cells, a total of 17,072 malignant epithelial cells were identified across 4 ESCC cell subtypes (Ep_c0:5080, Ep_c1:4110, Ep_c2:3968, and Ep_cycle:3914) (Fig. 2A), followed by enrichment analysis. By using FindAllMarkers and the Wilcoxon test, 10,161 significantly different marker genes were identified (Table S4). The marker genes of each subtype were identified based on the descending order of avg_log2FC. The expression of important marker genes for each subtype was visualized by bubble diagrams and heatmaps (Fig. 2B-D).

Fig. 2

Identification of tumor cell subtypes and evaluation of their marker gene expression. A The cell distribution of each subtype. B Bubble diagram of the top5 marker gene expression of subtypes. C Heatmap of the top10 marker gene expression of subtypes. D Bubble diagram exhibited the expression levels of cell-cycle genes for subtypes

To investigate the signaling pathways affected by the marker genes of different ESCC cell subtypes, GSEA was performed on each subtype, and the KEGG and GOBP results are displayed (Fig. 3A, B). The GSEA results based on KEGG and GOBP revealed that Ep_c0 was significantly enriched in epithelial cell differentiation and development, keratinization, and other pathways; Ep_c1 was significantly enriched in angiogenic, collagen metabolism, protein digestion and absorption, and other pathways; Ep_c2 was significantly enriched in drug metabolism, xenobiotic metabolism, cytokine‒cytokine receptor interactions, chemical carcinogenesis, fatty acid metabolism and other metabolic pathways; and Ep_cycle was significantly enriched in oocyte meiosis, the cell cycle, cell replication, DNA repairation, oxidative phosphorylation, clockgenes and other pathways. The GSVA results are presented in a heatmap (Fig. 3D). The pathways significantly upregulated in Ep_c0 included KRAS_SIGNALING_DN, IL2_STAT5_SIGNALING and APICAL_SURFACE. The pathways significantly upregulated in Ep_c1 included IL6_JAK_STAT3_SIGNALING, PROTEIN_SECRETION, UV_RESPONSE_DN and ANGIOGENESIS. The pathways significantly upregulated in Ep_c2 included FATTY_ACID_METABOLISM and PEROXISOME. The pathways significantly upregulated in Ep_cycle included MITOTIC_SPINDLE, ADIPOGENESIS and G2M_CHECKPOINT. Furthermore, the ssGSEA scores of each cell cluster in eight critical pathways of interest were calculated and compared (Fig. 3C). The results revealed that Ep_c1 had greater activity in ANGIOGENESIS than in other subtypes: Ep_c2 in FATTY_ACID_METABOLISM; Ep_cycle in OXIDATIVE_PHOSPHORYLATION; TGF_BETA_SIGNALING; ADHERENS_JUNCTION_ASSEMBLY and CLOCKGENES.

Fig. 3

Functional enrichment analysis of marker genes of tumor cell subtypes. A NES heat map base on KEGG enrichment analysis of tumor cell subtypes. B Top10 significantly enriched GOBP pathway in each subtype. C The ssGSEA enrichment scores of each subtype in eight critical pathways of interest. D HALLMARK pathway enrichment results of subtypes. (Wilcox. Test, *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.)

Identification and characterization of CAF subtypes

The NCRT samples were selected for clustering the CAF subtypes in the same way as the identification of tumor cell subtypes. After log normalization and dimensionality reduction, the obtained fibroblasts (n = 1218) were divided into 6 clusters (clusters 0–5) (Fig. 4A, B, C). According to marker gene expression (Figure S3), clusters 0 and 2 were identified as myCAFs (n = 678), cluster 1 as iCAFs (n = 290), cluster 3 as eCAFs (n = 145), cluster 5 as apCAFs (n = 44), and cluster 4 as the fib_ cycle (n = 61) (Fig. 4D). By using FindAllMarkers, all marker genes for each CAF subtype are shown in Table S5. The top 5 marker genes of each subtype are displayed in descending order according to avg_log2FC (Fig. 4E, F). The cell scale map (Fig. 4G) revealed that myCAFs were the subtype with the greatest number of cells in all the NCRT samples. Moreover, the expression of seven other marker genes that were not expressed or were expressed at lower levels in normal fibroblasts than in CAFs [28] was evaluated to better characterize activated CAFs (Fig. 4H, I). All the marker genes except SMN1 were highly expressed in myCAFs (Fig. 4H).

Fig. 4

Identification of CAF subtypes and evaluation of their marker gene expression. A, B The UMAP plot of the NCRT samples before/after batch effects corrected using the harmony. C The UMAP plot of six fibroblast clusters. D According to the marker gene expression, the distribution of CAF subtypes in samples were visualized by using the UMAP plot. E Heatmap of the top5 marker gene expression of each subtype. F The subtype distribution in the NCRT samples. G Bubble diagram of the top5 marker gene expression of each subtype. H, I Bubble diagram and distribution map of the 7 key marker gene expression in each subtype

Systematic characterization of the ligand‒receptor interactions between ESCC cells and CAF subtypes

The CellPhone DB was used to identify the ligand‒receptor pairs between ESCC cells and CAF subtypes in the NCRT samples that exhibited extensive communication (Fig. 5A). First, communication was analyzed using the chemokine gene family embedded in the “ktplot” package. Many significant interaction pairs, such as C-X-C motif ligand (CXCL)12-dipeptidyl peptidase-4 (DPP4) (Ep_c0-myCAFs), CXCL2-DPP4 (Ep_c0, c1, c2& cycle-myCAFs) and CXCL12- C-X-C chemokine receptor (CXCR4) (myCAFs-Ep_cycle), were detected between ESCC cell subtypes and myCAFs (Fig. 5B). Furthermore, the marker genes of fibroblasts and CAF subtypes were used to analyze their interactions. A wider range of ligand‒receptor interactions between ESCC cell subtypes and myCAFs was discovered, including interactions between collagen family members and integrins (Fig. 5C). Moreover, the WNT family also extensively interacts with secreted frizzled-related protein 2 (SFRP2) and frizzled related protein (FRZB), both of which are secreted frizzled-related proteins and are considered inhibitors of the WNT signaling pathway. In addition, a significant interaction between the platelet-derived growth factor (PDGF) family and its receptor was observed. When stimulated by tumor-associated mediators such as PDGF, CAFs are transformed into myCAFs, which regulate the progression of cancer cells by secreting a variety of cytokines and reshaping the extracellular matrix.

Fig. 5

Cell communication analysis. A Heat map of significant ligand-receptor pairs. B, C Bubble diagram of ligand-receptor interaction between tumor cell and CAF subtypes, based on the chemokine gene family and marker genes of CAF

Construction of a prognostic signature based on cell subtype interactions

The ssGSEA scores of each TCGA-ESCA sample were calculated based on the marker genes of the tumor cell and CAF subtypes. The ssGSEA scores for Ep_c0, Ep_c1, Ep_cycle, myCAFs, apCAFs, eCAFs and fib_cycle were significantly greater in tumor samples than in normal samples (p < 0.05) (Fig. 6A, B).

Moreover, univariate Cox regression analysis revealed negative correlations between the ssGSEA scores based on the marker genes of each subtype and prognosis, with no statistical significance (p > 0.05, hazard ratio (HR) < 1) (Table S6). Furthermore, according to the median ssGSEA score of each subtype, the TCGA-ESCA samples (n = 164) were separated into high- and low-score groups, and the relationships between the ssGSEA score and prognosis were analyzed. The results revealed that the high-score group of Ep_c1 (p = 0.0057) was associated with a better prognosis, whereas the high-score group of Ep_c2 was associated with a worse prognosis (p = 0.044) (Fig. 6C-F). Considering these findings collectively, ssGSEA scores, which are based on marker genes of Ep-c1 that are highly expressed in tumor samples and that are significantly negatively correlated with prognosis, were chosen for further analysis.

Finally, we selected the myCAF with the largest number of cells and a wide range of interactions for further investigation. The marker genes (n = 120) of myCAFs that significantly interact with Ep_c1 were chosen as candidate genes for constructing our prognostic signature. The ssGSEA score of each sample was calculated using the top 30 signature genes (AUC > 0.7) (Table S7) as the signature score, and the samples were then separated into high- and low-score groups by the median. K‒M survival analysis revealed that patients in the high-score group had a better prognosis than those in the low-score group did (Fig. 6G). Univariate Cox regression analysis revealed that the signature score was a significant independent prognostic factor for ESCC (p < 0.05) (Fig. 6H).

Fig. 6

A novel prognostic signature constructed based on cell subtype interaction. A, B Based on the marker genes of different subsets (A, tumor cell subtypes; B, CAFs subtypes), ssGSEA scores of TCGA-ESCA samples were calculated and compared between normal samples and tumor samples. C-G Kaplan–Meier curve result. (C, Ep_c0; D, Ep_c1; E, Ep_c2; F, Ep_cycle; G, myCAFs). H Univariate Cox analysis of signature scores and clinical characteristics. (Wilcox. Test, ns, P ≥ 0.05; *, P < 0.05; ****, P < 0.0001.)

Validating the application of the prognostic signature

To examine the robustness of our signature in predicting OS, we tested the prognostic value of the signature in a completely independent TCGA dataset with RNA-seq platform (GSE53625). The ESCC samples in the GSE53625 cohort were separated into high- and low-score groups by the median. K‒M survival analysis revealed that patients in the high-score group had a significantly better prognosis than those in the low-score group did in the GSE53625 cohort (Fig. 7A). Univariate Cox regression analysis confirmed that the signature score was a significant independent prognostic factor (p < 0.05) (Fig. 7B).

To predict the response of our signature to immunotherapy, we examined another dataset, GSE165252, derived from a clinical phase II PERFECT trial combining a PD-L1 inhibitor (atezolizumab) with NCRT in patients with resectable esophageal adenocarcinoma [27]. The samples in this clinical trial were categorized into three groups: the baseline group, the on-treatment group and the resection group. The signature scores of responders and non-responders in each group were calculated and compared. As a result, while the scores of responders were lower than non-responders in the baseline and on-treatment samples, the scores of those of responders were obviously greater than those of non-responders, with a statistically significant difference in the resection samples (Fig. 7C). In addition, as neoadjuvant therapy progressed, the proportion of non-responders with low scores increased, while the proportion of responders with high scores also increased. Ultimately, the non-responders accounted for 100% of the samples with low scores, whereas the responders represented 40% of the samples with high scores (Fig. 7D).

Fig. 7

Validation of the prognostic signature. A Kaplan–Meier curve of the signature scores in the GSE53625 cohort. B Univariate Cox analysis of signature scores and clinical characteristics in the GSE53625 cohort. C Differences in the signature scores between immune therapy non-responders and responders in the GSE165252 cohort. (a, baseline group; b, on treatment group; c, resection group). D Distribution of immune therapy response among different groups in the GSE165252 cohort. (a, baseline group; b, on treatment group; c, resection group) (Wilcox. Test, ns, P ≥ 0.05; *, P < 0.05.)

Immune infiltration analysis

To further investigate the predictive efficacy of our prognostic signature for immunotherapy, TSGA-ESCA samples were further analyzed. The immune landscape, as depicted in the heatmap, revealed that samples with high scores presented increased immune cell infiltration (Fig. 8A). Additionally, the relationships between signature genes and immune infiltration were investigated. Figure 8B shows that the signature genes were significantly positively correlated with most immune cells. These findings suggest that patients in the high-score group may derive greater benefits from immunotherapy, especially anti-PD-1 receptor blockers. Furthermore, on the basis of the expression of immune cell genes [33], we detected significant differences in the infiltration of the following immune cells between the high-score and low-score groups (p < 0.05): activated dendritic cells, CD56bright natural killer cells, central memory CD4 T cells, central memory CD8 T cells, effector memory CD4 T cells, effector memory CD8 T cells, gamma delta T cells, immature B cells, immature dendritic cells, macrophages, mast cells, MDSCs, memory B cells, natural killer T cells, plasmacytoid dendritic cells, regulatory T cells, T follicular helper cells, type 1 helper cells, and type 2 T helper cells (Fig. 8C).

Fig. 8

The immune infiltrations analysis. A Heatmap of immune cells infiltrations in high-and-low-score groups. B Correlation analysis between top30 characteristic genes and immune cells. C Comparison of proportions of 28 immune-related cells between high-and-low-score groups. (Wilcox. Test, ns, P ≥ 0.05; P < 0.05;**, P < 0.01;***, P < 0.001.)

The experiments of signature genes in clinical samples

To explore potential ESCC risk-related genes, ten genes involved in the signature were selected for further validation in twenty ESCC patients (n_NCRT=10, n_SA=10). Among the samples from the SA group, matrix metalloproteinase-2 (MMP2), lumican (LUM) and S100A13 presented significantly elevated expression levels in tumors, whereas DCN, fibulin 1 (FBLN1), NBL1, microfibril associated protein 4 (MFAP4) and lipoprotein receptor-related protein 1 (LRP1) presented significantly reduced expression levels in tumors (Fig. 9A-I). Moreover, in the NCRT group, MMP2, LUM and C1S exhibited significantly elevated expression levels in tumors, whereas DCN and NBL1 exhibited significantly reduced expression levels in tumors. Unfortunately, RNA expression of retinoic acid receptor responder 2 (RARRES2) was not detected in our ESCC samples after many attempts (Fig. 9A-I). These distinctions suggest that these genes may serve as innovative biomarkers for early ESCC diagnosis.

More notably, the difference in gene expression between non-NCRT tumor and adjacent tissues was distinct from that between post-NCRT tumor and adjacent tissues. The high expression of these genes in post-NCRT patients may be related to resistance to NCRT.

In this study, 3 genes (C1S, DCN & NBL1) of interest that are rarely reported in ESCC were selected to investigate their expression changes after NCRT or NCRT + immunotherapy. The results revealed that the expression levels and colocalization ranges of these genes varied among different tumor samples (NCRT + immunotherapy > NCRT > SA) (Fig. 9J-L).

Fig. 9

The expression of signature genes in clinical samples with or without neoadjuvant therapy. A-I The expression of MMP2 (A), LUM (B), DCN (C), FBLN1 (D), C1S (E), NBL1(F), MFAP4 (G), LRP1 (H) and S100A13 (I) in normal tissues and ESCC tissue of patients with or without NCRT. J-L Multiplex immunohistochemistry of pre- and post-therapy (NCRT or NCRT + immunotherapy ) samples using the antibodies and colors are as follows: DAPI (blue), C1S (green), DCN (red) and NBL1 (yellow). J, K & L were the multiple immunofluorescence results of SA sample, NCRT sample and NCRT + immunotherapy sample respectively. (two-sided unpaired Student’s t-test, ns, P ≥ 0.05; *, P < 0.05; **, P < 0.01; ***, P < 0.001.)

Discussion

Recently, the morbidity and mortality rates of esophageal cancer have significantly decreased, while patients with locally advanced ESCC continue to experience high recurrence rates following NCRT and surgery. High recurrence rates are often linked to resistance against NCRT, which is closely related to the heterogeneity of tumors [34, 35]. Advances in scRNA-seq can be used to characterize the heterogeneity of ESCC comprehensively and offer new strategies for precise diagnosis and treatment. Many studies have demonstrated a dynamic interaction between tumor cells and stromal cells, which play a significant role in tumor progression [36]. By elucidating the mechanisms of these interactions, it is possible to develop enhanced treatment strategies that simultaneously target multiple components of TME, ultimately improving patient outcomes [37]. CAFs are the predominant stromal cells in the TME and play a significant role in promoting tumor metastasis [38]. However, studies on the interactions between ESCC tumor cells and CAFs are still scarce, highlighting the importance of further exploration in improving our understanding of and therapeutic approaches for ESCC. Therefore, in this study, we conducted a comprehensive analysis of bulk RNA-seq and scRNA-seq data to identify distinct tumor cell and CAF subtypes and reveal their prognostic relevance, which in turn suggests a novel prognostic signature with excellent predictive efficacy for the immunotherapy response in ESCC patients.

First, our research on the microenvironment cell landscape provides critical insights into the cellular composition of the tumor microenvironment in ESCC and highlights the significant impact of NCRT on this composition. We compared the differences in the microenvironmental cell landscape between NCRT and SA ESCC samples. The proportion of T/NK cells in the SA samples was much greater than that in the NCRT samples, whereas the proportion of epithelial cells was significantly lower than that in the NCRT samples, which suggested a greater degree of malignancy and immunosuppression status in the NCRT samples. In a previous study on ESCC [24], the proportion of T cells observed in different treatment samples was SA > Ani-angiogenesis > Chemoradiotherapy > Adjacent normal > Chemotherapy, whereas the proportion of epithelial cells > Adjacent normal > Chemotherapy > Chemoradiotherapy > Ani-angiogenesis > SA > Chemotherapy, which is consistent with our results. Notably, the increased proportion of epithelial cells and decreased proportion of T/NK cells in the NCRT samples compared with the SA samples suggested a greater degree of malignancy and immunosuppression status in the NCRT samples. Moreover, the tumor cells in the NCRT samples accounted for 94.66% of the epithelial cells, which were probably derived from cell lineages resistant to NCRT. In addition, as a key factor in tumor development, the proportion of fibroblasts in NCRT samples was also significantly greater than that in SA samples, which also provides conditions for future studies. These findings indicate significant heterogeneity within the cell subtypes of NCRT and SA samples, which is likely associated with resistance and prognosis. This observation warrants further investigation into NCRT samples.

Next, we performed a comprehensive scRNA-seq analysis of ESCC samples with NCRT, which identified 5 tumor cell subtypes. Three of the subtypes were ultimately selected for univariate Cox regression, which revealed that the high-score groups of Ep_c0 and Ep_c1 were associated with better prognosis (p[Ep_c0] = 0.58, p[Ep_c1] = 0.0057), whereas the high-score group of Ep_c2 was associated with poor prognosis (p = 0.044), reflecting the internal heterogeneity of ESCC. The marker genes of Ep_c1 included spondin 2 (SPON2), membrane metallo-endopeptidase (MME), dynamin 3 opposite strand (DNM3OS), regenerating islet-derived 1 alpha (REG1A) and AC100801.1. SPON2 can promote the progression and metastasis of colorectal cancer [39], gastric cancer [40], and lung adenocarcinoma [41], but there are no relevant reports on its role in ESCC. MME suppresses the metastasis of ESCC by inhibiting the FAK-RhoA signaling axis [42]. The marker genes of Ep_c2 included carboxylesterase 1 (CES1), ATP-binding cassette, sub-family A, member 4 (ABCA4), regucalcin (RGN), glutathione-S-transferase mu-1 (GSTM1) and AC022706.1. CES1 belongs to a large mammalian serine esterase family [43]. When the activity of CES1 is blocked, lipid signaling interferes with the CES1-PPARα/γ-SCD axis, which sensitizes hepatocellular carcinoma (HCC) cells to cisplatin treatment [44]. Variation in ABCA4is associated with the neoadjuvant cytotoxic chemotherapy response in breast cancer [45]. RGN was found to act as a suppressor protein in carcinogenesis, and high expression of endogenous RGN was thought to have preventive and therapeutic effects on carcinogenesis [46]. However, a recent study suggested that high expression of RGN was associated with poor OS in patients with lung squamous cell carcinoma and that patients with low expression of RGN were more sensitive to chemotherapy drugs such as docetaxel and gemcitabine [47]. Taken together, although Ep_c1 and Ep_c2 have opposite prognostic significance, their marker genes both contain risk factors and protective factors, which provides evidence for the internal heterogeneity of tumors and a clue to therapy resistance.

Furthermore, our results revealed that the proportion of myCAFs was the highest among the five CAF subtypes in all the NCRT samples. It is also noteworthy that the expression of fibroblast activation protein (FAP) and podoplanin (PDPN) was greater in myCAFs than in other CAF subtypes. FAPs influence tumor growth via multiple mechanisms, including promoting proliferation, invasion, angiogenesis, epithelial-to-mesenchymal transition, stem cell promotion, immunosuppression and drug resistance [48]. In the TME, PDPN expression is upregulated in the tumor stroma, including CAFs and immune cells. CAFs play significant roles in extracellular matrix remodeling and the development of an immunosuppressive TME. Additionally, PDPN functions as a coinhibitory molecule on T cells, indicating its involvement in immune evasion [49, 50]. These reports suggest that myCAFs may be closely related to tumor immunosuppression and drug resistance.

Moreover, CellPhoneDB analysis revealed many ligand‒receptor interactions between myCAFs and ESCC cell subtypes. The interactions between myCAFs and Ep_c1 include mainly CXCL2-DPP4, PDGF family-PDGFR, and WNT family-SFRP2. CXCL2 is a chemokine from the CXC family whose expression in the tumor microenvironment attracts myeloid-derived suppressor cells (MDSCs) to the tumor, where chemokines that increase cancer cell survival are produced [51]. A novel study suggested that the expression of CXCL2 was greater in cisplatin-resistant bladder cancer cells than in their parent strains [52]. It has been reported that CXCL2 and CXCL10 can induce β-arrestin-1 recruitment to atypical chemokine receptor 2 (ACKR2), whereas ACKR2 can internalize and reduce the availability of CXCL10 in the extracellular space. The activity of CXCL10 toward ACKR2 was drastically reduced by DPP4 N-terminal processing [53]. These findings suggest a potential role for the CXCL2‒DPP4 interaction pair in the regulation of the TME. The functions of the PDGF/PDGFR pathway in tumorigenesis, tumor progression and metastasis, as well as its therapeutic targeting, have been reviewed extensively. SFRPs containing SFRP2 constitute a family of extracellular Wnt signaling antagonists [54]. Aberrant activation of Wnt signaling has been associated with the pathogenesis of virtually all human cancers [55]. Aberrant SFRP2 methylation is one of the major mechanisms for Wnt signaling activation in colorectal cancer (CRC). A recent study demonstrated that the beneficial association of tumor SFRP2 hypomethylation with non-chemoradiotherapy CRC was dependent on patient Body Mass Index (BMI), suggesting a possible tumor suppressor role for SFRP2 in overweight and obese patients [56]. Thus, the interactions between myCAFs and Ep_c1 play an important role in tumor progression and represent a potential resistance mechanism to NCRT in ESCC. By integrating previous studies with our results, it becomes evident that future experiments aimed at verifying the role of these interactions in resistance mechanisms are essential. Our systematic characterization of these significant interactions between myCAFs and ESCC cell subtypes highlights them as priority targets for subsequent studies investigating NCRT resistance mechanisms in ESCC.

Additionally, the marker genes of myCAFs were used to construct a prognostic signature that was significantly correlated with patient prognosis and response to immunotherapy. The immune infiltration analysis were performed to further explore the link between our prognostic signature and immunotherapy. And the results suggest that most immune cell populations are increased in the high-score samples (Fig. 8A). These identified immune cell populations include activated, immature, and plasmacytoid dendritic cells, CD56bright natural killer cells, as well as central and effector memory CD4/CD8 T cells (Fig. 8C). Many of these cell populations are mechanistically tied to anti-tumor immunity. Activated dendritic cells facilitate antigen presentation and limit tumor growth by stimulating CD8 + and CD4 + T cells that recognize tumor-derived neoantigens [57]. In addition, central memory and activated T cells are correlated with durable immune memory and improved responses to ICIs [58, 59]. Conversely, two highly plastic cell populations (neutrophils and Type 17 T helper cells) were found to be enriched in the low-score group without statistical significance (Fig. 8C), which suggests that they may primarily play immunosuppressive and tumor-promoting roles in this TCGA cohort. Neutrophils can directly promoting tumor progression, metastasis, and angiogenesis by releasing enzymes including myeloperoxidase (MPO), neutrophil elastase (NE), MMPs [60]. Type 17 T helper cell can promote tumor growth by means of angiogenesis (production of VEGF and angiogenic chemokines) and/or through immunosuppression (regulatory T cells conversion and myeloid-derived suppressor cell recruitment) [61]. The strong correlation between signature genes (e.g., DCN, FBLN2, CCDC80, C1S) and cytotoxic lymphocyte (T and NK cells) infiltration suggests that these biomarkers may reflect an immunologically “hot” tumor microenvironment. This provides a rationale for combining CAF-targeted therapies with PD-1 blockade in high-score patients, as CAF reprogramming could enhance pre-existing immune responses [62, 63].

Finally, in clinical samples, the top 10 signature genes with the highest AUC values were further analyzed. The qRT‒PCR results revealed that the expression trends of all the genes were the same in the different tumor samples (SA or NCRT) compared with the adjacent normal samples. The genes were differentially expressed between NCRT and adjacent normal samples, which might be related to resistance to NCRT. These genes included MMP2, LUM and C1S, whose expression was significantly upregulated, as well as DCN and NBL1, whose expression was significantly downregulated in tumor tissues. Cai et al. showed that reticulocalbin3 could promote the expression of MMP-2 by regulating the inositol triphosphate receptor 1 (IP3R1)-Ca²⁺-calcium/calmodulin-dependent protein kinase II-c-Jun signaling pathway, thereby promoting the progression and platinum resistance of ESCC [64]. Luo et al. demonstrated that carbon ion beams could regulate CDH1 and MMP2 downstream of the signal transducer and activator of transcription 3 (STAT3) pathway and inhibit ESCC cell metastasis, which activated the STAT3 signaling pathway [65]. These studies suggest that ESCC patients with high MMP2 expression may be resistant to platinum-based chemotherapy drugs and sensitive to radiotherapy, and such patients have a high possibility of recurrence after NCRT. Several previous studies have demonstrated that LUM and DCN can be used as biomarkers for the early detection of ESCC¹. In addition, NBL1 has been found to be downregulated in ESCA patients, which was associated with poor OS rates in patients relative to those in the high-expression group [68]. Ge et al. demonstrated that the downregulation of C1S inhibits ESCC cell proliferation through the Wnt1/β-catenin pathway and promotes ESCC cell apoptosis by regulating the expression of Bcl2, Bax and cleaved-caspase3 [69]. Combined with the results of the immune infiltration analysis, we selected C1S, DCN and NBL1 to evaluate their expression after different neoadjuvant therapies. The results revealed that the expression levels of these markers were significantly elevated in the NCRT samples compared with those in the SA samples, albeit slightly diminished relative to those in the NCRT + immunotherapy samples. These observations indicate that Ep-C1 cells and myCAFs exist before neoadjuvant therapy and evolve into a resistant population following NCRT while still exhibiting a favorable response to immunotherapy. In addition, these marker genes hold promise as potential biomarkers for ESCC and may serve as predictors of response to neoadjuvant therapy. Furthermore, the preliminary exploration on the biological roles of these key genes in clinical samples underscores their potential as priority targets for future investigations into the mechanisms of NCRT resistance in ESCC.

Targeting CAFs has emerged as a prominent area in current targeted therapies [70] and serves as a potential effective complement to the treatment of ESCC. Current tumor treatment strategies targeting CAFs mainly include direct interventions, such as clearing CAFs [71] and inhibiting their activity [72], as well as indirect approaches that affect CAFs functions by inhibiting downstream effectors [73]and targeting ECM [74]. Moreover, CAFs mediate immune evasion of tumor cells [75]. Thus, therapeutic strategies to modulate CAFs activity may improve the responsiveness of ESCC to immunotherapy. For instance, targeting specific CAF-derived signals or depleting CAFs can reduce immunosuppressive effects and promote a more robust anti-tumor immune response [76, 77]. This is of great significance as the tumor microenvironment, in which CAFs play a pivotal role, can significantly impact the efficacy of immunotherapeutic approaches in ESCC. By targeting CAFs, we can potentially reshape the TME to be more conducive to the action of the immune system. In fact, some clinical trials has proven that CAF-targeting drugs such as the angiotensin inhibitor-losartan [78] and depleting hyaluronic acid-PEGylated recombinant human hyaluronidase (PEGPH20) [79] can enhance the drug anti-tumor effect of chemotherapy. In this study, distinct CAF subtypes in ESCC samples with NCRT and tumor cell subtypes with extensive connections to them were identified. And a prognostic signature correlated with prognosis and response to immunotherapy was successfully established based on the signature genes of the myCAFs. To some extent, these results implies the potential role of CAF-targeted therapy in NCRT or immunotherapy resistance towards ESCC.

Recently, a large number of clinical studies have been conducted on CAF-targeted therapies. It is noteworthy that some of these targeted drugs have synergistic potential with ICIs, such as FGFR inhibitors and MMP9 inhibitors. Fexagratinib (AZD4547), a FGFR inhibitor, the anti-tumor efficacy of who have been demonstrated in a clinical trial [80], was reported to synergize with ICIs in the HCC [81]. However, Andecaliximab (ADX), a MMP9 inhibitor, who theoretically enhances the delivery of chemotherapeutic and immunotherapeutic agents, was not demonstrated significant anti-tumor efficacy combined with chemotherapy in a recent phase III clinical study [82]. These studies demonstrate that CAF-targeted therapies has potential to synergize with ICIs in cancer treatment. However, more research is needed in the future to further explore the effectiveness and specific mechanisms of CAF-targeted therapies. The combination of CAF-targeted therapies and ICIs may lead to a more comprehensive and potent anti-tumor immune response, offering new hope for improving the treatment outcomes of patients. Overall, we hope that our study can provide some hints for exploring potential therapeutic targets in ESCC, and further exploration of CAF-targeted therapies is crucial for formulating more effective ESCC treatment strategies in the future.

Although our study offers valuable insights, it is important to acknowledge certain limitations. First, the signature was constructed based on retrospective data sourced from public databases (GEO&TCGA); therefore, further prospective and multicenter cohorts of patients with ESCC are essential to mitigate any potential bias. Second, the biological roles of these key signature genes and the significant interactions between myCAFs and ESCC cell subtypes in the context of NCRT resistance in ESCC remain to be elucidated. Future research should prioritize extensive in vivo and in vitro studies to explore these mechanisms. This represents a crucial direction for the ongoing development of our research.

Conclusions

In our study, we conducted an extensive investigation into the populations of tumor cells and CAFs in ESCC samples with NCRT to identify many distinct clusters. Furthermore, the ligand‒receptor interactions between ESCC cell subtypes and myCAFs were fully characterized. Top 30 marker genes of myCAFs were used to construct a prognostic signature. Moreover, our signature was associated with the immune landscape and is suitable for predicting the responsiveness of ESCC patients to immunotherapy targeting PD-L1 blockade. Overall, this study could be used as a reliable predictor of the efficacy of NCRT and propose a potential mechanism underlying treatment resistance in ESCC patients.