Single-cell profiling of response to neoadjuvant chemo-immunotherapy in surgically resectable esophageal squamous cell carcinoma

Participants who had received immune therapy, such as monoclonal antibodies against programmed death-1 (PD-1)/PD-1 ligand (PD-L1) and cytotoxic T lymphocyte-associated antigen-4 (CTLA-4) or targeted therapies and radiotherapy were ineligible for this study. After screening, 22 patients who were histologically diagnosed with stage II to IV esophageal cancer between July 2020 and March 2021 and had not received the treatment mentioned above were included in this study. Relevant clinical information, such as sex, age, and treatment information, was retrieved from the database for each patient. The clinical stage was determined according to the American Joint Committee on Cancer (eighth edition) [32]. This study was conducted following the Declaration of Helsinki and was approved by the Institutional Ethics Committee of Tangdu Hospital, Fourth Military Medical University (No. 202102–23). Written informed consent was obtained from each patient before sample collection.

All 22 patients received neoadjuvant chemo-immunotherapy (NAT) comprising tislelizumab (BeiGene Co. Ltd, Shanghai, China) (N = 20) or camrelizumab (Jiangsu Hengrui Pharmaceuticals Co. Ltd) (N = 2), along with carboplatin/nedaplatin (Qilu Pharmaceutical Co. Ltd, Jinan, China) and albumin-bound paclitaxel (Qilu Pharmaceutical Co. Ltd, Jinan, China) for three cycles (3 weeks/cycle). Among these, 18 patients underwent curative-intent complete surgery following NAT.

Freshly excised samples were stored in the sCelLiVE® Tissue Preservation Solution (Singleron Biotechnology, Nanjing, China) on ice and transferred to the laboratory within 30 min after resection. All specimens were washed with Hanks Balanced Salt Solution (HBSS) three times, minced into 1-mm cubic pieces, and then digested with 3 mL sCelLiVE® Tissue Dissociation Solution (Singleron Biotechnology, Nanjing, China) at 37 °C for 15 min. The cell suspension was collected and filtered through a 40-μm sterile strainer. To remove red blood cells, the GEXSCOPE® red blood cell lysis buffer (RCLB) (Singleron Biotechnology, Nanjing, China) was added into the cell suspension at a ratio of 2:1 (RCLB: cell) and incubated at room temperature for 5 min. The mixture was then centrifuged at 300 × g for 5 min at 4℃. The supernatant was removed and the sediment was gently resuspended in phosphate-buffered saline (PBS) (HyClone, China).

scRNA-seq libraries were constructed as previously described [34]. In brief, single-cell suspensions (~ 2 × 10⁵ cells/mL) were loaded onto a microfluidic chip of the Singleron Matrix® Single Cell Processing System (Singleron Biotechnology, Nanjing, China). After the mRNA from each cell was labeled with a unique molecular identifier (UMI), single-cell suspensions were collected from the microfluidic chip, followed by reverse transcription and PCR amplification to obtain the scRNA-seq libraries. Individual libraries were diluted to 4 nM and paired-end sequenced on the NovaSeq 6000 sequencing system (Illumina, San Diego, USA).

The CeleScope (version 1.9.0, https://github.com/singleron-RD/CeleScope↗) pipeline was used to generate the gene expression matrix containing normalized gene counts versus cells per sample. The expression matrix was filtered for each sample dataset through the Seurat (v3.1.2) R toolkit [35]. The primary exclusion criteria are as follows: (1) cells with a gene count less than 200 or with a top 2% gene count; (2) cells with a maximum of 2% of UMI counts; (3) cells with mitochondrial content more than 50%; (4) genes expressed in less than five cells. The quality control data for scRNA-seq can be found in Additional file 2.

The filtered raw count matrix was normalized by dividing by the total counts per cell and transformed by calculating the natural logarithm. The top 2000 variable genes were selected using the “FindVariableFeautres” function in Seurat (v3.1.2) [35]. Principal component analysis was then performed, and the top 20 principal components were subjected to dimensionality reduction and subsequent cell type clustering using the “FindClusters” function. Batch effect correction was performed using Harmony v1.0 [36]. Cell clustering was achieved by the Louvain algorithm, setting the resolution parameter to 0.8, and visualized by the Uniform Manifold Approximation and Projection (UMAP) method [37]. Cell-ID (v0.1.0, https://github.com/RausellLab/CelliD↗) [38] was used for automated annotation of cell types, followed by manual correction of cell types by classical markers in the SynEcoSys database (Singleron Biotechnology, Nanjing, China) [39–41]. All Seurat functions were run with default parameters unless otherwise specified.

Inferred copy number variations (CNVs) were analyzed using the inferCNV tool (v1.2.1; https://github.com/broadinstitute/inferCNV↗) [42] using non-malignant cells as baselines. In brief, genes expressed in more than three cells were sorted by their genomic locations on each chromosome. A sliding window comprising 101 genes was used to smooth the relative expression on each chromosome to remove the effect of gene-specific expression. The relative expression values were centered at 1, and the ceiling of the relative expression values was set to 1.5 standard deviations from the residual normalized expression levels. The R package Pheatmap (v1.0.12; https://github.com/raivokolde/pheatmap↗) [43] was used to visualize patients’ inferred CNV profiles.

SLICE (https://research.cchmc.org/pbge/slice.html↗) was applied to determine the single-cell entropy (scEntropy) of cancer cells [44]. The differentiation of cancer cells was subsequently delineated using the scEntropy-directed pseudotime trajectory generated by the Monocle 2 method (v2.22.0, https://www.bioconductor.org/packages/release/bioc/html/monocle.html↗) [45]. The “FindVariableFeatures” in Seurat (v3.1.2) was used to select highly variable genes from clusters. DDRTree, a reversed graph embedding algorithm in Monocle 2, was used to reduce dimension and reconstruct the temporal and bifurcation structure of the database based on global gene expression levels. The trajectory was visualized by the “plot_cell_trajectory” function in Monocle 2.

Hotspot (v0.9.1, https://github.com/YosefLab/Hotspot↗) was used to identify functional gene modules that illustrate heterogeneity within cancer cell subsets [46, 47]. Briefly, we used the “danb” model and selected the top 500 genes with the highest autocorrelation z-score for module identification. Modules were then identified using the “create_modules” function with the parameters: min_gene_threshold = 15 and fdr_threshold = 0.05. Module scores were calculated by using the “calculate_module_scores” function. Jaccard similarity coefficients were used to evaluate the transcriptional similarity between 12 cancer modules and signatures of cancer cell clusters.

Single-cell regulatory network inference and clustering (SCENIC), a computational method for network inference and motif discovery allowing high-confidence prediction of key regulators and their direct target genes, was performed using the R package “SCENIC” (v0.1.5; https://github.com/aertslab/SCENIC↗) [48, 49]. The co-expression modules were run by GRNBoost. We downloaded the motifs database for Homo sapiens from https://pyscenic.readthedocs.io/en/latest/↗. The input matrix consisted of the normalized expression values of the cells of interest.

Differentially expressed genes (DEGs) were identified by the “FindMarker” function in Seurat based on the Wilcoxon rank sum test. Genes expressed in more than 10% of cells in a cluster and with an average log (fold change) greater than 0.25 and an adjusted P value less than 0.05 were identified as DEGs. The cutoff for DEGs remains the same unless otherwise specified.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed using the “clusterProfiler” package in R [50]. GO gene sets include three categories, namely biological process (BP), cellular component (CC), and molecular function (MF). Pathways with an adjusted P value less than 0.05 were considered significantly enriched.

The CellPhoneDB (v2.1.0), equipped with a known ligand-receptor pair database, was used to predict cell–cell interaction [51, 52]. The number of ineffective distributions of reciprocal exchanges was set to 1000 and otherwise followed the default settings of the software. The predicted interaction pairs with a P value less than 0.05 and average log (fold change) greater than 0.1 were considered significant.

The Cancer Genome Atlas (TCGA)-esophageal carcinoma (ESCA) dataset, which includes bulk RNA-seq data and overall survival information for 159 patients, was obtained from cBioPortal (https://www.cbioportal.org/↗). To evaluate the relationship between the target gene set and patient prognosis, single sample gene set enrichment (ssGSVA) was used to calculate the score of each gene set for samples as previously described [53]. Samples were assigned to either high or low-expression groups using the best score as the cutoff. Kaplan–Meier curves were used to show differences in the overall survival of patients between high and low-expression groups for each gene set using the “Survival” package in R (https://github.com/therneau/survival↗).

The unpaired two-sample Wilcoxon test was used to compare cell distribution within two groups of cells. An unpaired two-sample t-test was used to compare the mean gene expression or gene signatures within two groups of cells. The Jonckheere-Terpstra test was used to assess the statistical significance of the trends across IPR, MPR, and pCR groups. Log-rank test was used to compare the overall survival of subgroup patients. A two-sided P value of less than 0.05 was considered significant for all tests unless indicated otherwise (*P < 0.05, **P < 0.01, ***P < 0.001). All statistical analyses were performed in R (version 4.1.3).

Notably, four patients had an unknown pathological response status due to surgical cancellations. The rest of the 18 patients were categorized into three pathological response groups, including those with pathological complete response (pCR), major pathological response (MPR), and incomplete pathological response (IPR) (Fig. 1A; Additional file 1: Fig. S1A). The clinical characteristics of the remaining 18 patients with a definitive response to NAT are shown in Additional file 3: Table S1.

To describe the transcriptional expression atlas of ESCC at the single-cell level, we performed scRNA-seq on all 46 samples, including tumor and normal adjacent tissues before and after NAT (Fig. 1A). By comparing the expression of canonical cell-type marker genes, we identified 14 distinct cell subsets totaling over 250,000 cells, including cancer, immune (including T, B, plasma, mononuclear phagocytes (MP), erythrocyte and mast cells), stromal (including endothelial cells (EC), mural cells and fibroblasts), parenchymal (including basal, gland mucous cells, and gland ductal cells), and Schwann cells (Fig. 1B–D; Additional file 1: Fig. S1). In post-treatment ESCC tumors, we observed a notable reduction in the proportion of cancer cells, accompanied by an augmentation in the presence of T cells (Fig. 1B). These results indicate a dynamic change within the ESCC TME in response to treatment. Meanwhile, UMAP embeddings of pre-treatment ESCC tumors, stratified by pathological response, displayed extensive overlap, implicating a profoundly intricate and heterogeneous TME among these patients (Fig. 1E).

The dynamic characteristics and heterogeneity of malignant cancer cells were further investigated by trajectory analysis. To determine the start point of the inferred pseudotemporal cancer cell ordering, we first applied SLICE (Single Cell Lineage Inference Using Cell Expression Similarity and Entropy) to quantify single-cell entropy (scEntropy) of each cancer cell type [44]. As entropy inversely correlates with cell differentiation state, cells with lower scEntropy had more well-defined cell fates and functionalities compared to those clusters with higher scEntropy. In this context, clusters 1 and 11 exhibited a significantly lower entropy compared to others, marking the start point in trajectory mapping (Additional file 1: Fig. S2D). Subsequently, we used the “Monocle 2” method [56–60] to construct entropy-directed transitional trajectories of cancer cells (Fig. 2B). Notably, differentiation trajectory analysis revealed a branched structure of cancer cells, with several clusters of cancer cells positioned along the extended branch away from the start point, suggesting that these cell clusters were likely associated with higher functional uncertainty and greater differentiation potential (Fig. 2C). Interestingly, we observed that less differentiated clusters (higher entropy, such as clusters 2, 4, 5, 9, 10) were primarily identified in pre-treatment tumor samples of pCR and MPR patients, while more differentiated clusters (lower entropy, such as clusters 1 and 11) were predominantly found in post-treatment tumor samples (Fig. 2D). These observations suggest an interesting link between cancer cell differentiation and the pathological response to neoadjuvant chemo-immunotherapy. Specifically, ESCC tumors enriched with less differentiated cancer cells might be associated with a better pathological response to neoadjuvant chemo-immunotherapy, whereas tumors enriched with clusters of more differentiated cancer cells may resist the treatment. Likewise, the pseudotime faceted map of paired tumors categorized by response displayed a consistent pattern, with a better pathological response observed in tumors enriched with cancer clusters characterized by higher entropies (Additional file 1: Fig. S2E). Nevertheless, it is crucial to recognize that scRNA-seq inherently presents limitations in precisely calculating cell proportions, as demonstrated by the underrepresentation of cancer cells in post-treatment tumor samples from IPR patients. Consequently, these results should be approached with caution, and future validation is imperative.

Furthermore, we analyzed the transcriptomic changes of cancer cells associated with transitional states (Fig. 2E, F). Three phases of differentiation were identified. Cancer cells in state 3, which positively responded to NAT, exhibited significantly elevated expression levels of immune-related genes, including those implicated in neutrophil activation, antigen processing and presentation, and response to IFNγ (Fig. 2F). Conversely, cancer cells mapped to state 2 showed significantly increased expression of genes associated with mRNA catabolic processes. On the other hand, state 1 was predominantly occupied by cancer cells characterized by higher expression of genes involved in epidermal development. These results were consistent with the dynamic changes in gene expression of cancer cells during the transition (Additional file 1: Fig. S2F, G). Collectively, these findings suggest that state 3 likely encompassed cancer cells actively engaged in immune responses, whereas state 1 appeared to consist primarily of epithelial-like cells at the transient stage.

Hotspot analysis [46] was performed to explore informative gene modules characterized by distinct expression patterns and their correlation with pathological response to NAT. It is noteworthy that unsupervised clustering of cancer cells, as a data-driven method, facilitates the identification of subpopulations of cancer cells sharing similar expression patterns. On the other hand, hotspot analysis represents a hypothesis-driven approach, wherein specific gene sets and pathways are predefined based on prior knowledge or biological significance. Here, we identified 12 distinct gene modules corresponding to the 12 previously annotated cancer cell clusters in ESCC (Fig. 2G; Additional file 1: Fig. S2H). Cancer cluster 3 was predominantly observed in pre-treatment tumor samples of IPR patients and exhibited a strong correlation with module 3 (Fig. 2D; Additional file 1: Fig. S2H). However, it remains inconclusive whether the significant reduction of cluster 3 in post-treatment samples is correlated with the efficacy of NAT due to the limited representation of cancer cells in IPR_T_A samples (Fig. 1E). On the other hand, cancer cluster 1 with lower entropies was significantly correlated with module 1 (Jaccard similarity coefficient = 0.4). Consistent with the notion that cluster 1 cells displayed a stronger resemblance to normal epithelial cells, Gene Ontology (GO) enrichment analysis revealed that module 1 was dominated by genes associated with epidermal cell differentiation and development (Additional file 1: Fig. S2I). Notable genes within this category included EMP1, SCEL, SPINK5, ANXA1, SPRR3, PPL, CNFN, and TGM3. In contrast, gene modules associated with pathological complete response to NAT were primarily implicated in biological processes such as extracellular matrix organization (module 6), and activation of protein complex assembly (module 9) (Additional file 1: Fig. S2J). Regarding gene modules significantly contributing to MPR, the results of GO functional analysis unveiled a noteworthy enrichment of genes associated with epidermis development and response to type I interferon (module 7), response to hypoxia (module 10), and protein targeting to ER (module 4) (Additional file 1: Fig. S2K).

Given the unfortunate unavailability of survival data from the study cohort, we conducted an exploratory analysis utilizing the TCGA-ESCA dataset (N = 159). A Cox regression model was applied to investigate the relationship between cancer-associated gene modules and the patient’s overall survival (OS). Remarkably, higher expression of gene modules 6 and 12 was significantly associated with a more favorable prognosis of patients (P = 0.039 and 0.044, respectively; Additional file 1: Fig. S2L, M). In contrast, cluster 10-associated module 9 and cluster 5-associated module 10 were likely associated with a worse prognosis, though the difference did not reach statistical significance (Additional file 1: Fig. S2H, L). Nevertheless, it is worth noting that bulk RNA-seq data was used rather than scRNA-seq data and that the external dataset adheres to a different treatment regimen from the clinical settings in our analysis. Further exploration and validation of the association between cancer cell differentiation and immunotherapy-related patient prognosis are warranted in future studies.

We also employed SCENIC [49] to identify specific combinations of transcription factors that govern the expression of their respective target genes within cancer cells at the single-cell level. The top 5 regulons in cancer cells from patients exhibiting various responses to neoadjuvant anti-PD-1 combination therapy showed substantial differences in their expression patterns (Fig. 2H). Notably, pre-treatment pCR tumors (pCR_T_B) showed an upward trend in the expression of ZNF217 and MEF2A, whereas a downward trend was observed in the expression of TFCP2L1, ZSCAN31, and SOX2 compared to the other two groups (all P < 0.001, Jonckheere-Terpstra test) (Fig. 2I). Importantly, SOX2 has been extensively studied in the context of cancer, playing a pivotal role in tumor cell proliferation, cancer metastasis, and response to therapies [61]. Accumulating evidence suggests that SOX2 promotes cancer cell stemness, and its overexpression is associated with lymph node metastasis and poor prognosis in cancer patients [62]. Collectively, these findings suggest that the dynamic change in these transcription factors and their associated regulatory networks might be associated with the pathological response of patients who received NAT.

While PD-L1 expression has been employed to predict the efficacy of pembrolizumab in advanced EC [63], the predictive role of PD-L1 expression remains inconclusive for neoadjuvant therapy in EC patients. In our study, a total of 14 patients exhibited TPS scores < 50% or CPS scores < 50, thereby being classified as having PD-L1 low expression. Among these patients, 9 individuals (64.3%) exhibited pCR or MPR to NAT and five demonstrated IPR examined by conventional histologic approach (Additional file 3: Table S1). Dissecting the transcriptomes of cancer cells at single-cell resolution in low PD-L1 patients may provide significant clinical insights for understanding how these patients respond to neoadjuvant chemo-immunotherapy. In patients exhibiting low PD-L1 expression, a notable upregulation of genes associated with the protein targeting pathway was observed in tumors that achieved pCR. Conversely, in IPR tumors, there was an increased expression of genes associated with responding to oxidative stress (Additional file 1: Fig. S2N-P). These findings collectively suggest a plausible connection between the activation of these transcriptional signatures in cancer cells and the different pathological responses to NAT observed in patients with low PD-L1 expression.

Transcriptional analysis was performed to identify treatment response-related molecular signatures in CD8 + effector T cells at the single-cell level. Emphasis was placed on identifying DEGs and enriched pathways within pre-treatment ESCC tumors of pCR and IPR patients who exhibit the most pronounced variability in response to NAT (Fig. 3D). The expression levels of genes involved in the cellular response to IFNγ (such as MT2A, CCL4, and B2M) were notably upregulated in baseline pCR tumors compared to those in IPR tumors (Fig. 3E, F). Of note, a decreased expression of metallothionein (MT) family genes, including MT2A, MT1X, and MT1E, has been previously associated with resistance to anti-PD-1 therapy in ESCC [64]. While MT2A expression was significantly higher in MPR patients compared to those with IPR (P = 0.008), the difference did not reach statistical significance (Additional file 1: Fig. S3C-F). These findings suggest that the functional significance of CD8 + effector T cells in baseline MPR patients may resemble that in pCR patients rather than IPR patients. Assessing DEGs across various comparison groups, we found that upregulated genes linked to better pathological responses were mainly markers associated with cellular response to IFNγ (CCL4, MT2A, and B2M), neutrophil degranulation (B2M and HSPA8), regulation of I-kappaB kinase/NF-kappaB signaling (S100A4), and negative regulation of T cell apoptotic process (TSC22D3) (Fig. 3G). In contrast, downregulated gene signatures in baseline ESCC tumors of patients with better NAT-associated pathological response were primarily involved in the immune response-activating cell surface receptor signaling pathway (such as IGHG4, IGHA1, C3AR1, FCER1G, TYROBP, AAPL1, and LPXN) (Fig. 3E; Additional file 1: Fig. S3G). Among these genes, higher expression of C3AR1, FCER1G, and AAPL1 has been previously linked to an unfavorable prognosis in ESCC patients [65, 66].

To gain further insight into potential treatment strategies, we examined the top 30 ligand-receptor interactions between CD8 + effector T cells and cancer cells, as well as other T cell subsets. We observed a significantly lower number of cell–cell interactions in baseline tumors of IPR patients compared to those with pCR (Fig. 3H, I). In both IPR and pCR tumors, CCL5 was commonly expressed as a ligand on CD8 + effector T cells, with its receptor CCR1 expressed on NK cells and Treg cells. However, the identification of CCL5-CCR1 interactions within CD8 + effector T cells was observed only in baseline ESCC tumors of pCR patients. In contrast, CCL5-CCR5 interactions were widely expressed across various cell types in pCR tumors at baseline, except for cancer cells, NK cells, and innate lymphoid cells (ILC). Besides, CCL3 and CCL4 may also stimulate cell–cell communications through most immune cell subsets by binding to its receptor CCR5. Through assessment of the diverse expression patterns of immune checkpoints, we noticed a significant enrichment of interactions, including TIGIT-PVR and LGALS9-HAVCR2 in baseline IPR tumors, suggesting novel immunotherapy targets for ESCC patients who may not benefit from NAT (Fig. 3J, K).

Next, transcriptional changes in CD8 + effector T cells associated with pathological response to NAT were examined using pCR and IPR patients with low PD-L1 expression. Consistently, pCR patients with low PD-L1 expression exhibited higher expression levels of genes involved in the cellular response to IFNγ, along with lower expression of genes associated with immune response-activating signaling transduction pathway (Additional file: Fig. S3H-J). Overall, the recapitulation of GO enrichment analysis in low PD-L1 patients suggested that increased expression of IFNγ-responsive gene signatures and decreased expression of genes related to immune response-activating signaling pathway in CD8 + effector T cells might predict a favorable pathological response in locally advanced ESCC patients who underwent NAT. 1

We then investigated the interactions between NK cells and various cell types. The expression of LRPs involving NK cells and immune cells, such as CCL5-CCR4, CCL5-CCR1, and CCL4-CCR8, were exclusively observed in pCR patients before NAT but not in IPR patients (Fig. 4C, D). Furthermore, while the CCL5-CCR1 interaction was detected between CD8 + effector T cells and NK cells in both pre-treatment pCR patients and IPR patients, this interaction pair was also identified in pCR patients between CD8 + naïve T cells and NK cells, as well as within the NK cell population. A range of immune checkpoint pairs, such as TNFSF4-TNFRSF4, CTLA4-CD86, LGALS9-HAVCR2, and TIGIT-PVR, between NK cells and CD8 + effector/exhausted T cells, were significantly enriched in pre-treatment IPR patients, indicating additional potential targets for immunotherapy in nonresponders following the current treatment regimen (Fig. 4E, F).

Similar to what was found in CD8 + effector T cells, in low PD-L1 expression pCR patients, NK cells exhibited a markedly elevated expression of genes involved in the cellular response to IFNγ, as well as those implicated in the response to interleukin-1, such as CCL4, CCL3, CCL5, XCL1, UBB, ANXA1, and XCL2 (Additional file 1: Fig. S4D). Conversely, in IPR patients with low PD-L1 expression, there was a significant upregulation of genes associated with antiviral responses and those pivotal in RNA catabolic processes.

The examination of cell–cell interactions revealed that Treg cells displayed restricted interactions with other T cell subsets in IPR patients, whereas the communication was slightly intensified in pCR patients by the presence of CCL5-CCR1 and CCL4-CCR8 interactions (Fig. 5C, D). Immune checkpoint pairs, such as TIGIT-PVR, TNFSF4-TNFRSF4, LGALS9-HAVCR2, and CD70-CD27, were identified in baseline IPR patients as potential targets for immunotherapy (Fig. 5E, F).

Additional file 1. All supplementary figures (Fig. S1-Fig. S5).Additional file 2. scRNA-seq quality control data.Additional file 3:Table S1. Clinical characteristics of patients grouped by pathological response.Additional file 4. Cell count data.Additional file 5. Differential expression analysis data in cancer cells.Additional file 6. Differential expression analysis data in CD8+ effector T cells.Additional file 7. Differential expression analysis data in NK cells.Additional file 8. Differential expression analysis data in Treg cells.Additional file 9. Raw data used for CellphoneDB analysis.