Integrated single-cell transcriptome analysis reveals heterogeneity of esophageal squamous cell carcinoma microenvironment

Abstract

The tumor microenvironment is a highly complex ecosystem of diverse cell types, which shape cancer biology and impact the responsiveness to therapy. Here, we analyze the microenvironment of esophageal squamous cell carcinoma (ESCC) using single-cell transcriptome sequencing in 62,161 cells from blood, adjacent nonmalignant and matched tumor samples from 11 ESCC patients. We uncover heterogeneity in most cell types of the ESCC stroma, particularly in the fibroblast and immune cell compartments. We identify a tumor-specific subset of CST1⁺ myofibroblasts with prognostic values and potential biological significance. CST1⁺ myofibroblasts are also highly tumor-specific in other cancer types. Additionally, a subset of antigen-presenting fibroblasts is revealed and validated. Analyses of myeloid and T lymphoid lineages highlight the immunosuppressive nature of the ESCC microenvironment, and identify cancer-specific expression of immune checkpoint inhibitors. This work establishes a rich resource of stromal cell types of the ESCC microenvironment for further understanding of ESCC biology.

Fig. 1. Single-cell transcriptomic landscape of esophageal squamous cell carcinoma (ESCC).

A A schematic graph showing the study design. B UMAP (Uniform Manifold Approximation and Projection) visualization of the clustering of 41,237 cells from all 22 nonmalignant and tumor samples, color coded by either major cell type (left), sample type (middle) or patient origin (right). C Overlay of expression of representative marker genes for each cell type defined in (B). D The frequency of each cell type in nonmalignant and tumor samples (left), and in each of the 11 patients (right, an analysis restricted within tumor samples). Source data are provided as a Source Data file.

Fig. 2. Fibroblast heterogeneity in ESCC.

A UMAP visualization of the clustering of 12,126 fibroblast cells from all 22 nonmalignant and tumor samples. B The fraction of each fibroblast subset in nonmalignant and tumor samples (left), and in each of the 11 patients (right, an analysis restricted within tumor samples). C Dotplot showing the expression of top 10 most variable genes across each fibroblast subset. D Violin plots of the expression of representative ECM (extracellular matrix) genes and ACTA2 in each fibroblast subset. E Enrichment of GO (gene ontology) terms of each fibroblast subset (FDR-adjusted P < 0.001, Fisher exact test with multiple comparisons using topGO). Source data are provided as a Source Data file.

Fig. 3. CST1⁺ myofibroblasts are cancer-specific and predict survival outcome in ESCC.

A Box plots of the average levels of gene signature defined for F_3 in the TCGA bulk RNA-Seq and independent microarray datasets of ESCC samples. n = 60 paired tumor and nonmalignant samples in GSE53622, n = 119 paired tumor and nonmalignant samples in GSE53624, n = 80 tumor and 11 nonmalignant samples in TCGA. P values are calculated by two-tailed t test. B The expression of CST1 mRNA across myofibroblast (F_3), other fibroblasts and all the other cell types identified in ESCC. C Rank order of DE (differentially expressed) genes based on average logFC between tumor and nonmalignant samples within myofibroblasts; CST1 was identified as the most upregulated gene in tumor samples. N, nonmalignant; T, tumor. D Representative images of immunofluorescence double staining of both CST1 and COL1A1 in ESCC tumor and nonmalignant samples. Scale bar = 100 μm. E Quantification of the ratio of CST1⁺ cells out of COL1A1⁺ fibroblasts. F Representative images of IHC (immunohistochemistry) staining of CST1 in ESCC tumor and nonmalignant samples. Scale bar = 50 μm. G Quantification of the percent of CST1⁺ fibroblasts out of all stromal cells from IHC staining. H Kaplan-Meier curves of either overall survival or (I) disease-free survival of ESCC patients stratified by the abundance of CST1⁺ fibroblasts. J Kaplan-Meier curve of overall survival of ESCC patients stratified by the mRNA level of CST1 in an independent cohort. K Bar plots of the percentages of cells expressing CST1 mRNA in different scRNA-Seq datasets for different cancer types: ESCC (this study), lung, colon and head and neck. N.A., no data available from nonmalignant samples. In the box plots (A, E, G), the middle bar represents the median, and the box represents the interquartile range; whiskers indicate the maximum and minimum values. Dots are all the data points including outliers. P values are calculated by two-tailed Mann Whitney U test. Source data are provided as a Source Data file.

Fig. 4. Antigen-presentation fibroblast in ESCC.

A Violin plots of the expression of MHC II class genes and selected fibroblast marker genes in each fibroblast subset. B Representative images of immunofluorescence double staining of both HLA-DR and VIM in ESCC tumor and nonmalignant samples. Scale bar = 100 μm. C Quantification of the ratio of HLA-DR⁺ cells out of VIM⁺ fibroblasts. The middle bar represents the median, and the box represents the interquartile range; whiskers indicate the maximum and minimum values. Source data are provided as a Source Data file.

Fig. 5. Myeloid and T cell landscapes of ESCC.

A UMAP visualization of 17,923 immune cells extracted in silico from all 22 nonmalignant and tumor samples. B Dotplot showing the average expression of representative marker genes and percentages of expressed cells in each immune subset. C Violin plots showing the expression of selected genes for CD4⁺ T cell subsets or (D) CD8⁺ T cell subsets across nonmalignant (green) and tumor (blue) samples. E Volcano plot of the differentially expressed genes between the two CD8⁺ T cell subsets (CD8_1 and CD8_2 specific genes are highlighted in blue and red, respectively). F A circos plot showing the higher overall number of significant interacting pairs estimated by CellPhoneDB (P < 0.05) between myeloid and T cell subsets in tumor (blue) and nonmalignant (green) samples. G Estimated mean interaction scores for specific interactions (PDCD1-CD274, CTLA4-CD80/CD86) from indicated cell types in tumor and nonmalignant samples. N, nonmalignant; T, tumor. Source data are provided as a Source Data file.

Fig. 6. The cellular heterogeneity of myeloid compartment in ESCC.

A Violin plots showing expression of representative markers in either DC (dendritic cell) or (B) monocyte and macrophage subsets. C Dotplot showing the expression of top 10 differentially expressed genes across monocyte and macrophage subsets. D Scatter plot showing the correlation of M1 and M2 gene signatures in individual macrophage subsets using all myeloid cells as background. E The expression ratio of M2 over M1 gene signatures. F Histogram of the distribution of myeloid cells along the Diffusion-1 component from either the PBMC (peripheral blood mononuclear cells), tumor or nonmalignant samples. G Diffusion component analysis of the myeloid compartment displayed in 3D plot showing the first 3 diffusion components colored by either sample origin or (H) annotated subsets. The 1st diffusion component reflected a trajectory from blood monocytes to tissue monocytes/macrophages, the 2nd diffusion component reflected the activation of monocyte-derived macrophages in tissue, and the 3rd diffusion component reflected a trajectory from pDCs (plasmacytoids dendritic cell) to more mature DC subsets. I–K 3D plots showing the correlation of representative genes with each diffusion component. Source data are provided as a Source Data file.

Fig. 7. Validation of candidate cell subsets by multiplexed immunofluorescence (IF) staining.

A–B Representative images of multiplexed IF staining of tumor samples using Panel-1 (A) and Panel-2 (B). Scale bar = 100 μm. C–F Quantification of CST1⁺ fibroblasts (C), exhaustive CD8⁺ T cells (D), Treg cells (E), M2-like macrophages (F) in nonmalignant (green) and tumor (blue) samples. The number of regions quantified are shown in the parentheses. Data are presented as mean values ± SEM. P values are calculated by two-tailed Mann Whitney U test. G–H Scatter plots showing the positive correlation (G) between exhaustive CD8⁺ T cells and M2-like macrophages (n = 120), and (H) between CST1⁺ fibroblasts and Treg cells (n = 120). Source data are provided as a Source Data file.