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

Abstract

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.

Fig. 1. Overview of ESCC ecosystem characterized by scRNA-seq.

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).

Fig. 2. Correlation between intra- and inter-tumor heterogeneity of epithelial cells.

a tSNE plot of 44,122 epithelial cells, colored by cluster. Only tumors with at least 100 epithelial cells are shown. b Pie charts showing the proportions of cells originating from each patient detected in each cluster, colored by patient. Clusters are ordered by the proportion of the most dominant patient and are classified into Group 1 cluster or Group 2 cluster. c Stacked histograms showing proportions of epithelial cell clusters for each patient, colored by cluster as in (a). d Schematic showing the methods for quantifying the intra-tumor (solid line) and inter-tumor (dashed line) heterogeneities of epithelial cells for each patient based on the patient average (yellow star) and global average (red star). e Scatterplot showing the positive correlation between the intra- and inter-tumor heterogeneities. Two-sided Spearman correlation coefficient and P value are indicated. Shaded region indicates 95% confidence interval for the correlation. f Boxplot showing the relative inter-tumor heterogeneity levels for patients with Group 1 cluster (green) or without Group 1 cluster (yellow). Boxplots show the median (central line), the 25–75% interquartile range (IQR) (box limits), the ±1.5 times IQR (Tukey whiskers).

Fig. 3. Identification of eight common expression programs of epithelial cells in ESCC tumors.

a Heatmap showing pairwise correlations of 274 modules derived from 52 tumors. The common expression programs across tumors are aggregated into clusters. b Heatmap showing expression of genes within each program across single cells. Randomly selected 500 cells for each program are shown. Colors above columns correspond to cell state (see Methods). c Heatmap showing differences in pathway activities between program cells and non-program cells for each program scored by GSVA. Each column is normalized by z-score to indicate the relative pathway activities. d Heatmap showing odds ratios assessing for each pair of programs (rows, columns) if they are co-occurrent (≥2, red) or exclusive (≤0.5, blue) than expected by chance (P < 0.05). P values are derived from pairwise Fisher’s exact test. e Boxplot showing the proportion of program cells (hereafter refers to as program score) for each tumor (N = 52) among 8 expression programs, sorted by the median program score. 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 Clustered heatmap showing the normalized program scores for all programs in each tumor.

Fig. 4. Characterization of immunosuppressive ESCC tumor microenvironment.

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). P values 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. 1b legend 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.

Fig. 5. Characterization of specific trajectories for fibroblast and pericyte differentiation.

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. 1b legend 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.

Fig. 6. Correlation of gene expression levels in the Mucosal program with ESCC survival.

a Heatmap showing two-sided Spearman correlation coefficients between program scores and proportions of TME compositions as well as T cell function associated scores (hereafter refers to as T cell scores) in tumor samples. P < 0.05 is used as the cutoff value. b Bubble plot showing relationships between mutation status of WES-derived frequently mutated genes and program scores, proportions of TME compositions as well as T cell scores. Compositions that significantly increase or decrease in mutated tumors (two-sided Wilcoxon test, P < 0.05) are indicated by purple or green, respectively. Circle size represents P value. MUT, mutant; WT, wild type. c Associations between program gene expression levels and ESCC survival in Discovery cohort (N = 139) and Validation cohort 1 (N = 94). The significances of genes in two cohorts are indicated by color. P_log-rank < 0.05 (two-sided) is used as the cutoff value. d Kaplan–Meier (KM) curves of Discovery cohort (N = 139) grouped by the mucosal index (MI). e KM curves of Validation cohort 1 (N = 94) grouped by the MI. f Representative images of IHC of the proteins produced by the 3 Mucosal genes in Validation cohort 2. Shown are samples with MI 1 (12-9_A4 and 12-8_B4) and MI 2 (12-4_C4). Scale bar, 50 μm. g KM curves of Validation cohort 2 (N = 226) grouped by the MI. h The relationships between the proportion of epithelial cells expressed AGR2, CXCL17 or MUC20 (total expression >0) and Mucosal score in tumor samples. Two-sided Spearman correlation coefficient and P value are indicated. i Differences in selected pathway activities by GSVA between tumors with high or low expression level of Mucosal genes in Discovery cohort and Validation cohort 1, respectively. Shown are t values from the linear model. P values in (d, e, g) are derived from log-rank tests.