A scRNA-Seq census of human intramucosal ESCC identifies 12 distinct cell populations and differentiation relations
To create a comprehensive single-cell landscape of human intramucosal ESCC, we applied scRNA-seq using paired samples of endoscopically resected intramucosal ESCC, para-ESCC esophageal tissue from endoscopically resected specimen, and PBMC. Intramucosal ESCC lesion and para-ESCC esophageal tissue were identified using iodine staining on the endoscopically resected specimen. Namely, the non-staining area indicated early ESCC whereas deep-staining area represented non-ESCC paratumor tissue, respectively. A 1mm*1mm sample from the lesion and a 1mm*1mm sample from the paratumor area was taken from each resected specimen in order to minimize the influence on pathological assessment (Fig. 1A). Five pairs of samples were included after pathological confirmation, and demographic data indicated in Table 1. After quality control to remove low-quality cells, a total of 16,4715 cells were retained for further analyses. Among them, 48,711 cells derived from intramucosal ESCC samples (6736 ~ 11794, median: 9912/patient), 55,191 cells derived from paratumor samples (8417 ~ 13876, median: 10127/patient), and 60,813 cells derived from PBMC (10422 ~ 14256, median: 11951/patient). A median of 1634 genes per cell was detected.
Table 1
Demographic and pathological characteristics of patients.
| Age | Sex | Tissue size(cm) | Iodine uncolored area(cm) | Microscopic lesion(cm) | Stage (depth of tumor invasion) |
Patient01 | 55 | Male | 3*2.5 | 2.5*1.5 | 3*2 | T1a (Muscularis mucosae) |
Patient02 | 65 | Female | 3.8*2.3 | 2.5*1.3 | 3*l.2 | T1a (Lamina propria mucosa) |
Patient03 | 72 | Male | 3.9*5 | 2.7*2 | 3*2.5 | Tis |
Patient04 | 75 | Male | 9*8 | 8*7 | 4.8*4.3 | T1a (Lamina propria mucosa) |
Patient05 | 68 | Male | 5.8*2.8 | 5.2*2.2 | 4.8*1.5 | T1a (Lamina propria mucosa) |
Based on the expressions of canonical markers, we classified cells into 12 clusers by using t-distributed stochastic neighbor embedding (tSNE) (Fig. 1B), including megakaryocyte cells, cancer-associated fibroblasts (CAFs), dendritic cells, plasma cells, endothelial cells, myofibroblasts, proliferating cells, B cells, T cells, monocytes, epithelial cells, and natural killer (NK) cells. Marker genes of the identified 12 major cell types were cell type specific and consistent with classic, well-established markers for each respective cell population. (Fig. 1C, Supplementary Figure S1A).
Noteably, the proportion of immune cells in the tumor tissue increased significantly compared with the paratumor tissue, such as B and T cells. In contrast, epithelial cells are the dominant cluster in the paratumor tissue (Fig. 1D). What’s more, an increase in B cell proportion was observed in both the tumor and paratumor tissues from patient 3 (Fig. 1E), who was identified as high-grade intraepithelial neoplasia rather than carcinoma in situ (Table 1). In intramucosal ESCC, the average proportion of infiltrating B cells was significantly higher than that in the progressive ESCC8 (Fig. 1F). However, a decrease in the proportion of infiltrating B cells in progressive ESCC compared to other cancer types was reported in previous studies19, indicating the positive effects of tumor-infiltrating B cells.
Transcriptomic inter-tumor heterogeneity and transitional cell status of epithelial cells in human intramucosal ESCC
To investigate how normal esophageal epithelium develops into carcinoma, various tools were introduced to assess the genetic alterations and functional changes in epithelial cells harvested from tumor and paratumor tissues. A total of 28397 epithelial cells were identified and were classified into 7 subclusters designated as E1 ~ E7, and the tSNE plot revealed the subcluster distribution among sample pathology (Fig. 2A). E6 ~ E7 was enriched in tumor samples, while a higher fraction of E1, E2 and E3 clusters was observed in paratumor samples (Fig. 2B,C). GO terms were enriched in normal epithelial function such as epidemis development and cornification in E1, E2 and E3 subclusters. E4 was enriched in ATP metabolic process, and E5 ~ E7 were enriched in protein folding and RNA splicing (Fig. 2D). Although there was no distinct difference in gene expression profile across cell subclusters, a gradual increase in the expression level of oncogenes with a gradual decrease in tumor supressors could be identified. For example, TMPRSS11B, MUC21, CRNN were highly expressed in the upper right corner of the tSNE plot characterized by E1 subcluster, while the expression of SOX4, MDK and IFI6 was higher in the upper left corner characterized by E6 and E7 subclusters (Fig. 2E). DSC2, RT4 and IL1RN were highly expressed in subclusters other than E6 or E7 (Supplementary Figure S1B). Interestingly, when assessing the proportion of epithelial cells that contributed to each cluster by each patient sample, we found that E7 subcluster did not present in patient 03 (Fig. 2F). Moreover, the pathological diagnosis of this patient was high grade intraepithelial neoplasia, while all the others were carcinoma. This finding further confirmed the hypothesis of E7 as a subcluster of malignant cells. What’s more, the proportion of E7 also demonstrated positive correlations with an increase in the depth of tumor invasion (Supplementary Table S1) rather than the size of microscopic lesion, indicating the potential relationship between E7 and tumor invasion. Notably, the suggested malignant E7 subcluster and potentially malignant E4, E5 and E6 subclusters also presented in paratumor tissues(Supplementary Figure S2B). These results suggested the paratumor tissues may have shared an extensive overlap with intramucosal ESCC.
CNV analysis of epithelial cells was performed and Bayes normalization was applied (Fig. 2G, Supplementary results). Epithelial cells originated from tumor tissues exhibited a relative variation in CNV profiles compared to epithelial cells originated from paratumor tissues(Supplementary Figure S2A). E1 subcluster showed lower variation in CNV compared to other subclusters, suggesting the benign nature of this subcluster. The CNV signal at the anterior end of chromosome 6 was highest in E7 subcluster and also increased in E5 and E6 subclusters. E7 subcluster exhibited a unique gain in CNV at the posterior end of chromosome 20 and the middle of chromosome 3.
We then performed pseudotemporal and PCA analysis and found two evolution fates of esophageal epithelial cells during ESCC tumorigenesis initiated from proliferative Epi4, and then through an intermediate state characterized by E5, to two different ternimal states of normal differentiated E1/E2/E4 or malignant E6/E7, respectively (Fig. 2H). This is also similar to the previously reported cell fates of epithelial cell status transitions in mouse ESCC model9. Besides, expression of several genes of E7 sphingolipid metabolism changed along the trajectory (Fig. 2I).
Metabolic characteristics of epithelial cells in human intramucosal ESCC
Aberrant metabolism is a major hallmark of cancer. Abnormal cancer metabolism, such as aerobic glycolysis and increased anabolic pathways, plays an important role in tumorigenesis, metastasis, drug resistance, and cancer stem cells20.
To understand the metabolic landscape of epithelial cells in intramucosal ESCC, we first computed the score of all 76 active metabolic pathways in all E1-E7 clusters and selected 24 representing metabolic pathways by using the scMetabolism15. We found significant metabolic change trend from E1 to E7. E7(a group of characteristic cells prone to malignancy) showed significantly different characteristics compared with other epithelial cell clusters.
The level of glycolysis, pyruvate metabolism, oxidative phosphorylation, amino acid metabolism (e.g., Tyrosine metabolism) and lipid metabolism (e.g., Steroid biosynthesis) were lower in E7 than other epithelial clusters (Fig. 2J), which is different from progressive cancer 21. Similar to the progressive ESCC, a distinct increase in glycan biosynthesis and metabolism(e.g., Glycosaminoglycan degradation /biosynthesis, other types of O − glycan biosynthesis) was observed in E7.
In order to explore the relationship between E7 and tumor progression, we analyzed the correlationship between the representing metabolic pathaways and tumor phenotypes by computing the classical phenotypic scores of cancer cells(Fig. 2K). Through the correlation of phenotypic score and metabolic pathways, we found that markers of E7 were strongly associated with metastasis and angiogenesis, which is not significant in other phenotypes such as cell proliferation.
Sphingolipid metabolism and glycan biosynthesis/metabolism(e.g., glycosaminoglycan degradation /biosynthesis) showed strong positive correlation with both metastasis and angiogenesis, which is consistent with previous studies22. We further found marker genes such as KDSR, DESG1,PLPP1 and NEU1 after searching for the marker genes related to sphingolipid metabolism in E7. The expression of KDSR was increased along the malignant trajectory, and NEU1 were significantly highlited in the E7 cells as mentioned above (Fig. 2I).. The results of PCR demonstrated that the knockdown of KDSR potentially inhibited EMT of ESCC in vitro by observing a decrease in the expression of snail and slug at the RNA level, indicating a potential target for the development of anti-tumor strategies (Supplementary Figure S2C).
In conclusion, all the results above demonstrated that early ESCC shared both similarities and differences with progressive cancer and ehanced sphingolipid metabolism in E7 showed a strong correlationship with tumor matastasis and angiogenesis.
Exhausted T cells and suppressive immune microenvironment in human intramucosal ESCC
As the vast majority of T cells derived from PBMC, clustering of T cells were independently performed for PBMC, tumor tissues and paratumor tissues, respectively, so as to avoid losing signals from tissue samples.
For PBMC, ten clusters of T cells were characterized with distinctive signatures, namely CD8+ CTL, CD4+ Naïve T cells, CD8+ translational T cells, CD4+ TEMs, γδ T cells, CD4+ Cytotoxic T cells, NKT, CD4+ TCMs, NK-like CD16+ T cells, and Tregs (Fig. 3A). For tissue samples, Four T cell clusters presented both in tumor tissues and paratumor tissues (Fig. 3B, C, Supplementary results), including Tregs, CD8+ CTLs, PD1+CD4+ T cells, CD4+ TMs (Memory T cells). Among them, PD1+CD4+ T cells has previously been reported not to be present in paratumor tissues of advanced ESCC samples23. T cells from tumor tissues and paratumor tissues also exclusively exhibited distinct clusters. CD8+ TMs, MAIT (mucosal-associated invariant T) cells and CD4+ TNs (Naïve T cells) presented in paratumor tissues while CD8+ TEXs (Exhausted T cells), CD4+ISG+ T cells and CD4+ TNs only presented in tumor tissues(Fig. 3D). Additionally, it is worth noticing that CD8 + TEXs expressed a high level of GZMK rather than its typical marker CXCL13, suggesting that CD8 + TEXs in intramucosal ESCC tissues are mainly in an intermediate state, and might go through a ‘Tn to IL7R + Tm to GZMK + TEM/EX’ pathway to become terminal TEX 24. Notebly, The proportions of T cell clusters in paratumor tissues in one of the patient 03 was significantly different from other patients, with a higher proportion of CD4+ TMs, PD1+CD4+ T cells, and MAIT cells, whereas a decreased proportion of CD8+ TMs (Supplementary Figure S3B).
Taken together, the abundance of CD8 + TEXs and enrichment of Tregs and PD1 + CD4 + T cells suggested exhausted and suppressive immune microenvironment in human intramucosal ESCC.
Metabolic charateristics of extausted T cells in human intramucosal ESCC
Recent studies on immunometabolism have demonstrated a close connection between metabolic programing and the specific immune functions they support during both health and disease 25.
To understand the metabolic landscape of T cells in intramucosal ESCC microenvironment, we computed the score of 24 representing metabolic pathways using scMetabolism in all T clusters. We found effector T cells and exthausted T cells (cells normally marked with CD28 family inhibitory receptors), acted abnormally after activation(Fig. 3F).
Exhausted CD8 + T cells were only found in tumor. As mentioned above, exhausted CD8 T cell had a significantly lower metabolic level, including glycolysis, oxidative phosphorylation, TCA cycle, etc (Fig. 3F), which indicated the failure of activation. In contrast, the steroid biosynthesis of these cells was specifically up-regulated, which may help tumor cells escape26. Different from exthausted CD8 + T cells, we observed the enhancement of glycolysis/gluconeogenesis and oxidative phosphorylation in intramucosal ESCC infiltrated PD1 + CD4 + T cells, with the increase of other types of O − glycan biosynthesis, etc (Fig. 3F), indicating the tumor infiltrated PD1 + CD4 + T cells still maintained high metabolic activity, despite the expression of PD1.
To explore the association of metabolic characteristics with phenotypic changes of PD1 + CD4 + T cells, we also carried out the correlation analysis here by computing the classical phenotypic scores of T cells such as interleukin signaling pathway, antigen processing and presentation, inflammation, and immune suppression(Fig. 3G). Interestingly, interleukin signaling pathway is sigificently negatively accociated with several metabolic pathways such as glycolysis, citrate cycle (TCA cycle) and pentose phosphate pathway. It is supposed that PD1 + CD4 + T cells are more active in catabolism rather than anabolism, which makes the cytokines secreted by CD4 + T cells insufficient, resulting in abnormal immune response and inhibition of cell proliferation and growth pathway27.
More importantly, we found that PD1 + CD4 + T cells were associated with significently more numbers in tumors compared to the adjacent tissues, and the top 20 markers of these cells in tumor tissues (logfc > 0.95, P < 0.001) were associated with a lower probability of overall survival in ESCC (Fig. 3H). However, there was no significant correlation between the top markers of PD1 + CD4 + T cells in adjacent tissues and the survival of esophageal cancer. Among them, we observed the specific high expression of CXCL13 and CXCR5 in PD1 + CD4 + T cells in tumor tissues (Fig. 3I).
Characterization of myeloid cell differentiation
Seven monocyte/macrophage cell clusters were categorized (Fig. 4A). Mono1(with high expression of VCAN) and Mono2(with high expression of PADI4) both showed similarity with the classical CD14 + CD16- monocytes. Classically activated macrophages (M1, Macro1, with high expression of CSF1R), alternatively activated macrophages (M2, Macro3, with high expression of ATF3), and myeloid-derived suppressor cells (MDSCs) (Fig. 4B), macrophage clusters were identified using published gene signatures of myeloid cells8, 9. Macro2 may represent the intermediate state, both enriched monocyte, M1 and M2 signature. Macro4 was denoted as tissue resident macrophages (TRM). Interestingly, PADI4, the marker gene of Mono2 subcluster, was strongly associated with a high probability of overall survival in ESCC (Fig. 4C), suggesting it as a prognostic biomarker in ESCC. Human PADI4 was first described in HL-60 cells and was found to participate in differentiation of HL-60 cells into granulocytes and monocytes/macrophages28.
Eight DC clusters were characterized according to discriminative markers, including cDC1, CD1C + LYZ- DC, CD1C + LYZ + DC, pDC, cycling DC, LAMP3 + DC, S100A9 + DC, and IFN induced DC (Fig. 4D, Fig. 4E, Supplementary results). Notebly, cycling DC(subcluster in cell cycle) was enriched in tissues compared with PBMC (Fig. 4F). LAMP3 + DC was recently recognized by several studies 29–31 as a tumor-specific DC subcluster, which was the most activated DC subset with potential migration capacity in tumors. Although some of the subclusters showed different proportions in PBMC and tissues, only LAMP3 + DC was significantly different in quantity between tumor and paratumor tissues (Fig. 4F). Accordingly, LAMP3 + DC was greatly enriched in human intramucosal ESCC specimen, compared to paratumor tissues and PBMC. Among the 8 DC subgroups, LAMP3 + DCs expressed various immune-relevant genes encoding cytokines (IL7R, IL15, IL32), chemokines and chemokine receptors(CCR7, CCL19, and CCL22) (Fig. 4E,G), exhibiting potential to regulate multiple subtypes of lymphocytes32.
Distinct fibroblast subpopulations in human intramucosal ESCC ecosystem
1,279 cancer associated fibroblasts (CAFs) were clustered into 5 subpopulations (Fig. 5A). All 5 subclusters expressed high levels of fibroblast markers such as ACTA2 (a-SMA), PDGFRa, DCN, FN1 and LUM (Fig. 5B), while each subcluster displayed distinct transcriptomic signatures(Fig. 5C, Supplementary Figure S4, Supplementary results). Interestingly, although the GO terms enriched for iCAF1 and iCAF2 have common pathways related to extracellular matrix organization and extracellular structure organization (Fig. 5D), their pathways related to inflammatory response were opposite. iCAF1 participated in the negative regulation of immune system process, while iCAF2 paricipated in complement activation and complement-dependent cytotoxicity. Therefore, the role of iCAF1 was further identified as inhibition in immune modulation and iCAF2 as activation. Moreover, iCAF1 was mainly enriched in ESCC tissues, whereas proportion of iCAF2 was relatively higher in adjacent tissues, which was also consistent with the suppressed immune status in tumorgenesis (Fig. 5E). The active metabolism of iCAF1 was featured with high level of lipid metabolism through scMetabolism(Supplementary Figure S5), including glycosphingolipid biosynthesis and glycerophospholipid metabolism. It was proposed that the ehancement of lipid metabolism of CAFs promotes tumor metastasis33. mCAF1 and mCAF2 expressed high levels of extracellular matrix signatures and GO terms(Fig. 5D). eCAF mainly expressed epithelium-specific marker genes and GO terms, which was in consistent with the EMT-like CAFs in a previous report 34. ScMetabolism also displayed a significant increase in the strong metabolic activity of eCAF(Supplementary Figure S5). Therefore, we speculate that this cluster of cells may be related to the proliferation and biosynthesis of tumor related fibroblasts, and may communicate with tumors. Notably, when combining all CAF subgroups together, GO analysis showed that CAFs derived from tumor specimen manifested different pathway enrichment compared with CAFs derived from paratumor specimen(Fig. 5F). CAFs derived from tumor specimen was associated with higher signals in extracellular matrix organization, extracellular structure organization,antigen processing, and response to oxygen levels. Previous reports have also highlighted a subcluster called vCAF in malignant tumors 34, 35, characterized by microvasculature signature genes such as CD146 (MCAM), MYH11 and RGS5, as well as inflammatory chemokines such as IL-6 and CCL8, with functions mainly in response to hypoxia and mesenchymal cell proliferation. However, this subcluster was not found in human intramucosal ESCC, which may be related to the biobehavior of this very early stage of ESCC.
E7 cells interact with CAFs via the MDK-NCL pathway
To explore the interactions among T cells, epithelial cells and CAFs, we conducted intercellular interaction analyses based on ligand-receptor pairs (Fig. 6A ~ C). Different from progressive ESCC, the interactions between different cell types decreased a lot because of the early stage of the tumor21. In the network of the intramucosal ESCC, interactions between the epithelial cells and CAFs were predicted to be most significant. An increasing intercellular interaction in the MDK (Midkine)-SDC2 (Syndecan-2), MDK-NCL (Nucleolin) and MDK-LRP1(LDL Receptor Related Protein 1) axis was observed between epithelial cells and CAFs, and the interactions between E7 cluster and eCAFs/iCAFs were observed to be the most prominent. In the MDK-NCL axis, E7 showed a particularly strong intercellular interaction with eCAF and mCAF2 (Fig. 6D).
Gypia database was then explored for the expression level of MDK, NCL and LRP1. The results suggested the expression level of MDK and NCL in human esophageal cancer was higher than that in paratumor tissue (Fig. 6E). In correlation analysis using data from TCGA and GTEx, the expression level of MDK was significantly correlated with NCL in human esophageal cancer and paratumor esophageal tissues(Fig. 6F). However, there was no significant difference for LRP1. We then further performed IHC in ESCC specimen of paratumor tissue, high grade dysplasia, intramucosal ESCC and progressive ESCC (Fig. 6G, Supplementary Table S2). Expectedly, with the progressing of ESCC, the expression of MDK and NCL in tumor tissue increased step by step, and the close location indicated the interaction between the CAFs and tumor.
Collectively, these results indicated that MDK was upregulated in esophageal epithelial cells during malignancy change and interacted with NCL in mCAF and eCAF to promote the proliferation of CAFs. The enriched MDK-NCL signal could serve as a sign of CAF activation to stimulate downstream pathways for facilitating tumor invasion, and therefore might be a potential early biomarker of ESCC progression .