Dissecting esophageal squamous-cell carcinoma ecosystem by single-cell transcriptomic analysis

Esophageal squamous-cell carcinoma (ESCC), one of the most prevalent and lethal malignant disease, has a complex but unknown tumor ecosystem. Here, we investigate the composition of ESCC tumors based on 208,659 single-cell transcriptomes derived from 60 individuals. We identify 8 common expression programs from malignant epithelial cells and discover 42 cell types, including 26 immune cell and 16 nonimmune stromal cell subtypes in the tumor microenvironment (TME), and analyse the interactions between cancer cells and other cells and the interactions among different cell types in the TME. Moreover, we link the cancer cell transcriptomes to the somatic mutations and identify several markers significantly associated with patients’ survival, which may be relevant to precision care of ESCC patients. These results reveal the immunosuppressive status in the ESCC TME and further our understanding of ESCC.

Esophageal cancer ranks the sixth leading cause of cancer-related death worldwide1 with esophageal squamous-cell carcinoma (ESCC) accounting for over 90% of total esophageal cancer in China. ESCC has poor prognosis with the 5-year survival rates around 20–30% probably due to the difficulty in early diagnosis and the lack of effective therapy1,2,3. Deciphering the complexity of ESCC tumor and its microenvironment is therefore fundamental for establishing early diagnosis and creating effective and precision treatment. In recent years, whole-genome/exome sequencing analysis on bulk tumor tissues has suggested that mutations of certain genes, such as TP53 and NOTCH1, or aberrant pathways, such as cell cycle and PI3K-AKT, may drive the development of ESCC4,5,6,7. However, it seems that genomic alterations are not completely attributed to the initiation and progression of ESCC and failed to translate into clinical use because such alterations are also seen in aging but normal epithelium8,9.

In addition, the tumor microenvironment (TME) compositions including various immune cells and nonimmune stromal cells also play determinant roles in tumor growth or decline10,11,12. It has been shown that immune cells and other stromal cells in TME rarely have genomic alterations13,14,15,16; however, they might interact with cancer cells to affect tumor progression and anticancer treatment17,18,19. For example, the checkpoint blockade immunotherapy can only benefit limited part of cancer patients, due to the cancer-immune interplays20,21,22. The failure of several anti-angiogenic therapies for cancer has been believed to be due to the metabolic adaptation and reprogramming of cancer cells and the abnormality of endothelial cells and their interaction with pericytes23.

Transcriptomic aberrancies are also critical events in promoting cancer development and progression24,25,26. In the last decades, many transcriptomic studies on ESCC have been reported, but previous studies were relied on bulk-tissue-based cDNA microarrays or RNA-sequencing analysis4,5,27,28,29. The bulk-tissue-based transcriptome analysis detects a mixed transcriptome derived from both cancer cells and TME components and usually generates massive marker genes in diversified pathways, which do not tell the dynamic status of various cells and their interactions during the development of cancer in the heterogeneous tumor tissues. In this regard, transcriptome analysis using bulk tumor tissues seems not suitable for the purpose of deciphering various unknown molecular events at the transcriptome level from different cells involving in the progression of ESCC.

High-throughput single-cell RNA sequencing (scRNA-seq), a valid method developed in recent years, has been proved to enable the dissection of heterogeneous tumors and deciphering the interaction between cancer cells and their microenvironment components30,31,32,33,34. Elucidating the transcriptome characteristics of cancer cells and the microenvironment components and their interactions, which are largely unknown in ESCC, is basic and fundamental in further understanding the cancer and developing effective early diagnosis and treatment strategies.

In the present study, we have performed scRNA-seq on a large scale of cells derived from ESCC tumors obtained from 60 patients to decode the transcriptome alterations associated with the cancer. We have also performed whole-exome/genome sequencing on bulk tissues from the same ESCC tumors and incorporated genomic data for analysis. We have dissected 8 expression programs from malignant epithelial cells and decomposed the TME compositions into 42 functional subtypes. By integrating these results, we have established a primary association framework between cancer cells and various noncancerous cells in the TME, which contributes to ESCC progression and prognosis.

Overview of ESCC ecosystem characterized by scRNA-seq

To decipher the cell composition within the ESCC tumors, we performed scRNA-seq and T cell receptor (TCR)-seq (for CD45+ cells only) on 60 ESCC tumor and 4 adjacent normal tissue samples obtained from 60 individuals using 10X Genomics platform (Fig. 1a; Supplementary Fig. 1a; Supplementary Data 1a and 2). After quality controls, we retained single-cell transcriptome of 208,659 cells including 97,631 CD45− and 111,028 CD45+ cells. Following regressing against read depth and mitochondrial read count, we performed graph-based clustering of the combined CD45− and CD45+ dataset and annotated the clusters using established marker genes (Fig. 1b, c; Supplementary Fig. 1b, c). We identified 8 main cell populations: epithelial cells (N = 44,730), fibroblasts (N = 37,213), endothelial cells (N = 11,267), pericytes (N = 3102), and fibroblastic reticular cells (FRC; N = 1,319) from CD45− dataset and T cells (N = 69,278), B cells (N = 22,477) and myeloid cells (N = 19,273) from CD45+ dataset. The epithelial cells exhibited high expression level of classic epithelial markers including EPCAM, SFN, and cytokeratins, and high genomic instability as demonstrated by severe copy number variations (CNVs) inferred from the transcriptome dataset and bulk whole-genome sequencing (WGS) results (Supplementary Fig. 1d–g), indicating that most epithelial cells are malignant. We then systematically analyzed the ecosystem compositions of each ESCC and found that the cell-type proportions were highly variable within and across patients, with some variations being associated with tissue type and tumor stage (Supplementary Fig. 1h). ESCC tumors had more epithelial cells and pericytes but fewer fibroblasts than adjacent normal tissues (P < 0.05; Fig. 1b; Supplementary Fig. 1i; Supplementary Data 3 and 4). The proportion of B cells was significantly less in stage II/III ESCC tumors compared with that in stage I tumors (P = 0.013; Supplementary Fig. 1j). These results suggest that ESCC ecosystems are highly heterogeneous, which should be further deciphered in the following analysis.

a Scheme of the overall study design. b, c tSNE plots of 97,631 CD45- cells (b) and 111,028 CD45+ cells (c), colored by cell type (left) and tissue type (top right). The proportions of cell types in normal and tumor tissues are shown on the bottom right (colors as in the left panel). Cell types that significantly increase (solid line filled with color) or decrease (dashed line) in tumor tissue (two-sided Wilcoxon test, P < 0.05) are indicated by links between normal and tumor bars. P values: 0.001 (epithelial cell), 0.001 (fibroblast), 0.010 (pericyte).

Full size image

Correlation between intra- and inter-tumor heterogeneity of epithelial cells

We then investigated whether and how the expression states varied among epithelial cells. By analyzing the transcriptome patterns of epithelial cells in 52 tumors that had >100 epithelial cells, we found that they could be divided into 38 clusters (Fig. 2a). We then analyzed the patient contributions to each cluster and found that 24 clusters had cells (≥75%) from sole individual patients while other 14 clusters had the expression patterns shared by multiple patients (Fig. 2b). We defined the 24 clusters as Group 1 clusters (with ≥75% personal cells) and the remaining as Group 2 clusters. We found Group 1 clusters showed increased pathway activities of cell proliferation and EMT, while Group 2 clusters presented activation of immunity-related pathways such as the complement, inflammatory, and IL2/STAT5 signaling (Supplementary Fig. 2a). We also compared the epithelial composition in each ESCC and found that ESCC from 21 patients were mostly composed of Group 1 clusters (≥60%) and ESCC from 31 patients were comprised of multiple Group 2 clusters (Fig. 2c). These results indicate a pervasive intra- and inter-tumor heterogeneity in ESCC.

We used principal components to measure the intra- and inter-tumor heterogeneity levels of epithelial cells. The intra-tumor heterogeneity is defined by the dispersion of cells from average within each sample, while inter-tumor heterogeneity is characterized by the distance of cells between each sample and global average as illustrated in Fig. 2d (see “Methods”). We integratedly quantified the levels of intra- and inter-tumor heterogeneity and found that they were positively correlated (r = 0.64, P = 7.3 × 10–7; Fig. 2e; Supplementary Fig. 2b, c). Among patients with the same intra-tumor heterogeneity levels, those having the Group 1 cluster had relatively higher inter-tumor heterogeneity levels than those without Group 1 cluster (P = 1.2 × 10–7; Fig. 2e, f; Supplementary Fig. 2d). The relative heterogeneity level was associated with age and drinking status of ESCC patients (Supplementary Fig. 2e). The relative heterogeneity level in ESCC with age >60 was significantly higher than that in ESCC with age ≤60 (P = 0.041), and it was also higher in ESCC from non-drinkers than that in ESCC from drinkers (P = 0.048).

Identification of eight common expression programs of epithelial cells in ESCC tumors

We developed a meta-cluster approach to uncover co-expressed programs among malignant cells from multiple samples. We first clustered the epithelial cells within each sample and generated a total of 274 intra-tumor subclusters in 52 ESCC tumors. We then defined an expression module for each subcluster with preferentially expressed 30 genes and subsequently aggregated the 274 modules into multiple recurrent expression programs using hierarchical clustering based on their expression profiles. We identified 8 expression programs with different functions and cell status (Fig. 3a, b; Supplementary Fig. 2f, g; Supplementary Data 5) and defined cells expressing ≥70% of genes in a given program as program cells (Supplementary Fig. 2h). We selected the most activated pathways in the program cells by comparison with non-program cells to analyze their functions (Fig. 3c). The mucosal immunity-like (Mucosal) program was characterized by the expression of genes associated with innate immune response (e.g., S100P) and mucosal defensive mechanisms including mucosal chemokine (e.g., CXCL17) and mucus production (e.g., AGR2and MUC20). Mucosal cells showed activation of the complement, TNFα signaling, apoptosis, and inflammatory response pathways. The stress responses (Stress) program consisted of immediate early genes (e.g., EGR1, JUN, and FOS) that are activated in response to widespread cellular stimuli and displayed upregulation of TNFα signaling, UV response, p53, and apoptosis pathways. The antigen presentation (AP) program had increased expression of major histocompatibility complex (MHC) class II molecules (e.g., CD74, HLA-DPA1, and HLA-DRA/B1/B5) that are involved in initiating adaptive antitumor immune responses. Immunity-related pathways, such as the allograft rejection, IFN-γ, IFN-α response, and the complement pathway, were activated in AP cells, probably indicating the reactivity to tumor neoantigens. The cell cycle (Cycling) program was characterized by high expression of genes involved in cell proliferation (e.g., CENPW, CKS1B, and BIRC5) and presented activation of the E2F targets, G2M checkpoint and MYC targets pathways, suggesting tumor cell proliferation. Two programs were associated with epithelial differentiation (Epi1/2). The Epi1 program was characterized by the expression of stress keratins (KRT6, KRT16, and KRT17) that are associated with keratinocyte hyperproliferation and therefore may play a role in enhancing tumorigenesis and tumor growth35,36,37. Epi1 cells displayed upregulation of pathways involved in cell growth and proliferation, such as mTORC1 signaling and PI3K/AKT/mTOR signaling pathways, as well as metabolic pathways including oxidative phosphorylation and glycolysis. The Epi2 program had the overexpressed genes related to the terminal differentiation such as envelope proteins (SPRR1A/1B) and calprotectin (S100A8/9), apical surface, the PI3K/AKT/mTOR signaling, the complement, and p53 pathways. The mesenchymal cell-like properties (Mes) program consisted of genes such as VIM and SPARC and showed activation of epithelial-mesenchymal transition (EMT) and angiogenesis pathways. Finally, the oxidative stress or detoxification (Oxd) program was characterized by the expression of multiple peroxidases and reductases (e.g., GPX2 and AKR1C1) involved in the defense against oxidative damage.

We performed the dependency analysis to infer the pairwise interactions among the expression programs and identified 5 significant co-occurring program pairs as well as 7 mutually exclusive program pairs with odds ratios ≥2 and ≤0.5, respectively (all P < 0.05) (Fig. 3d). The Epi1 was highly co-occurred with the Epi2, probably reflecting the need of both Epi1 and Epi2 for epithelial cell differentiation. The Epi1 program was also co-occurred with the Cycling program, probably reflecting the activation of cell growth and proliferation pathways in Epi1 cells. The Epi2 program was highly co-occurred with the Mucosal program, suggesting that mucosal immunity is present in well-differentiated epithelial cells. The Mucosal program was exclusive of the Cycling and Mes programs, while the Mes program was exclusive of the Epi1, Epi2, and Oxd programs. The AP program was exclusive of the Epi2 and Oxd programs. We further measured the activity of these programs in each ESCC by quantifying the program scores as proportion of corresponding program cells (Fig. 3e). Patients were then hierarchically clustered based on the program scores (Fig. 3f). We observed that most ESCC tumors contained the epithelial cells that expressed more than one program and that the Epi1 scores were highly variable among tumors with median score = 0.38 (ranging from 0 to 0.97). Patients with high Epi1 score usually had high Cycling score. The Epi2 and Mucosal programs had very small activities with the median scores being 0.03 among patients. Together, these results indicate that epithelial cells in ESCC tumors had 8 expression programs and the activity of each program was highly heterogeneous among tumors.

Characterization of immunosuppressive ESCC tumor microenvironment

Since T cells are the most abundant tumor-infiltrating lymphocytes (TILs) and highly heterogeneous in the TME, dissecting T cell populations would allow the accurate characterization of the immune status in tumors. We re-clustered 69,278 T cells and identified 9 distinct T cell subtypes and phenotypically related natural killer (NK)/natural killer T (NKT) cells (Fig. 4a; Supplementary Fig. 3a, b; Supplementary Data 6a). T cell subtypes included naïve T (TN) cells, T helper 17 (TH17) cells, follicular helper T (TFH1/2) cells, regulatory T (Treg) cells, memory T (TMEM-CD4/CD8) cells, effector T (TEFF) cells and exhausted T (TEX) cells. Compared with TFH2 cells, TFH1 cells were more canonical as they expressed CXCR5 that helps TFH cells localize to B cell follicles and promote the differentiation and maturation of B cells. We found that ESCC tumors were rich in Treg and TEX cells but poor in TN, TMEM, and TEFF cells compared with adjacent normal tissues, indicating an immunosuppressive status in the TME (Fig. 4a; Supplementary Fig. 3c; Supplementary Data 3 and 4). Stage II/III tumors had fewer TH17 and TFH1 cells but more TEXcells than stage I tumors (Supplementary Fig. 3d). We calculated the cytotoxicity score and exhaustion score to quantify properties of CD8+ T cells (Supplementary Fig. 3e, f; see “Methods”) and found that TEX cells had the highest cytotoxicity and exhaustion scores. Most of TEX cells were likely tumor-reactive T cells because of high levels of CD39 (ENTPD1) and CD103 (ITGAE) and very low level of KLRG1 (Supplementary Fig. 3g)38,39,40,41. We also generated a Treg score to quantify the activity of CD4+ T cells (Supplementary Fig. 3h). We found stage II/III ESCC tumors had significantly elevated exhaustion and Treg scores compared with stage I tumors, indicating that the immunosuppressive status in the ESCC TME is worse with tumor progression (Fig. 4b).

a tSNE plots of 69,278 T cells, colored by cell type (left) and tissue type (top right). The proportions of cell types in normal and tumor tissues are shown on the bottom right (colors as in the left panel). See Fig. 1b legend for the line indications. b Boxplots showing exhaustion score (left) and Treg score (right) of T cells are differentially distributed in patients at stage I (N = 16) and stage II/III (N = 44). Pvalues of two-sided Wilcoxon test are shown. Boxplots show the median (central line), the 25–75% interquartile range (IQR) (box limits), the ±1.5 times IQR (Tukey whiskers), and all data points, among which the lowest and the highest points indicate minimal and maximal values, respectively. c tSNE plot of T cells colored by clone size. d Bar plot showing fractions of proliferating cells in each T cell subtype. The cell numbers of T cell subtypes are labeled on the right. e tSNE plots of 19,273 myeloid cells, colored by cell type (left) and tissue type (top right). The proportions of cell types in normal and tumor tissues are shown on the bottom right (colors as in the left panel). See Fig. 1blegend for the line indications. f Violin plots showing the expression distribution of selected genes involved in immune suppression and angiogenesis in the monocyte and TAM clusters. g Violin plots showing the expression distribution of selected genes involved in DC maturation and immunosuppression in the three DC clusters. h Boxplots showing the expression levels of PD-L1 and PD-L2 in tDC are higher in tumor tissues (N = 60) compare with normal tissues (N = 4; two-sided Student’s t-test). The data presented as mean ± SEM. i, j 100,000 autologous CD8+ T cells isolated from human patients’ peripheral blood and tDCs were co-cultured in a 96-well plate at a ratio of 10:1 with Human T-Activator CD3/CD28. Anti-PD-L1 was administered at 100 μg/ml. T cell number was counted by hemocytometer at day 1 and day 2 (i). IL-2 and IFN-γ levels of the coculture supernatant at day 2 were measured by BD Cytometric Bead Array (CBA) kit (j). In (i, j), P values are derived from two-sided Wilcoxon test; *P < 0.05, **P < 0.01, ***P < 0.001. The data presented as mean ± SD and mean ± SEM in (i) and (j), respectively. Data in (i, j) are representative of three independent experiments for each group. k Clustered heatmap showing two-sided Spearman correlation coefficients between proportions of immune cell subtypes in tumor samples. P < 0.05 is used as the cutoff value.

Full size image

We profiled TCR repertoires to analyze the lineage structure of T cells in ESCC and found that the TCR clonotype compositions and proliferative cell proportions were highly diverse across different T cell subtypes (Fig. 4c, d; Supplementary Fig. 3i). The TEX cells had the highest level of clonal expansion (clone size ≥3 cells; Supplementary Fig. 3i) and the highest proportion of proliferation (nearly 10-fold greater than proliferative TEFF cells) (Fig. 4d). By measuring the proportion of clonotypes in a given T cell subtype (primary phenotype) shared with another T cell subtype (secondary phenotype), we found that CD8+ cell subtypes had a high degree of clonal sharing and extensive transitions existed among TEFF, TEX, and TMEM subtypes (Supplementary Fig. 3j). CD4+ T cell subtypes, except for TFH1 and TFH2 cells, had a low degree of clonal sharing and null sharing with Treg cells. These results indicate that TEX cells were highly proliferative and dynamic in the immunosuppressive TME throughout ESCC progression. Using the scTCR-seq and scRNA-seq data from 4 patients with matched tumor and adjacent normal samples, we found 23.9% of the T-cell clones in tumors were shared with normal samples, suggesting the parallel expansion of TEFF and TMEM-CD8 in these sites (Supplementary Data 7)42.

We also profiled B cells (N = 22,477) in both ESCC tumor and adjacent normal tissues and generated 5 subtypes designated as resting B, activated B, germinal center B (GC B1/B2) and plasma cells (Supplementary Fig. 4a, b; Supplementary Data 6b). High MKI67 and TOP2Aexpression in GC B2 cells suggested that these cells were active in proliferation (Supplementary Fig. 4b). The presence of GC B, FRCs and TFH cells indicates the existence of tertiary lymphoid structures (TLSs) which plays an antitumor role in the otherwise immunosuppressive TME43. These B cell subtypes provided key snapshots of B cells activation, proliferation, and differentiation in ESCC tumors.

Identification of tDC as a major player in the immunosuppressive ESCC microenvironment

Since tumor-infiltrating myeloid cells (TIMs) have been shown to be fundamental in regulating both innate and adaptive immune responses and facilitating tumor angiogenesis, invasion and metastasis44,45, we re-clustered all myeloid cells (N = 19,273) to look at their status and potential roles in ESCC. We found that TIMs in ESCC had 11 subtypes comprising 4 categories: monocytes (Mono01-03), tumor-associated macrophages (TAM01-04), mast cells (Mast) and dendritic cells (DC) (Fig. 4e; Supplementary Fig. 4c; Supplementary Data 6c). DC had three distinct subtypes with divergent roles: conventional DC (cDC), tolerogenic DC (tDC) and plasmacytoid DC (pDC). We found that monocytes and TAMs expressed not only genes associated with immunosuppression (e.g., TGFB1 and COX2), but also genes involved in angiogenesis (e.g., VEGFA, CXCL8, MMP9 and MMP12) (Fig. 4f). VEGFAwas upregulated in monocytes, while MMPs were mainly expressed in TAMs.

Mature cDC had high expression of MHC class II molecules (e.g., HLA-DRA), costimulatory factors (e.g., CD80 and CD86), and proinflammatory cytokines (e.g., TNF and IL1B) that may activate T cells and other antitumor immunity (Fig. 4g). Similar to mature cDC, tDC also had high expression of MHC class II molecules and costimulatory factors, but null expression of the proinflammatory cytokines TNF and IL1B, suggesting semi-mature phenotype of these cells. In addition, we found that tDC expressed the highest level of CCR7, which essentially contributes to both immunity and immunotolerance46,47,48; pDC specifically expressed ICOSLG, a gene involved in the expansion of Treg cells (Supplementary Fig. 4d). Furthermore, among all immune cell subtypes, tDC had the highest expression of the immune checkpoint genes (IDO1, PD-L1 and PD-L2) (Fig. 4g; Supplementary Fig. 4d–f). We also found that ESCC tumors were rich in tDC and the expression of PD-L1/L2 in tDC was significantly higher in tumors compared with adjacent normal tissues (Fig. 4h; Supplementary Fig. 4g; Supplementary Data 3 and 4). Ligand–receptor interaction analysis showed that tDC had stronger interactions with multiple subtypes of T cells than that of cDC, pDC, or TAM in terms of immunosuppression (Supplementary Fig. 4h). Coculture experiment of tDC isolated from ESCC tissue of one patient and autologous CD8+ T cells isolated from human peripheral blood sample with the stimulation by anti-CD3/28 showed that tDCs significantly inhibited CD8+ T cell proliferation (Fig. 4i). Furthermore, CD8+ T cells co-cultured with tDCs had diminished production of IL-2 and IFN-γ. We found that when anti-PD-L1 antibody was included in the coculture, the reduced CD8+ T cell proliferation and suppressed IL-2/IFN-γ production were rescued, suggesting that the suppressing effect of tDC on CD8+ T cell proliferation and activation is likely mediated through PD1 and PD-L1 interaction (Fig. 4i, j). Collectively, these results suggest that tDCs among TIMs play a crucial role in the immunosuppressive ESCC TME.

Identification of immune cell interwoven entities in the ESCC microenvironment

We then investigated the interactions and associations among different types of immune cells in tumor samples. Cells were clustered into multiple groups with different association patterns consistent with their related functions (Fig. 4k). We found that TN, TFH1, and B cells (activated, resting, and GC B cells) as well as TH17, TEFF, and TMEM-CD8 cells were aggregated into two major groups that contribute to immune responses. We also identified an immunosuppressive group consisting of Treg, TEX and plasma cells and found negative correlations between suppressive and effective immune cell subtypes. Furthermore, we found that the proportions of CXCR5+ TFH1 and GC B cells were highly correlated (r = 0.66, P = 2.4 × 10−9; Fig. 4k; Supplementary Fig. 4i), which is in line with the known fact that these two cell types are directly interactive in the germinal center49. Together, these results demonstrate a broad interaction among these immune cell subtypes that may play an important role in the construction of the immunosuppressive ESCC microenvironment.

Characterization of specific trajectories for fibroblast and pericyte differentiation

Fibroblasts and pericytes that can differentiate into myofibroblasts via pericyte-fibroblast transition are two important components of TME that may facilitate tumor invasion and metastasis23,50,51. Therefore, we profiled fibroblasts and pericytes to investigate the complexity and dynamic relationship of these cells in the ESCC TME. We found that they can be clustered into 9 subtypes as normal mucosa fibroblasts (NMF), normal activated fibroblasts (NAF1/2), cancer-associated fibroblasts (CAF1–4), pericytes and vascular smooth muscle cells (VSMC) (Fig. 5a; Supplementary Fig. 5a, b; Supplementary Data 6d). These fibroblast subtypes had distinct patterns of marker gene expression and pathway activities as shown in Fig. 5b and Supplementary Fig. 5c. NMF expressed a panel of genes encoding the protease inhibitors (e.g., SLPI and PI16) and genes involved in the complement, coagulation, peroxisome, and apoptosis pathways. NAF1/2 expressed the wound healing response-related genes (e.g., IGF1, C7, and APOD). CAF1 and CAF2 expressed proinflammatory chemokines (e.g., CXCL1 and CXCL6) and other cytokines involved in recruiting immune cells showing an activated inflammatory status. Interestingly, while CAF1 expressed high levels of CCL2/11 and CXCL14, CAF2 cells expressed activated-pericyte-specific cytokines CXCL5/8 and CSF352 and had the most significant upregulation of genes in inflammatory response pathway among all fibroblast subtypes. CAF3 showed an expression pattern resembling CAF4 (i.e., myofibroblasts) and expressed lower levels of myofibroblasts’ hallmark genes (e.g., TAGLN, ACTA2, and ACTG2) than CAF4. Both CAF3 and CAF4 had upregulation of apical junction and EMT pathways. Among all fibroblast subtypes, CAF4 had the highest activities of glycolysis and angiogenesis, two hallmarks of cancer. Furthermore, we found that each subtype had the distinct collagen repertoire; for example, NMF had lower expression of all collagen genes than CAFs, while CAF4 expressed the highest levels of COL1A1/2, COL3A1, COL5A1/2, and COL6A1/2/3 (Supplementary Fig. 5d).

a tSNE plots of 40,315 fibroblasts, colored by cell type (left) and tissue type (top right). The proportions of cell types in normal and tumor tissues are shown on the bottom right. See Fig. 1blegend for the line indications. b Dotplot showing expression status of selected fibroblast markers in each cell cluster. c Diffusion map of fibroblasts. Dot represents single-cell colored by cluster. The two main differentiation branches are indicated with solid arrows. d PAGA analysis of fibroblast subtypes. Line thickness corresponds to the level of connectivity. e Boxplots showing proportion of CAF3 (left) and CAF4 (right) in normal tissues (N = 4), stage I (N = 16) and stage II/III (N = 44) tumors (two-sided Wilcoxon test). Boxplots show the median (central line), the 25–75% interquartile range (IQR) (box limits), the ±1.5 times IQR (Tukey whiskers), and all data points, among which the lowest and the highest points indicate minimal and maximal values, respectively. f tSNE plots of 11,267 endothelial cells, colored by cell type (left) and tissue type (top right). The proportions of cell types in normal and tumor tissues are shown on the bottom right. See Fig. 1b legend for the line indications. g Heatmap showing differences in pathway activities scored per cell by GSVA between different endothelial clusters. h Violin plots showing the expression distribution of selected genes in the endothelial cell clusters. i Dotplot showing selected ligand–receptor interactions. The means of the average expression level of PDGFB in TECs and interacting receptor, PDGFRA and PDGFRB, in fibroblast subtypes are indicated by color. j Spearman’s correlation between the mean expression of PDGFB in TECs and the ratio of pericyte to endothelial cell in all samples. k Spearman’s correlations between the mean expression of PDGFB in TECs and the proportion of CAF2 (left) and CAF4 (right) in all samples, respectively. l Clustered heatmap showing two-sided Spearman correlation coefficients between proportions of immune and nonimmune cell subtypes in tumor samples. P < 0.05 is used as the cutoff value.

Full size image

Diffusion map and principal component analysis showed two distinct developmental trajectories of fibroblasts (Fig. 5c; Supplementary Fig. 5e; see “Methods”). In one branch, we detected a sequential continuum of phenotypes from NMF to NAF1/2, CAF1, CAF3 and CAF4, but in another branch, we found that pericytes, CAF2 and CAF4 formed a continuous and sequential path. We inferred the connectivity structures between the subtypes using partition-based graph abstraction (PAGA) analysis and found that CAF1 were connected to NAFs and CAF4, while CAF2 were the linkage between pericytes and CAF4, which was compatible with the results of diffusion map analysis (Fig. 5d). We also found that TGF-β pathway that is pivotal in the fibroblast-to-myofibroblast transition53 was activated in CAF1/2/3 (Supplementary Fig. 5c). These results suggest that both normal fibroblasts and pericytes may differentiate into myofibroblasts.

Interestingly, we found that the proportions of NMF and NAF1 were more prevalent in normal tissues than tumor tissues, while CAF2/3/4, VSMC, and pericytes were predominant in tumor tissues compared with normal tissues (Fig. 5a; Supplementary Fig. 5f; Supplementary Data 3 and 4). Moreover, we found that stage II/III ESCC tumors had more CAF3 and CAF4 but fewer CAF1 than stage I tumors (Fig. 5e; Supplementary Fig. 5g). These results indicate an accumulation of CAFs, especially the myofibroblasts, in the ESCC microenvironment that might play crucial roles in tumor progression.

Identification of high PDGFB level in TECs and its role in pericyte-myofibroblast transition

Since blood vessels and endothelial cells (EC) often exhibit aberrant phenotypes and functions in the TME, which might restrict them to respond to anti-angiogenic therapy, we therefore addressed these issues. We clustered EC (N = 11,267) into 6 subtypes including 3 normal EC (NEC1–3) and 3 tumor EC (TEC1–3) (Fig. 5f; Supplementary Data 6e). NEC1/2 expressed vein endothelial markers (e.g., CPE and ACKR1) while NEC3 expressed artery markers (e.g., SEMA3G and GJA5)54 (Supplementary Fig. 5h). While NEC1 were prevalent in normal tissues, TECs were exclusively present in ESCC tumors and TEC2/3 had increased expression of ANGPT2 that disrupts pericyte–EC interactions to enable angiogenesis23(Fig. 5f, Supplementary Fig. 5h, i; Supplementary Data 3 and 4). We found that TECs had specifically activated pathways involved in cell proliferation, angiogenesis, TGF-β signaling, EMT and energy metabolism (Fig. 5g). Notably, TECs had lower expression levels of genes associated with antigen presentation (e.g., CD74 and HLA-DQA1) and cell adhesion (e.g., ICAM1 and VCAM1) than NECs, which has been known to impede immune cell infiltration55, but had higher expression of molecules associated with angiogenesis (e.g., VEGFR1/2/3 and PDGFB) than NECs (Fig. 5h).

Since PDGFB-PDGFRB signaling pathway play a pivotal role in pericyte-myofibroblast transition and PDGFB is mainly secreted by endothelial cells51, we analyzed its role in the interaction of TECs and fibroblasts especially pericytes. We found a strong interaction between TECs and pericytes based on the ligand–receptor interaction analysis (Fig. 5i). The ratio of pericytes to ECs was positively correlated with the levels of PDGFB in TECs (r = 0.52, P = 8.6 × 10−6; Fig. 5j), in line with the recruitment role of PDGFB for pericytes. In addition, the expression levels of PDGFB in TECs were correlated with the proportions of CAF2 (r = 0.64, P = 1.2 × 10−8) and CAF4 (r = 0.46, P = 1.2 × 10−4; Fig. 5k). Furthermore, SCENIC analysis showed that several ETS family transcription factors (EHF, ELF3, and ETS2) downstream of MAPK were activated in CAF2 that can be induced by the PDGFB-PDGFRB pathway (Supplementary Fig. 5j)56. These results suggest that TECs in the ESCC TME may promote tumor progression by facilitating tumor angiogenesis via inducing pericytes to myofibroblasts transition.

Interactions among immune and nonimmune stromal cells in the ESCC microenvironment

We next explored the associations between immune and nonimmune stromal cells by clustering the related cell subtypes according to their association patterns, resulting in the benign and malignant groups (Fig. 5l; Supplementary Fig. 5k–m; Supplementary Data 3; see “Methods”). The benign group comprised the NMF, NAF1, CAF1, and NECs subtypes, while the malignant group comprised CAF2−4, TECs, and VSMC. We observed strong associations between benign group and effective immune cells, including resting B, GC B, TN, TFH1, cDC and activated B cells; however, the malignant group tended to correlate with suppressive immune cells, such as Treg, TEX, TAM, tDC, and pDC. Treg cell proportions were positively correlated with CAF2−4 and TEC1/2 proportions but negatively correlated with NMF, NAF1 and NEC2 proportions. TEX cell proportions were positively correlated with VSMC and CAF3 proportions. In addition, the proportions of TH17 cells that produce IL-17, a cytokine associated with pericyte activation52, were correlated with the proportions of CAF2 (r = 0.36, P = 0.004; Supplementary Fig. 5n), which is in line with the developmental trajectories of pericytes, suggesting that CAF2 may be activated pericytes. These extensive interactions among immune and nonimmune stromal cells might shape an immunosuppressive ESCC TME that promoted tumor progression. We also found that some clinical characteristics of patients such as sex, smoking status, and drinking status were significantly associated with the proportions of TME cell types (Supplementary Fig. 5o).

Associations of epithelial expression programs and genomic alterations with TME compositions

We next wanted to look at the interactions between malignant epithelial cells and other cells in the ESCC TME (Fig. 6a; see “Methods”). We found that ESCC tumors with higher Mucosal program scores had more TN, TFH1, NK/NKT, GC B, and cDC but fewer TEX, plasma, tDC, and pDC. Such ESCC tumors also possessed more NEC2 but fewer CAF2−4. ESCC tumors with higher AP program scores exhibited more TEFF cells, resting B cells and NEC2 but fewer CAF4 and TEC2. However, ESCC tumors with higher Cycling program scores had fewer TMEM-CD8 cells, NMF, NAF1/2, CAF1, and NEC1/2 but more CAF2/3 and TEC2/3 as well as higher cytotoxicity, exhaustion, and Treg scores. ESCC tumors with high Epi1 scores had reduced infiltrating effective immune cells (TFH1 and GC B) and benign stromal cells (NMF, NAF1/2, and NEC1/2). Conversely, such ESCC tumors contained more malignant CAF3/4 and TEC1 and had higher cytotoxicity and Treg scores.

We next sought to know whether the genomic alterations may affect the ESCC ecosystem. We detected 11 putative driver genes mutated in ESCC by whole-exome sequencing (WES) and among them the most frequently mutated genes were TP53 (89.1%) and NOTCH1 (28.3%) (Supplementary Fig. 6a). By analyzing the differences in epithelial expression program scores and TME cell-type proportions associated with TP53 or NOTCH1 mutation status, we found that ESCC tumors with the mutant TP53 had higher Cycling and Epi1 program scores and higher immune suppression levels than tumors with the wild-type TP53 (Fig. 6b; Supplementary Fig. 6b). The TP53 mutations were also associated with higher cytotoxicity and Treg score and increased TAMs but decreased TN cell and monocyte infiltration. In contrast, ESCC tumors with NOTCH1 mutations had lower Epi1 and Epi2 program scores but higher Mes scores, reflecting the role of mutant NOTCH1 signaling in epithelial differentiation (Fig. 6b; Supplementary Fig. 6c). Tumors with NOTCH1 mutations also had increased TEFF cell infiltration than tumors without NOTCH1 mutations, indicating that the mutations might have produced neoantigens. By integrative analyzing the mutations of TP53 and NOTCH1, we found sole mutation of TP53 was associated with a reduced proportion of TN, TFH1, GC B, resting B, monocyte, and NEC2 and increased proportion of plasma, TAM and TEC1 as well as higher cytotoxicity and Treg score (Supplementary Fig. 6d). For tumors with TP53 mutation, the mutation of NOTCH1 presented higher infiltration of TEFF cells but reduced Epi1 and Epi2 scores than wild-type NOTCH1. The mutational signatures were also associated with the heterogeneity level of malignant cells and the TME compositions (Supplementary Fig. 6e, f). These results suggest that somatic mutations of some driver genes in esophageal epithelial cells may drive the alterations of epithelial expression programs and the TME compositions that facilitate tumor progression.

Correlation of gene expression levels in the Mucosal program with ESCC survival in patients

Since epithelial expression programs may reflect the genomic alteration information and drive remodeling of ESCC TME compositions, we therefore analyzed the association between epithelial expression programs and the disease outcome. We found that the expression levels of the identified 14 genes were significantly associated with survival time in ESCC patients (Discovery cohort, N = 139; all P log-rank < 0.05; Fig. 6c, Supplementary Data 1b and 8). Among the 14 associated genes, the AGR2, CXCL17 and MUC20 in the Mucosal program were also significantly associated with survival time in a validation cohort of 94 ESCC patients (Validation cohort 1; all P log-rank < 0.05; Supplementary Fig. 7a, b; Supplementary Data 1c and 8). AGR2 encodes an endoplasmic reticulum protein essential for intestinal mucus production that mediates fundamental innate immunity57. CXCL17produces a mucosa-associated chemokine that supports innate immunity and the homeostasis and integrity of the mucosa58. MUC20 encodes a mucin protein that is part of an insoluble mucous barrier59. These three genes are specifically expressed in epithelial cells so that can be measured by bulk-tissue RNA-seq (Supplementary Fig. 7c). We further categorized the expression level of these 3 genes as high (>median) or low (≤median) and integrated them as the Mucosal Index (MI) for further analysis. We found that in Discovery cohort and Validation cohort 1, patients with high MI (MI ≥ 2) had significantly longer median survival time than patients with medium MI (MI = 1) (P = 0.057 and P = 0.039, respectively) or low MI (MI = 0) (P = 3.7 × 10−4 and P = 0.001, respectively) (Fig. 6d, e). The correlations of the expression levels of 3 genes and ESCC survival time were further validated in another cohort consisting of 226 patients by immunohistochemical staining of protein levels in ESCC tumors (Validation cohort 2; Fig. 6f; Supplementary Data 1d) and the results also showed that higher expression levels of their proteins (expressed as H-score) were correlated with significantly longer ESCC survival time (all P log-rank < 0.05; Fig. 6g; Supplementary Fig. 7d–f). We found that the proportions of epithelial cells expressing AGR2/CXCL17/MUC20 were correlated with the Mucosal program scores in scRNA-seq data (r = 0.77, P = 2.4 × 10−11; Fig. 6h), suggesting that these three genes could be collectively used to predict the Mucosal program activity.

We further investigated Mucosal program-correlated epithelial aberrations in bulk-tissue RNA-seq data of discovery cohort and validation cohort 1 using GSVA (Fig. 6i). Among hallmark pathways differentially expressed between tumors with the high or low expression level of Mucosal genes, the EMT and angiogenesis pathways were significantly downregulated in Mucosal-high tumors in both cohorts, suggesting suppressed migration and proangiogenic abilities of tumors with high Mucosal activities. Tumor suppressor p53-related pathway was significantly enriched in Mucosal-high tumors in both cohorts, suggesting that Mucosal-high tumors are relatively less malignant. We also found the E2F targets and G2M checkpoint pathways were downregulated in Discovery cohort, suggesting suppressed cell proliferation of tumors with high Mucosal activity. These results collectively point to a relative less malignant phenotype and better TME of Mucosal-high ESCC (Fig. 6a), which might explain why patients with mucosal-high ESCC had a better prognosis.

In the present study, we have performed scRNA-seq in 208,659 cells from ESCC tumor and their adjacent normal samples collected from 60 patients and deciphered the phenotypes and compositions of the ESCC ecosystem. We found that samples can be reliably classified into two groups according to the different adoptions of shared clusters for epithelial cells in regard to the sample compositions and these differences can be further quantified using PCA-based heterogeneity scores. We have dissected 8 common expression programs from malignant epithelial cells and decomposed the TME compositions into 42 functional subtypes including 26 immune cell subtypes and 16 nonimmune stromal cell subtypes. We have also elucidated the possible interactions between cancer cells and TME cells and the interactions among different cell types in the TME. Moreover, we have dissected the relationships between ESCC ecosystem components and genomic alterations or disease outcomes. These results provide fundamental information for deeply understanding the mechanism for pathogenesis and progression of ESCC, which in turn may benefit the development of approaches for precision cares of ESCC patients.

The present study has provided several significant findings. Firstly, we have identified an immunosuppression status in ESCC TME featured by a large amount of TEX, Treg, and myeloid cell infiltration. Specially, we have revealed that in the TME, tDC express the highest levels of immune checkpoint genes (PD-L1/2 and IDO1) that have been shown to induce the T cell anergy and produce Treg cells60,61. Our experimental results showed that tDC can suppress the activation of CD8+ T cells, suggesting that they may play an important role in evoking immunosuppression. Given the central role of dendritic cells in controlling immunity homeostasis, measurement of tDC in ESCC tumors might predict whether the tumor is vulnerable to immunotherapy and targeting both tDC and T cells may be necessary in immunotherapy of ESCC, as reported in prostate cancer for a synergistic effect62. Furthermore, we did not identify the existence of recently reported MR1-restricted T cells, which have been reported to recognize a broad range of cancer cells63. This negative result might further support that ESCC tumors have the immunosuppressive TME.

Secondly, we have deciphered the complicated compositions of nonimmune stromal cells and their dynamic transition in the ESCC TME. We have discovered two hidden intermediate phenotypes of fibroblasts (CAF1 and CAF2) jointly manifesting two developmental trajectories among them. These transitionary fibroblast subtypes are activated and produce various proinflammatory cytokines that play important roles in ESCC progression. We have also found that the characteristic high expression of PDGFB in TECs is related to the proportions of pericytes, CAF2, and myofibroblasts in the ESCC TME, likely indicating that TECs play a part in reshaping the TME. These findings imply that targeting PDGFB signaling pathway may inhibit angiogenesis via blocking TECs to recruit pericytes and repressing the pericyte-myofibroblast transition51,64,65, which offers a new therapeutic option for cancer treatment. One limitation is that the differentiation trajectories for fibroblasts and pericytes and the role of PDGFB in ESCC TME need to be verified by further functional studies.

Thirdly, we have identified several different cellular interactions in the TME, such as the interaction between TFH1 and GC B cells, which produces an antitumor immune environment by the formation of TLS, the interaction between TEC2 and TN cells, which impedes infiltration of immune cells to tumor tissues, and the interaction between TH17 cells and CAF2 mediated by cytokines secreted from TH17 cells, which enhances CAF2 cell activity. CAF2 are the intermediate activated phenotype in the pericyte-myofibroblast transition and TH17 cells are implicated in activation of pericytes52, suggesting that TH17 cells may play an important role in pericyte-myofibroblast transition. We have found that the amount of cell subtypes that are phenotypically or functionally relevant is highly correlated and these cells can be grouped as malignant nonimmune stromal cells (e.g., CAFs and TECs) and suppressive immune cells (e.g., Treg, TEX, and TAM). These results suggest that cells in the ESCC TME are interwoven together for promoting tumor survival and development. Interestingly, we have found an extensive cooperation between active cancerous epithelial cells, exhibited by the expression programs, and malignant nonimmune cells or suppressive immune cells. The proliferative tumors had high amount of CAFs and TECs in their TME, providing aberrant extracellular matrix, cytokines and vascular niche for tumor cells to be progressive and invasive19. These results present a high-resolution picture of entire and detailed cellular interactions in the ESCC tumor ecosystem, which would be significant in developing new therapeutic strategies.

Furthermore, we have also linked the single-cell transcriptomes to the genomic alterations detected by bulk-tissue WES. The results have shown that the mutations in NOTCH1 are associated with decreased expression of the Epi1/2 programs but elevated expression of the Mes program in epithelial cells, indicating that the mutations may impair the differentiation of epithelial cells, a well-known hallmark of cancer. However, ESCC having the NOTCH1mutations seem to have more infiltrated TEFF cells, suggesting that tumors with the mutations have more effective immunity. In contrast, the TP53 mutations confer ESCC to be more proliferative and immunosuppressive and have an effect on the expression of the Epi1 program opposite to NOTCH1 mutations. It has been shown that the NOTCH1 mutations are less frequent in ESCC tumors than in aged normal esophageal tissues8,9. This phenomenon can be explained by scRNA-seq results showing that tumors with the TP53 mutations have higher transcriptomic expressions related to proliferation and immunosuppression than tumors with the NOTCH1 mutations. The proliferative advantage of cancer cells may result in accumulation of TP53 mutations in ESCC, whereas the differentiation-related NOTCH1mutations are accumulated in aged noncancerous esophageal epithelium.

Most importantly, we have found the expression levels of the mucosal immunity-like (Mucosal) program are significantly associated with ESCC survival time in patients. Because the high expression levels of the Mucosal program in tumor are strongly associated with the amount of infiltrating effective immune compositions such as TFH1, GC B, and cDC, one may expect that patients with high expression of the Mucosal program would have higher antitumor immunity and thus have better prognosis. Notably, GC B cells are the major cell type in TLS, which have been shown to be able to promote immunotherapy response in several types of malignancies66,67,68. We have identified that within the program the expression levels of CXCL17, AGR2 and MUC20 are the best markers associated with ESCC survival. Although the underlying mechanisms for AGR2 and MUC20 are currently not evident and warrant further investigation, CXCL17 has been reported as an antitumor factor since it can recruit dendritic cells into tumors, which may enhance antitumor immunity58. These results may be clinically relevant in precision and individualized care of ESCC patients. Furthermore, the production of genes in the Mucosal program, such as CXCL17, might be useful in treatment of cold ESCC tumors.

In summary, the present study provides an atlas of the ESCC ecosystem based on the large-scale single-cell RNA-seq results. We have depicted a comprehensive picture of interactive network among cancerous epithelial cells, stromal cells, and various infiltrating immune cells in the ESCC tumors and their TME. We have also linked the transcriptome in cancerous epithelial cells to their somatic genome alterations and identified several significant markers associated with patients’ survival. These results deepen and extend our understanding of the complexity of ESCC tumors.