Introduction

Esophageal cancer (EC) is a major cause of cancer-related deaths globally, with an estimated 604,100 new cases and 544,076 deaths in 2020.1 In Eastern Asia, esophageal squamous cell carcinoma (ESCC) accounted for the majority of EC.

Recent evidence has suggested that patients with ESCC could gain survival benefits from multimodal neoadjuvant therapies.2–5 After neoadjuvant chemoradiotherapy (NCRT), pathological complete response (pCR) rates of 48.6% and 43.2% were observed in the CROSS trial (37 ESCC cases) and the NEOCRTEC5010 trial (451 ESCC cases), respectively.6 7 Compared with esophageal adenocarcinoma, ESCC was associated with a higher pCR rate after NCRT, although this was not an independent prognostic factor for long-term survival.8 In recent years, immunotherapy has revolutionized the treatment of cancer patients. Immune checkpoint blockade (ICB), such as anti-PD-1/PD-L1 and anti-CTLA-4, can activate T lymphocytes in patients with various cancers, which was first observed in melanoma in 2011.9–11 As expected, several clinical trials have been initiated to evaluate the safety and efficacy of ICB in patients with ESCC in both metastatic and curative settings. In the KEYNOTE-181 trial, ICB showed an advantage over chemotherapy in terms of overall survival (OS) and objective response rate 12; and we conducted the single-arm trial of PALACE-1 (20 ESCC cases, NCT03792347), and reported a 55.6% pCR rate after NCRT combined with ICB.13 Despite these encouraging results, not all patients benefit. Unfortunately, the response to treatment with ICB was not always correlated with the expression of PD-L1.14 15 An immuno-oncogenic classification was proposed lately, and EGIC1 group showed improved survival in ESCC patients with chemotherapy plus ICB treatment.16 Several studies explored the effects of neoadjuvant therapy on the tumor microenvironment (TME) in ESCC, showing that NCRT may augment the antitumor immune response by activating CD4+ and CD8+ T lymphocytes.17–19 A recently published paper reported that ICB therapy could promote the clonal revival of pre-existing and new infiltrating T cells, as revealed by RNA-seq and T cell receptor sequencing.20 However, a comprehensive study of immune TME in ESCC after neoadjuvant therapy at protein level was urgently required to explore the potential mechanism of treatment response and the biomarkers to inform treatment decisions.

Mass cytometry by time-of-flight (CyTOF) is a high-dimensional proteomics method that enables simultaneous analysis of more than 50 markers at the single-cell level using metal isotope-labeled antibodies. This technique is very useful in the context of tumor immunology, especially in studies dissecting TME alterations involved in the response to ICB.21

Here, based on the neoadjuvant cohort of ESCC patients, including our PALACE-1 trial of NCRT combined with ICB, we performed CyTOF in 26 ESCC samples from patients who underwent surgery with or without different types of neoadjuvant therapy. Our results revealed dynamic cellular phenotypes in the immune TME of ESCC reshaped by neoadjuvant therapies.

MethodsPatients and samples

Surgically resected specimens were collected from patients with resectable locally advanced ESCC from August 2021 to March 2023. The clinicopathological characteristics of all the patients are listed in online supplemental table S1 and S2. Treatment response was evaluated using the AJCC (American Joint Committee on Cancer)/NCCN (National Comprehensive Cancer Network) modified Ryan tumor regression grading (TRG) system.22

ESCC tumors were dissected immediately after specimen collection, and the viable tissue was washed with ice-cold normal saline. Normal esophageal mucosa was obtained from a site at least 5 cm away from the tumor margin and managed in a manner similar to the tumor samples. Surgical specimens were collected and stored in a tissue storage solution (Miltenyi Biotec, Germany). The specimens were shipped between 2°C and 8°C within 24 hours.

Sample digestion and preparation

Normal esophageal mucosae and ESCC tumors were dissociated into single-cell suspensions by combining mechanical dissociation with enzymatic digestion according to the instructions of the Tumor Dissociation Kit (Miltenyi Biotec, Germany). Cells were pelleted by centrifugation at 400×g for 10 min, washed, and resuspended in ice-cold PBS. Cell counting was performed using a hemocytometer.

Mass cytometry staining and data acquisition

A total of 3×106 cells per sample were collected from the single-cell suspension for further staining. Each sample was stained for viability with 0.25 µM cisplatin (Fluidigm, USA) in 100 µL PBS on ice for 5 min. The cells were then washed twice with the FACS buffer. After blocking with a 50 µL blocking mix for 20 min on ice, 50 µL surface antibody cocktail was added to each sample. Surface staining was performed on ice for 30 min. They were washed twice with FACS buffer and incubated in 200 µL Fix and Perm Buffer (Fluidigm, USA) containing 250 nM nucleic acid Intercalator-Ir (Fluidigm, USA) overnight at 4°C. The fixed cells were then washed with 1×Permeabilization Buffer (eBioscience, USA), and stained with an intracellular antibody cocktail on ice for 30 min. Stained cells were washed and incubated with combinations of up to three out of six barcoding palladium metals (103Pd, 104Pd, 105Pd, 106Pd, 108Pd, 110Pd, Fluidigm, USA) for 30 min at room temperature, washed, and pooled. Pooled samples were resuspended in deionized water and analyzed on a Helios mass cytometer (Fluidigm, USA). Details of the panel and antibodies are presented in online supplemental table S5.

Mass cytometry data analysis

A total of 26 ESCC tumor samples and 6 normal esophageal mucosa samples (normal group) were analyzed: 14 tumor samples from patients who received neoadjuvant therapy (neo group), while the other 12 did not (naive group). Among six normal tissues, only two of the six normal tissues were matched to the tumor samples in the naive group, as the other four lost their matched tumor samples during the quality control process. According to the types of neoadjuvant therapies, the neo group was divided into three subgroups: neo-c group (neoadjuvant chemotherapy, n=5), neo-cr group (NCRT, n=5) and neo-cri group (NCRT combined with ICB, n=4).

After the raw data were debarcoded and normalized, the FlowJo software V.10.4.0 (BD Biosciences, USA) was used to exclude the debris, dead cells, and doublets. Preprocessed data were then imported into R (V.4.1.2) and analyzed using the package cytofWorkflow (V.1.18.0). A total of 100,000 cells were selected from each sample by random sampling for subsequent analysis, whereas samples with cell numbers less than 100,000 were not downsized (online supplemental table S6). Following quality control, the cytofWorkflow function ‘cluster’, which is based on the FlowSOM algorithm, was used for unsupervised clustering with corresponding lineage and function markers. The clusters were then manually merged and annotated according to their specific marker expression. The data were visualized using the t-distributed stochastic neighbor embedding (t-SNE) algorithm via the cytofWorkflow function ‘runDR’ and ‘plotDR’, illustrating the distribution of clusters and marker expression. Marker expression of each cluster was demonstrated by heatmap and dotplot via the cytofWorkflow function ‘plotExprHeatmap’ and ggplot2 (V.3.3.5), respectively. In addition, some cell subsets were identified by traditional gating strategies using FlowJo, because of the lack of accuracy of the clustering algorithm. We did not include markers related to tumor cells, and they are thus clustered in cells that are negative for any lineage markers. Data from certain cell subsets were extracted for subsequent analyses.

Single-cell RNA-seq data analysis

The processed single-cell RNA-seq (scRNA-seq) data of ESCC with annotation files were retrieved as indicated in previous publications.23 24 The two datasets were treated as independent studies and were analyzed separately. The Seurat package (V.4.1.0) was applied for the reanalysis of the scRNA-seq data. Cell clusters were identified and visualized using the Seurat function ‘DimPlot’ by the coordinates of uniform manifold approximation and projection (UMAP) in the annotation files. Cluster marker expression was illustrated using the Seurat function ‘FeaturePlot’ and ‘DotPlot’. For dividing NK cell subsets without given annotation, variable genes were detected by the Seurat function ‘FindVariableFeatures’ and the expression data were then scaled by the Seurat function ‘ScaleData’. Principal components (PCs) were calculated using the Seurat function ‘RunPCA’. Top PCs were selected for clustering using the Seurat function ‘FindNeighbors’ and ‘FindClusters’, and the dimensionality reduction UMAP was calculated by the Seurat function ‘RunUMAP’. Cluster identification and marker demonstration were performed as described previously. Differential expression genes of CD16+ NK cells were analyzed by the Seurat function ‘FindMarkers’ and the upregulated genes (log2 fold change>0.25, p<0.05) were used for pathway enrichment analysis through the online toolkit Metascape.25 The gene ratio of pathway enrichment analysis equals the number of differentially expressed genes in a pathway set divided by the total gene numbers in that pathway set.

Flow cytometry staining and analysis

A total of 20 ESCC tumor samples (naive group, n=6; neo group, n=14) and 6 normal esophageal mucosa samples (normal group) were analyzed. Specifically, the neo group was consisted as follows: neo-c group (n=4), neo-cr group (n=5) and neo-cri group (n=5).

Each sample (3×106 cells) was incubated with 50 µL blocking mix for 20 min on ice. The diluted solution of surface antibodies and Fixable Viability Dye eFluor 780 (BioLegend, USA) was prepared, and each sample was resuspended with 100 µL solution and incubated at 4°C for 30 min. Following staining, cells were fixed using 1×IC fixation buffer (Invitrogen, USA) and incubated at 4°C for 40 min. Data were acquired on a LSRFortessa X20 cytometer (BD Biosciences, USA) and analyzed using FlowJo software.

Chemokine antibody array

Two ESCC tumor samples and two normal esophageal mucosa samples were collected independently from the mass cytometry samples. All tissue samples were homogenized mechanically and further analyzed using the Human Chemokine Array C1 according to the manufacturer’s instructions (RayBiotech, USA).

Immunohistochemistry

Forty-six formalin-fixed paraffin-embedded (FFPE) specimens of ESCC were retrospectively collected (naive group, n=10; neo group, n=36). The neo group was consisted as follows: neo-c group (n=9), neo-cr group (n=13) and neo-cri group (n=14). Their clinicopathological characteristics are listed in online supplemental table S3. A 4 µm section was stained with primary antibodies of CCR4/CCR6 system as specified in online supplemental table S5. H-score was evaluated by two experienced pathologists independently based on five random fields captured at 20×magnification. H-score was calculated as the formula: H-score=percentage of positive cells multiplied by intensity score (0–3).26

Immunofluorescence

Multiplex immunofluorescence analyses were conducted on paraffin-embedded tissue sections using Opal 3-Plex Detection Kit (Ayoka, USA) as described in the supplier’s protocol. Images were captured using ECLIPSE C1 microscope and NIKON DS-U3 system (NIKON, Japan).

The Cancer Genome Atlas data analysis

RNA-seq counts data and clinical information of 95 ESCC tumor samples and 11 normal esophageal mucosa samples from The Cancer Genome Atlas (TCGA) were downloaded from UCSC Xena. Normalization was conducted using the transcripts per million (TPM) method. Survival analysis was conducted using the R packages survival (V.3.2–13) and survminer (V.0.4.9). Kaplan-Meier survival curve was analyzed by the survival function ‘survfit’ and plotted by the survminer function ‘ggsurvplot’, and the p value was calculated by log-rank test. Dichotomous stratification could be based on single gene expression or gene signature scores, which were calculated by mean value of log2(TPM+1) of genes in a certain signature gene set.27

Prognosis model establishment

Seven genes in the CCR4/CCR6 chemokine system were used to establish a survival model. Ridge-penalized Cox (Ridge-Cox) regression, a machine learning method, was applied as all the genes were supposed to be preserved in the model establishment. The regression was computed using the R package glmnet (V.4.1-3). The regression fit was obtained using the glmnet function ‘cv.glmnet’ with fivefold cross-validation and λ.min was selected as the λ value in use. The predictive power of the model was demonstrated by a Kaplan-Meier survival curve. The regression fit was also adopted in the validation dataset from the PERFECT trial for the model score calculation. The median model score was used to separate high-risk and low-risk patients.

Statistical analysis

All statistical analyses were performed using GraphPad Prism (V.8.0) and R (V.4.1.2). We employed Student’s t-test to compare two groups, while for comparisons involving three or more groups, we utilized a one-way analysis of variance (ANOVA) followed by Tukey’s test for post hoc analysis. Statistical details are provided in the legends. A p<0.05 was statistically significant (*p<0.05, **p<0.01, ***p<0.001, ****p<0.0001). Bars in the graphs represent mean±SEMs.

ResultsSingle-cell mass cytometry analysis of ESCC in patients with or without neoadjuvant therapies

To elucidate the changes in the TME after neoadjuvant therapy in patients with ESCC, we collected 26 ESCC tumor samples and 6 normal mucosa samples for mass cytometry by CyTOF analysis using a 42-antibody panel (figure 1A). Following the restriction of quality control, we pooled all samples for analysis (online supplemental figure S1A). We obtained 524,924, 1,200,000, and 1,174,601 cells in the normal, naïve, and neo groups, respectively. The clinical features were comparable between naïve and neo groups, or between neo-c, neo-cr and neo-cri groups (online supplemental table S4).

Figure 1

Experiment overview and single-cell annotation of esophageal squamous cell carcinoma revealed by mass cytometry. (A) Study schema (created with BioRender.com). (B) t-SNE plot of single cells colored by major cell types. (C) t-SNE plots demonstrating the expression of marker genes, CD45, CD31, FAP and αSMA. (D) Frequency of major cell types in normal, naïve and neo groups. (E) Box plots showing the proportion of immune cells, CAFs, endothelial cells and unknown cells among normal, naïve and neo groups. P values were derived from one-way ANOVA, Tukey’s test; *p<0.05. Error bars: mean±SEM. ANOVA, analysis of variance; CAFs, cancer-associated fibroblasts; scRNA-seq, single-cell RNA-sequencing.

To characterize the cellular composition of the TME, we performed graph-based clustering using lineage markers. Four major cell subsets were identified: immune cells (n=914,328) with CD45 expression, cancer-associated fibroblasts (CAFs, n=734,496) with FAP or αSMA expression, endothelial cells (n=436,707) with CD31 expression and unknown cells (n=813,994) without lineage markers, because of the limited antibodies in the CyTOF panel (figure 1B,C, online supplemental figure S1B). The unknown cells could encompass tumor cells, normal epithelial cells, smooth muscle cells, neurons and any other cells that are negative for lineage markers. The proportions of immune cells, CAFs, endothelial cells and unknown cells were almost comparable among the normal, naïve, and neo groups (figure 1D,E, online supplemental figure S1C).

CAFs are one of the critical cell types in the TME that determines the oncological outcomes.28 To further investigate the CAF heterogeneity, we subdivided the CAFs into PDGFRα+ CAFs (n=161,541),29 CD73+ CAFs (n=339,850),30 31 which expressed the protein encoded by gene Ecto-5' nucleotidase (NT5E) and converting extracellular AMP to adenosine, and unknown CAFs (n=233,105) without clear feature markers (online supplemental figure S2A,B). PDGFRα+ CAFs produced a high level of the inflammatory cytokine CXCL8, similar to inflammatory CAFs described in a previous report (online supplemental figure S2B).32 The CAFs subset composition did not change after neoadjuvant therapy, regardless of the therapeutic strategy (online supplemental figure S2C,D).

Identification of 11 major immune cell types in ESCC microenvironment

Dynamic remodeling of the immune microenvironment during treatment may affect the clinical outcomes of ICB, chemotherapy and radiotherapy.33 34 Here, we performed graph-based clustering using lineage markers of immune cells to assign 11 cell types (figure 2A,B, online supplemental figure S3A). The corresponding markers of T cells (n=494,006), B cells (n=91,746), and NK cells (n=27,921) were CD3, CD19 and CD56, respectively (figure 2B,C, online supplemental figure S3A). Myeloid cells were initially clustered as a heterogeneous collection of cells with differential expression levels of CD68, CD14, CD11b, CD11c, and HLA-DR. To overcome the inability of the clustering algorithm to distinguish conventional dendritic cells (cDCs) from macrophages, traditional gating was adopted for cDCs recognition (CD14-CD11c+HLA-DR+, n=93,919), and the remaining myeloid cells were labeled as macrophages (n=136,475) (online supplemental figure S3B). Manual gating was used due to the lack of resolution for NKT cells as well, and they were separated from T cells as CD3+CD56+ cells (n=43,390) (online supplemental figure S3B). Granulocytes, which were often dropped out by scRNA-seq studies, here were thoroughly depicted. Neutrophils (n=10,918) and eosinophils (n=7055) were characterized by CD66b expression, while neutrophils were additionally marked by CD16. Basophils (n=1613) and mast cells (n=5138) both expressed FceRIa, but basophils showed preferential expression of CD123. In addition, plasmacytoid dendritic cells (pDCs, n=2147) were identified by the coexpressed markers CD123 and HLA-DR.

Figure 2

Characterization of immune cell clusters in ESCC revealing decreased dendritic cell percentage by neoadjuvant therapies. (A) t-SNE plot of immune cells colored by major immune cell types. (B) Dotplot showing the expression of marker genes that define each immune cell type. (C) t-SNE plots demonstrating the expression of marker genes, CD3, CD19, CD16, CD11b, CD14, CD11c, CD66b, CD123 and FceRIa. (D) Frequency of immune cell types in normal, naïve and neo groups. (E) Box plots illustrating the proportion of T cells, B cells, NK cells, macrophages, cDCs and pDCs among normal, naïve and neo groups. (F) Kaplan-Meier survival curves for overall survival of the ESCC patients in TCGA database according to the gene signatures of pDCs. P values in (E) were derived from one-way ANOVA, Tukey’s test; *p<0.05, **p<0.01. Error bars: mean±SEM. P values in (F) were derived from log-rank test. Baso, basophils; cDCs, conventional dendritic cells; Eos, eosinophils; Mac, macrophages; Mast, mast cell; Myeloid, myeloid cells; Neu, neutrophils; pDCs, plasmacytoid dendritic cells

pDCs frequency was diminished by neoadjuvant therapies and correlated with poor prognosis

We evaluated the frequency alterations in major immune cell types among the normal, naïve, and neo groups. As a result, cDCs and pDCs increased in the naïve group compared with the normal group and then significantly reduced after neoadjuvant therapy (figure 2D,E, online supplemental figure S3C). To investigate the functional role of cDCs and pDCs in ESCC, we calculated their gene signature scores in TCGA database, and explored whether the signature could differentiate the OS rates of ESCC (online supplemental figure S3D).35 36 cDCs signatures did not correlate with the survival rate of ESCC patients, but the pDCs signature separated the survival curve quite well (p=0.019) (figure 2F, online supplemental figure S3E). It was postulated that pDCs produce less type I interferon, promote a suppressive microenvironment and initiate tumor progression.37 The three neoadjuvant regimens diminished the pDC subset to a similar extent (online supplemental figure S4A).

Antitumor immunity of T cells was boosted by neoadjuvant therapies

The identified T cell included CD8+ T cells (n=227,192), regulatory T cells (Tregs, CD4+CD25+, n=58,268), conventional CD4+ T cells (Tconvs, CD4+CD25−, n=197,635) and γδ T cells (n=10,911) (figure 3A,B). Tregs were highly abundant in ESCC tumors relative to normal mucosa and were depleted by neoadjuvant therapy (figure 3C). In contrast, the proportion of CD8+ T cells was increased by neoadjuvant therapy, whereas the percentages of Tconvs and γδT cells were not significantly affected. Comparing the different neoadjuvant therapies, the neo-cr and neo-cri treatment groups had a stronger benefit in improving CD8+ T cell infiltration than the neo-c group (figure 3D). Different neoadjuvant strategies yielded a statistically similar result for Tregs depletion, but with a modestly higher Tregs frequency in the neo-c group (figure 3D).

Figure 3

Neoadjuvant therapies strengthening the anti-tumor immune responses of T cells. (A) t-SNE plot of T cells colored by T cell lineages. (B) t-SNE plots of T cells showing the expression of marker genes, CD4, CD25, CD8, TCRgd. Box plots illustrating the proportion of Tconvs, Tregs, CD8+ T cells and γδ T cells among (C) normal, naïve and neo groups and among (D) normal, naïve, neo-c, neo-cr and neo-cri groups. Box plots showing the percentage of naïve, TCM, TEM and TEMRA relative to the total number of Tconvs or CD8+ T cells among (E) normal, naïve and neo groups and among (F) normal, naïve, neo-c, neo-cr and neo-cri groups. P values were derived from one-way ANOVA, Tukey’s test; *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001. Error bars: mean±SEM. Tconvs, conventional CD4+ T cells; Tregs, regulatory T cells; TCM, central memory T cells; TEM, effector memory T cells; TEMRA, terminally differentiated effector memory T cells.

Tconvs and CD8+ T cells were assessed for their naïve and memory status, using the markers CCR7 and CD45RA. Naïve T cells were CCR7+CD45RA+, central memory T cells (TCM) were CCR7+CD45RA−, effector memory T cells (TEM) were CCR7−CD45RA−, and terminally differentiated effector memory T cells (TEMRA) were CCR7−CD45RA+ (online supplemental figure S5A,B). Both Tconvs and CD8+ T cells showed a decrease in TCM frequency and an increase in TEM frequency on neoadjuvant therapy, reversing the condition of high TCM/TEM ratio in ESCC tumors (figure 3E). This suggests a shift from TCM to TEM after neoadjuvant treatment. Specifically, chemoradiotherapy alone or in combination with ICB reformed the TEM/TCM ratio of Tconvs, indicating an induction of antigen exposure and T cell reactivation, while chemoradiotherapy alone significantly changed the TEM/TCM ratio of CD8+ T cells (figure 3F).

The chemokine receptors CXCR3, CCR4, and CCR6 were used for T helper cell subset classification. The gating strategies are shown (online supplemental figure S5C). Th2 and Th9 cells were enriched in the TME, but none of the T helper subsets were remodeled by neoadjuvant therapies (online supplemental figure S5D). Tissue resident memory T cells (Trm) were evaluated in Tconvs and CD8+ T cells by the marker CD103, but differences were not observed between naïve and neo groups (online supplemental figure S5D). Exhausted CD8+ T cells were marked by PD-1 expression, and the percentage of PD-1+CD8+ T cells and PD-1 expression level showed no significant difference after neoadjuvant therapy (online supplemental figure S5E). Moreover, the functional markers PD-1 and CD73 in Tregs did not change after neoadjuvant therapy (online supplemental figure S5E).

CD16+ NK cells were promoted by neoadjuvant ICB treatment

NK cell infiltration was significantly enhanced by chemoradiotherapy combined with ICB, and there was an increasing trend with chemoradiotherapy alone (figure 4A), verified by flow cytometry data (online supplemental figure S6A). CD16 (FCGR3A), an activation marker of NK cells, clearly divided NK cells into two clusters, CD16+ (n=20,149) and CD16− NK cells (n=7772) (figure 4B). To further validate these two clusters in scRNA-seq data, we obtained single-cell transcriptomic profiles of NKs from two published studies,23 24 including GSE145370 and GSE160269, and graph-based clustering grouped NK cells into NCR3+ NK cells, KLRC+ NK cells, and STMN1+ proliferating NK cells (figure 4C, online supplemental figure S6B). We found that CD16 was significantly expressed in NCR3+ NK cells, consistent with CD16+ NK cells in our CyTOF results, whereas KLRC+ NK cells and STMN1+ proliferating NK cells corresponded to CD16− NK cells (figure 4D,E, online supplemental figure S6C,D). To further investigate the functional roles of CD16+ and CD16− NK cells, we performed an enrichment analysis of differentially expressed genes (log2 fold change>0.25, p<0.05) between CD16+ and CD16- NK cells and found that CD16+ NK cells had superior cytotoxicity and enhanced activation status (figure 4F, online supplemental figure S6E). In addition, KLRC+ and STMN+ NK cells had high exhaustion scores and STMN+ NK cells were the most proliferative, as indicated by their G1/S and G2/M scores (online supplemental figure S6F). We demonstrated that ESCC tumor tissues had a diminished percentage of CD16+ NK cells compared with normal esophageal mucosa in the GSE145370 dataset23 (figure 4G), and ESCC with higher pathological stages had an even lower CD16+ NK cell ratio in the GSE160269 dataset24 (online supplemental figure S6G). A significant decrease or decreasing trend of the CD16+ NK cell ratio was observed in the naïve group (figure 4H,I). Comparing the naïve group to the three strategies of neoadjuvant treatment, we observed that neo-cri group remarkably enhanced CD16+ NK cell proportion, to which the proportional change in neo-cr group was less significant (figure 4H,I).

Figure 4

Neoadjuvant ICB combined with chemoradiotherapy enhancing the activation and cytotoxicity of NK cells. (A) Box plot showing the proportion of NK cells among normal, naïve, neo-c, neo-cr, and neo-cri groups. (B) t-SNE plot of NK cells colored by NK cell clusters. (C) UMAP plot of NK cells colored by NK cell clusters reanalyzed using the GSE145370 dataset. (D) UMAP plot showing the expression of FCGR3A in different subsets of NK cells from the GSE145370 dataset. (E) Dotplot illustrating the expression of marker genes of different NK cell clusters reanalyzed using the GSE145370 dataset. (F) Pathway enrichment analysis of upregulated genes in CD16+ NK cells using scRNA-seq data from the GSE145370 dataset. (G) Scatter plot demonstrating the proportion of CD16+ NK cell subset divided by the total NK cell number derived from the GSE145370 dataset. Box plots showing the frequency of CD16+ NK cell subset relative to the total NK cell number among normal, naïve, neo-c, neo-cr and neo-cri groups, demonstrated by (H) CyTOF and (I) flow cytometry analysis. P values in (A, H and I) were derived from one-way ANOVA, Tukey’s test; p values in (G) were derived from two-sided paired Student’s t-test; *p<0.05; **p<0.01; ****p<0.0001. Error bars: mean±SEM. ANOVA, analysis of variance; UMAP, uniform manifold approximation and projection

CCR4+CCR6+ macrophages were attenuated by neoadjuvant therapies

The percentage of total macrophages did not change after neoadjuvant therapy (figure 2E); therefore, we evaluated the M1/M2 status of macrophages based on CD86 and CD206 expression. However, these subsets were comparable among the normal, naïve, and neo groups (online supplemental figure S7A,B). Using graph-based clustering with macrophage-related markers, we found that chemokine receptors CCR4 and CCR6 were highly coexpressed in a rather independent macrophage subset (figure 5A). This new cell type was validated using immunofluorescence as CD68+CCR4+CCR6+ macrophages (figure 5B). CCR4+CCR6+ macrophages were modestly increased in the naïve groups and were significantly reduced after neoadjuvant therapy (figure 5C). Different neoadjuvant strategies (neo-c, neo-cr, and neo-cri) equally reduced the frequency of CCR4+CCR6+ macrophages (online supplemental figure S7C). We obtained similar results in additional samples throughflow cytometry analysis (online supplemental figure S7D–F). Furthermore, we explored the spatial relationship between CCR4+CCR6+ macrophages, T cells, and tumor cells. Our findings revealed a close proximity between CCR4+CCR6+ macrophage and T cells in the vicinity of tumor nests (online supplemental figure S7G), suggesting a potential regulatory role for this macrophage subset.

Figure 5

CCR4 and CCR6 centered signaling in ESCC reshaped by the neoadjuvant therapies. (A) t-SNE plot of macrophages demonstrating the coexpression pattern of CCR4 and CCR6. (B) Multiplex immunofluorescence showed the presence of CD68+CCR4+CCR6+ macrophages. (C) Box plots illustrating the proportion of CCR4+CCR6+ macrophages among normal, naïve and neo groups. (D) Heatmap showing the chemokine production of CCL3, CCL5, CCL17, CCL20, and CCL22 between normal esophageal mucosae and ESCC tumors detected by chemokine antibody array. (E) Box plots showing the expression levels of CCL3, CCL5, CCL17, CCL20, and CCL22 between normal esophageal mucosae and ESCC tumors from the TCGA ESCC dataset. Scatter plots comparing the H-score of CCL3, CCL5, CCL17, CCL20, CCL22, CCR4, and CCR6 between tumor regression grade (TRG) 0/1 and TRG 2/3 patients in (F) neo_cr and (G) neo_cri group. P values in (C) were derived from one-way ANOVA, Tukey’s test; p values in (E–G) were derived from two-sided unpaired Student’s t-test; *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001. Error bars: mean±SEM. ANOVA, analysis of variance; ESCC, esophageal squamous cell carcinoma; TCGA, The Cancer Genome Atlas.

Chemokines that paired with CCR4 and CCR6 were enriched in ESCC microenvironment

To investigate the regulatory factors of CCR4 and CCR6, we detected their ligands, CCL3, CCL5, CCL17, and CCL22 to CCR4 and CCL20 to CCR6,38 in the ESCC microenvironment using a chemokine antibody array. We found that all of the chemokine ligands to CCR4 and CCR6 were upregulated in the tumors (figure 5D, online supplemental table S7). The higher expression of CCL3, CCL5, CCL17, CCL20, and CCL22 in the ESCC tumors was confirmed by the transcriptomic results from TCGA database (figure 5E).

We then investigated the sources of each chemokine by analyzing the scRNA-seq data. CCL3 was generally produced by monocytes/macrophages, and CCL5 was mainly generated by CD8+ T cells and NK cells (online supplemental figure S8A–C). CCL20 was mainly expressed by monocytes/macrophages, and epithelial cells in the non-immune compartment (online supplemental figure S8A–D). Interestingly, CCL17 and CCL22 were both produced solely by a special cDCs subset, mregDCs (mature DCs enriched in immunoregulatory molecules) (online supplemental figure S8A–C). The mregDCs had both maturation and immunoregulatory features, characterized by the marker genes LAMP3, PD-L1 (encoded by CD274), and CCR7, and were also revealed by our analyses (n=1641) (online supplemental figure S9A,B).35 39 40 The percentage of mregDCs in the naïve group was higher than that in the normal group but only showed a modest decline in the neo group (online supplemental figure S9C,D). Except for mregDCs, cDCs were classified into cDC1s (n=74,073) and cDC2s (n=18,160) by the expression of FceRIa, but were not significantly changed by neoadjuvant therapies (online supplemental figure S9E–G). Also, we performed immunohistochemistry (IHC) analysis of the CCR4/CCR6 system in ESCC specimen with or without neoadjuvant treatment, and evaluated H-score of each protein. The protein expression of CCR4/CCR6 ligands and receptors was not much changed between the naïve group and three neo groups (online supplemental figure S10A,B). But in neo_cr and neo_cri groups, those proteins had a significantly lower expression in patients with a TRG of 0 or 1, which indicated a fair treatment response, than those with a TRG of 2 or 3 (figure 5F).

Establishment of a predictive model for the response of neoadjuvant therapies based on CCR4/CCR6 chemokine system

Considering that neoadjuvant therapy remodeled cell types associated with the CCR4/CCR6 chemokine system, we assumed that the gene signature of CCR4/CCR6 system could be related to the prognosis of patients with ESCC. To achieve this, we calculated TCGA data to verify the prognostic value of a seven-gene signature (CCR4, CCR6, CCL3, CCL5, CCL17, CCL20, and CCL22), and it was able to stratify the patients into two groups with significantly different OS (figure 6A). To improve the predictive value of the seven-gene signature, we used Ridge-penalized Cox (Ridge-Cox) regression with fivefold cross-validation and established a model (figure 6B,C). Survival prediction using the modified model signature showed clear classification (figure 6D). We further applied this predictive model to a test dataset derived from a clinical phase II trial of esophageal adenocarcinoma (PERFECT trial). In this study, patients received NCRT combined with anti-PD-L1 blockade.41 As a result, while the model scores of baseline samples were comparable between responders and non-responders, the scores were enhanced slightly during neoadjuvant treatment in the non-responders rather than the responders (online supplemental figure S10C). Furthermore, we analyzed the paired samples separately. Most of the non-responders experienced a significant increase in model scores compared with the on-treatment to the baseline samples, but the model scores of responders did not change (figure 6E). Based on the IHC analysis of post-treatment ESCC samples, we calculated a model H-score using the same coefficients as the transcriptomic model. Similarly, the model H-score was significantly higher in the patients with TRG of 2 or 3 than those with a TRG of 0 or 1, either in neo_cr or neo_cri group (figure 6F). In a nutshell, this indicates that the dynamic elevation of the CCR4/CCR6 chemokine-based model may reflect the unresponsiveness of neoadjuvant therapies and proposes a potential mechanism of treatment resistance.

Figure 6

CCR4/CCR6 chemokine-based model predicting the response to the neoadjuvant therapies of ESCC. (A) Kaplan-Meier survival curves for overall survival of the ESCC patients in TCGA database according to CCR4/CCR6 chemokine system. (B) Ridge coefficient profiles of se