Single-cell transcriptomic analysis in a mouse model deciphers cell transition states in the multistep development of esophageal cancer
单细胞转录组分析在小鼠模型中揭示了食管癌多步发展中的细胞过渡状态

Bulk RNA-seq and scRNA-seq of mouse esophageal samples
小鼠食管样本的 bulk RNA-seq 和 scRNA-seq

To explore the transcriptomic alterations at various pathological stages during ESCC tumorigenesis, 8-week-old female C57BL/6 mice were treated with 4NQO for 16 weeks, which resulted in five recognizable precancerous and cancerous lesions in the esophagus, i.e., INF, HYP, DYS, CIS, and ICA (Fig. 1a; Supplementary Fig. 1a). We examined mice receiving 4NQO and found they all developed the expected lesions in the esophagus at different time points of experiment. The ESCC number per animal (mean ± SD) at week 26 was 6.0 ± 3.6. We first performed conventional RNA-seq on mini-bulks of normal epithelial samples obtained by laser-capture microdissection (LCM) from control mice at different ages of 1, 2, 8, and 25 months and various precancerous or cancerous lesions from 4NQO-treated mice. Principal component analysis (PCA) of the differentially expressed genes showed that the expression programs in mice exposed to 4NQO were substantially different from that in control mice; however, there existed some overlaps in the expression profiles between 4NQO-exposed and non-exposed mice. For instance, the transcriptomic profile for stage INF in 4NQO-exposed mice was similar to that in control mice aged 25 months (Fig. 1b), indicating that conventional RNA-seq could not precisely clarify the path of malignant cell transition during the development and progression of ESCC due to high intra-tissue heterogeneity.
为了探索食管鳞状细胞癌（ESCC）发生过程中不同病理阶段的转录组学改变，使用 4-硝基偶氮苯（4NQO）处理 8 周龄的雌性 C57BL/6 小鼠 16 周，结果在食管中观察到五种可识别的癌前和癌变病变，即 INF、HYP、DYS、CIS 和 ICA（图 1a ；补充图 1a ）。我们检查了接受 4NQO 处理的小鼠，发现它们在实验的不同时间点均在食管中发展出了预期的病变。第 26 周每只动物的 ESCC 数量（均值±标准差）为 6.0±3.6。我们首先对通过激光捕获显微切割（LCM）从对照小鼠不同年龄（1、2、8 和 25 个月）获取的正常上皮样本以及 4NQO 处理小鼠的各种癌前或癌变病变中获得的小批量样本进行了常规 RNA-seq 分析。差异表达基因的主成分分析（PCA）显示，暴露于 4NQO 的小鼠的表达程序与对照小鼠存在显著差异；然而，暴露于 4NQO 和未暴露的小鼠之间存在一些重叠的表达谱。 例如，在 4NQO 处理的小鼠中，INF 期的转录组学谱型与 25 个月大的对照小鼠的谱型相似（图 1b ），这表明由于组织内高度的异质性，常规 RNA-seq 无法精确阐明恶性细胞在 ESCC 的发生和发展过程中过渡的路径。

We therefore conducted a time-ordered single-cell transcriptomic profiling on various esophageal lesions in 4NQO-exposed mice. For stage NOR or stage INF, we used the whole esophagus from 30 or 20 mice. For other stages, we used the lesion foci and the sample numbers were 32 HYP from 25 mice, 24 DYS from 17 mice, 24 CIS from 20 mice and 25 ICA from 23 mice, respectively (Fig. 1c). A total of 66,089 cells including 29,975 CD45⁺ immune cells and 36,114 CD45⁻ non-immune cells were obtained across various pathological stages (Fig. 1c; Supplementary Fig. 1b). The median unique molecular identifier (UMI) per cell was 7748 in immune cells and 9370 in non-immune cells (Supplementary Fig. 1c); with a median signal detect ability of 1936 genes for immune cells and 2620 genes for non-immune cells (Supplementary Fig. 1c), respectively. Based on the expressions of canonical markers, we classified immune cells into T cells, B cells, myeloid cells and natural killer cells (Fig. 1c; Supplementary Fig. 1d, e) and identified four clusters of non-immune cells including epithelial cells, fibroblasts, endothelial cells and myocytes by using t-distributed stochastic neighbor embedding (tSNE) (Fig. 1c; Supplementary Fig. 1f, g).
因此，我们在 4NQO 处理的小鼠不同食管病灶中进行了时间顺序的单细胞转录组测序。对于 NOR 阶段或 INF 阶段，我们使用了 30 只或 20 只小鼠的整个食管。对于其他阶段，我们使用了病灶区域，样本数量分别为：32 个 HYP（25 只小鼠），24 个 DYS（17 只小鼠），24 个 CIS（20 只小鼠）和 25 个 ICA（23 只小鼠）（图 1c ）。总共获得了 66,089 个细胞，包括 29,975 个 CD45 ⁺ 免疫细胞和 36,114 个 CD45 ⁻ 非免疫细胞（图 1c ；补充图 1b ）。免疫细胞的中位唯一分子标识符（UMI）为 7748，非免疫细胞为 9370（补充图 1c ）；免疫细胞的中位信号检测能力为 1936 个基因，非免疫细胞为 2620 个基因（补充图 1c ）。根据经典标记物的表达，我们将免疫细胞分类为 T 细胞、B 细胞、髓系细胞和自然杀伤细胞（图 1c ；补充图 1d, e 并通过 t 分布随机邻居嵌入（tSNE）识别了四种非免疫细胞簇，包括上皮细胞、成纤维细胞、内皮细胞和肌细胞（图 1c ；补充图 1f, g ）。

Identifying epithelial cell types during ESCC tumorigenesis
在 ESCC 肿瘤发生过程中识别上皮细胞类型

To discover how normal esophageal epithelium develops into invasive carcinoma, we next examined the expression alterations and functional changes in epithelial cells during the transition from normal to precancer or cancer. We identified 1756 epithelial cells across all six stages that were classified into six subtypes designated as EpiC 1 to EpiC 6 (Fig. 2a; Supplementary Fig. 2a and Supplementary Table 1). Through the analysis of pathway activities (Fig. 2b; Supplementary Fig. 2b), we found that EpiC 1 (n = 339) had higher expression of genes (e.g., Birc5, Mki67, Top2a, and Ube2c) in mitosis and proliferation^8–11 compared with other cell clusters, probably reflecting the ability of routine self-renewal of the esophageal epithelium. EpiC 2 (n = 463) had significantly higher expression of genes such as Gstm1, Gstm2, Adh7, and Aldh3a1, which were involved in the xenobiotic detoxification, likely due to the 4NQO exposure^12–14. EpiC 3 (n = 421) had substantially higher expression of genes in response to extracellular carcinogen stimuli, pro-inflammation and Hif-1 signaling (e.g., Jun, Junb, Zfp36, and Atf3)^15–18. EpiC 4 (n = 110) had a higher expression program for keratinization (e.g., Krt4, Krt13, S100a8, and S100a9)^19,20. EpiC 5 (n = 213) featured by the upregulation of genes in the epithelial–mesenchymal transition (EMT) pathways (e.g., Mmp13 and Itga6)^21,22 and angiogenesis (e.g., Vegfa)²³. EpiC 6 (n = 210) had significantly upregulated expression in the genes controlling cancer invasiveness (e.g., Mmp14 and Ecm1)^24,25 and metastasis (e.g., Tm4sf1 and Vim)^26,27. Besides, EpiC 5 cells also showed an expression pattern resembling the combination of EpiC 3, EpiC 4, and EpiC 6.
为了发现正常食管上皮如何发展为浸润性癌，我们接下来检查了从正常到癌前或癌变过程中上皮细胞的表达变化和功能变化。我们识别出了六个阶段中的 1756 个上皮细胞，并将其分为六种亚型，分别命名为 EpiC 1 到 EpiC 6（图 2a; 补充图 2a 和补充表 1 ）。通过通路活性分析（图 2b ；补充图 2b ），我们发现 EpiC 1（n=339）在有丝分裂和增殖（例如，Birc5、Mki67、Top2a 和 Ube2c）基因的表达上高于其他细胞簇，这可能反映了食管上皮的常规自我更新能力。EpiC 2（n=463）在异物解毒基因（例如，Gstm1、Gstm2、Adh7 和 Aldh3a1）的表达上显著高于其他细胞簇，这可能归因于 4NQO 暴露。EpiC 3（n=421）在对细胞外致癌刺激、促炎和 Hif-1 信号（例如，Jun、Junb、Zfp36 和 Atf3）的基因表达上显著增加。 EpiC 4（n=110）在角化方面具有更高的表达程序（例如 Krt4、Krt13、S100a8 和 S100a9） ^19,20 。EpiC 5（n=213）在上皮-间质转化（EMT）途径相关基因（例如 Mmp13 和 Itga6） ^21,22 和血管生成相关基因（例如 Vegfa） ²³ 的上调表达中表现出特征。EpiC 6（n=210）在控制癌症侵袭性（例如 Mmp14 和 Ecm1） ^24,25 和转移（例如 Tm4sf1 和 Vim） ^26,27 的基因表达上调方面显著。此外，EpiC 5 细胞还表现出类似于 EpiC 3、EpiC 4 和 EpiC 6 的表达模式。

Three subtypes, EpiC 1–3, were all present in normal esophageal epithelium and at the various stages of 4NQO-induced abnormal epithelia, although the proportions of cells at each stage were different (Fig. 2c), suggesting that these subtypes of epithelial cells were the background components of esophageal tissues. EpiC 4 and EpiC 5 appeared from INF and HYP, respectively, while EpiC 6 was present uniquely at ICA as the major cell subtype (Fig. 2c). These cell subtype dynamics well concorded with the transitional role of EpiC 5 and the malignancy of EpiC 6 and prompted us to examine the cell status transition more closely. A single-cell diffusion map based on the gene expression similarity revealed a clear cell status transition across all pathological stages (Fig. 2d; Supplementary Fig. 2c). Gene expression vectors pointed from stage NOR cells to the centroid of cells from each stage. It was clear that major directional changes happened before the HYP and ICA stages, whereas the transition from stage HYP to stage CIS was mainly the expansion along the same direction. Together, these results indicated an important role of INF to HYP transition in ESCC development.
三个亚型 EpiC 1-3 在正常食管上皮和 4NQO 诱导的异常上皮的各个阶段中均存在，尽管每个阶段的细胞比例不同（图 2c ），这表明这些上皮细胞亚型是食管组织的背景成分。EpiC 4 和 EpiC 5 分别从 INF 和 HYP 出现，而 EpiC 6 仅在 ICA 阶段作为主要细胞亚型存在（图 2c ）。这些细胞亚型的变化与 EpiC 5 的过渡作用以及 EpiC 6 的恶性程度高度一致，促使我们更仔细地研究细胞状态的转变。基于基因表达相似性的单细胞扩散图显示，在所有病理阶段中均存在明显的细胞状态转变（图 2d ；补充图 2c ）。基因表达向量从 NOR 阶段的细胞指向每个阶段细胞的质心。很明显，在 HYP 和 ICA 阶段之前发生了主要的方向变化，而从 HYP 阶段到 CIS 阶段的转变主要是沿着相同方向的扩展。综上所述，这些结果表明 INF 到 HYP 的转变在食管癌的发展中起着重要作用。

Further analysis revealed that the six representative genes selected from each epithelial cluster were differentially distributed among disease stages (Fig. 2e). Aldh3a1 and Atf3 were expressed across all stages and the levels were significantly higher at stage INF than that at the advancing stages. High level of S100a8 appeared at stage HYP and covered all precancerous and ICA stages whereas the highest levels of Itga6 and Mmp14 were seen at stage ICA although their expressions were also detected in cells across all stages. The dynamic expressions of these genes at protein level in mice esophageal tissues with different disease stages were compared by using immunohistochemistry and the results were generally in line with their RNA expression despite of some disparity (Fig. 2f; Supplementary Fig. 2d). The abrupt upregulation of S100a8 in cells at stage HYP suggests a dramatic transition related to immune response.
进一步分析发现，从每个上皮簇中选择的六个人代表基因在不同疾病阶段中的分布情况不同（图 2e ）。Aldh3a1 和 Atf3 在所有阶段均表达，且在 INF 阶段的表达水平显著高于进展阶段。S100a8 在 HYP 阶段的表达水平较高，并覆盖了所有癌前病变和 ICA 阶段，而 Itga6 和 Mmp14 在 ICA 阶段的表达水平最高，尽管它们在所有阶段的细胞中也有检测到。通过免疫组化比较不同疾病阶段小鼠食管组织中这些基因在蛋白水平上的动态表达，结果与 RNA 表达情况总体一致，尽管存在一些差异（图 2f; 补充图 2d ）。S100a8 在 HYP 阶段细胞中的突然上调表明与免疫反应相关的剧烈转变。

Identifying cell fates of epithelial cell status transitions
识别上皮细胞状态转换的命运

We performed pseudotime and PCA analysis and found two evolution fates of esophageal epithelial cells during ESCC tumorigenesis both starting from EpiC 1 cells that had the lowest pseudotime value. Some cells transformed from proliferative EpiC 1 to normal differentiated EpiC 4 while other cells transformed to malignant EpiC 6, processing through EpiC 2 to EpiC 5 cells (Fig. 3a; Supplementary Fig. 3a, b). The evolution of EpiC 1 to EpiC 6 was mainly along component 1. Gene set variation analysis (GSVA) of component 1 revealed a significant enrichment of genes related to cell invasiveness, EMT and angiogenesis (Fig. 3b, c; Supplementary Fig. 3c), which was concordant with the expression programs of EpiC 6 cells (Fig. 2b). As EpiC 6 cells appeared only at the ICA stage, these results implied that component 1 might be the underlying molecular mechanism for malignant transition of the esophageal tissues (Supplementary Fig. 3d).
我们在假时间分析和 PCA 分析中发现，在 ESC 癌变过程中，从具有最低假时间值的 EpiC 1 细胞开始，存在两种进化命运。一些细胞从增殖的 EpiC 1 转变为正常分化 EpiC 4，而其他细胞则转变为恶性 EpiC 6，经过 EpiC 2 到 EpiC 5 细胞（图 3a ；补充图 3a, b ）。从 EpiC 1 到 EpiC 6 的进化主要沿着成分 1。成分 1 的基因集变异分析（GSVA）揭示了与细胞侵袭性、EMT 和血管生成相关的基因显著富集（图 3b, c; 补充图 3c ），这与 EpiC 6 细胞的表达程序一致（图 2b ）。由于 EpiC 6 细胞仅在 ICA 阶段出现，这些结果表明成分 1 可能是食管组织恶性转化的潜在分子机制（补充图 3d ）。

We then examined whether the alterations of any transcription factors (TFs), well-documented ESCC-related mutation, or methylation dysregulation were included in the oncogenic evolution along component 1. We found that the expression levels of Pitx1, Trp53, and Bclaf1, which are known tumor suppressor TFs^28–30, were substantially decreased while TFs promoting EMT (e.g., Ets1 and Snai3)^31,32, immunosuppression (e.g., Eomes)³³, cell migration and invasion (e.g., Creb3 and Elk3)^34,35 were significantly upregulated. We also found two previously reported methylation-regulated genes in human ESCC, Rab25³⁶, and Met³⁷, were substantially down- or upregulated, respectively, along component 1. Furthermore, the expression levels of Notch1, a tumor suppressor frequently mutated in human ESCC³⁸, displayed a gradual decrease over time along the evolution of component 1 (Fig. 3d; Supplementary Fig. 3e).
我们随后检查了任何转录因子（TFs）的改变，这些 TFs 是已知与食管鳞状细胞癌（ESCC）相关的突变，或甲基化失调是否包括在沿成分 1 的致癌进化过程中。我们发现已知的肿瘤抑制 TFs Pitx1、Trp53 和 Bclaf1 的表达水平显著降低，而促进上皮-间质转化（EMT）的 TFs（例如，Ets1 和 Snai3） ^31,32 、免疫抑制（例如，Eomes） ³³ 、细胞迁移和侵袭的 TFs（例如，Creb3 和 Elk3） ^34,35 的表达水平显著上调。我们还发现两个之前在人类 ESCC 中报道的甲基化调节基因 Rab25 ³⁶ 和 Met ³⁷ 在沿成分 1 的进化过程中分别显著下调或上调。此外，肿瘤抑制基因 Notch1 在人类 ESCC 中经常发生突变 ³⁸ ，其表达水平在成分 1 的进化过程中逐渐下降（图 3d ；补充图 3e ）。

Genes in component 2 were mainly related to cell proliferation (e.g., Birc5, Mki67, and Top2a). Cell cycle regression analysis and GSVA scoring revealed a correlation between proliferative ability and value of component 2 (Fig. 3e, f). We found that EpiC 6 cells had the highest expression of proliferation-related genes and the highest proportion of cells at G2/M status as compared with other clusters except EpiC 1 cells (Fig. 3e, f; Supplementary Fig. 3f–h), indicating increased proliferative ability of this cell cluster. On the other hand, EpiC 4 cells were at the G1 and S status (Supplementary Fig. 3g, h), showing reduced proliferative capacity. Together, these unsupervised analyses depicted two clear cell fates during ESCC carcinogenesis and indicated that attacked by the carcinogen, a portion of esophageal epithelial cells went into oncogenic route (from EpiC 1 and EpiC 2 to EpiC 6) but not into normal differentiation route (from EpiC 1 and EpiC 2 to EpiC 4). To exclude the impact of potential uncertainty of pathologically staging of early lesions such as hyperplasia and dysplasia, we also performed a sensitive analysis by excluding cells in week 20 or week 22. As a result, the cell state transition was similar to that without the exclusion (Supplementary Fig. 3i, j).
组件 2 中的基因主要与细胞增殖相关（例如，Birc5、Mki67 和 Top2a）。细胞周期回归分析和 GSVA 评分揭示了增殖能力和组件 2 值之间的关联（图 3e, f ）。我们发现，EpiC 6 细胞在与增殖相关的基因表达和 G2/M 期细胞比例方面最高，与其他簇相比，除了 EpiC 1 细胞（图 3e, f ；补充图 3f–h ），表明该细胞簇的增殖能力增强。另一方面，EpiC 4 细胞处于 G1 和 S 期（补充图 3g, h ），显示出降低的增殖能力。综上所述，这些无监督分析描绘了 ESCC 致癌过程中两种明确的细胞命运，并表明致癌物质攻击时，一部分食管上皮细胞进入了致癌途径（从 EpiC 1 和 EpiC 2 到 EpiC 6），而不是正常分化途径（从 EpiC 1 和 EpiC 2 到 EpiC 4）。为了排除早期病变如增生和异型增生的病理分期潜在不确定性的影响，我们还进行了敏感性分析，排除了第 20 周或第 22 周的细胞。 因此，细胞状态的转换与不排除情况下的结果相似（补充图 3i, j ）。

We further investigated the key pathways driving epithelial cells from stage INF to stage HYP and found that non-malignant EpiC 1–4 epithelial cells had significantly elevated expressions of the genes involving in carcinogen detoxification such as Nqo1, Aldh3a1, and Gstp1 (Fig. 3g), which reflected normal cellular response to the damage induced by 4NQO. The continuous damage might induce immune response via the stimulator of interferon genes (STING), AIM2 and NF-κB signaling because the expression levels of Tmem173, Aim2, Nfkb1, and Nfkb2 were significantly elevated (Fig. 3g). In addition, we observed substantial differences in the expressions of NF-κB downstream genes in epithelial cell clusters (Fig. 3h). EpiC 3 cells had an increased expression of some transcription-related genes such as Junb and Fos, whereas EpiC 5 cells showed higher expression levels of inflammation genes (e.g., Cxcl1, Il6, Tnf, and Csf3) compared with other EpiCs, implying that they had crosstalks with immune cells. EpiC 5 cells also had high level of Hif1a, which may lead to metabolism remodeling that promoted cell transformation. As malignant cells, EpiC 6 showed strong activation of survival-related genes (e.g., Bcl2 and Xiap) (Fig. 3g, h). These results indicated that under 4NQO attack, some epithelial cells might endure the damage and activate the inflammatory response, resulting in consequent cell survival and transformation.
我们进一步研究了从 INF 阶段到 HYP 阶段上皮细胞的关键通路，发现非恶性 EpiC 1-4 上皮细胞中涉及致癌物解毒的基因（如 Nqo1、Aldh3a1 和 Gstp1）的表达显著升高（图 3g ），这反映了这些细胞对 4NQO 诱导的损伤的正常细胞反应。持续的损伤可能通过干扰素基因刺激器（STING）、AIM2 和 NF-κB 信号通路引发免疫反应，因为 Tmem173、Aim2、Nfkb1 和 Nfkb2 的表达水平显著升高（图 3g ）。此外，我们观察到上皮细胞簇中 NF-κB 下游基因的表达存在显著差异（图 3h ）。EpiC 3 细胞中一些与转录相关的基因（如 Junb 和 Fos）的表达增加，而 EpiC 5 细胞则表现出更高的炎症基因（如 Cxcl1、Il6、Tnf 和 Csf3）的表达水平，表明它们与免疫细胞存在相互作用。EpiC 5 细胞中 Hif1a 的水平也较高，这可能导致代谢重编程，促进细胞转化。 作为恶性细胞，EpiC 6 显示出强烈的生存相关基因（例如，Bcl2 和 Xiap）激活（图 3g, h ）。这些结果表明，在 4NQO 攻击下，一些上皮细胞可能耐受损伤并激活炎症反应，从而导致随后的细胞生存和转化。

Transcriptomic changes in fibroblasts promote tumorigenesis
成纤维细胞的转录组变化促进肿瘤发生

We next examined the transcriptomic alterations of cells in the tissue microenvironment and identified eight fibroblast clusters (FibCs) by different gene expression patterns (Fig. 4a; Supplementary Fig. 4a, b) that showed similar expression pathways and dynamics to epithelial cells (Fig. 4b–d). FibC 1 was the dominant composition in stage NOR and then markedly decreased with the processing of tumorigenesis (Fig. 4c, d). FibC 3 cells were the major fibroblasts in stage INF and had a strong interferon response and MHC-mediated antigen presentation activity (Fig. 4b), representing the initial immune response to tissue damage. FibC 4 and FibC 6, the major fibroblasts in stage DYS, had the elevated expression of genes involved in the TGF-β, STAT5 induced by IL-2 or STAT3 induced by IL-6 signaling pathways. FibC 8 cells increased along the tumorigenic stages with high expressions of genes in Myc and angiogenesis pathways. The replacement of FibC 3 by FibC 4 and FibC 6 during INF to HYP transition further confirmed the shift of immune response during early ESCC development. Specifically, beginning from stage HYP, fibroblasts actively recruited immune cells through increasing the expressions of complements and chemokines (Fig. 4e). For example, we found that FibC 6 and FibC 8 had significantly elevated expressions of various chemokine-related genes (e.g., Ccl7 in FibC6 and Cxcl12 in FibC 8, Supplementary Fig. 4c), which were well-known molecules that recruit immune cells. The dynamic changes of fibroblast clusters suggested that the immune response is modulated during ESCC tumorigenesis.
我们接下来检查了组织微环境中细胞的转录组变化，并通过不同的基因表达模式识别出了八种成纤维细胞簇（FibCs）（图 4a; 补充图 4a, b ），这些簇显示出与上皮细胞相似的表达途径和动态变化（图 4b–d ）。在 NOR 阶段，FibC 1 是主要组成成分，然后随着肿瘤发生的发展显著减少（图 4c, d ）。在 INF 阶段，FibC 3 是主要的成纤维细胞，具有强烈的干扰素反应和 MHC 介导的抗原呈递活性（图 4b ），代表了对组织损伤的初始免疫反应。在 DYS 阶段，FibC 4 和 FibC 6 是主要的成纤维细胞，这些细胞中与 TGF-β、IL-2 诱导的 STAT5 信号或 IL-6 诱导的 STAT3 信号通路相关的基因表达升高。FibC 8 细胞在肿瘤发生阶段中表达 Myc 和血管生成途径相关基因的水平升高。INF 到 HYP 过渡期间 FibC 3 被 FibC 4 和 FibC 6 取代进一步证实了早期食管癌发展中免疫反应的转变。 特别是在从 HYP 阶段开始，成纤维细胞通过增加补体和趋化因子的表达积极招募免疫细胞（图 4e ）。例如，我们发现 FibC6 和 FibC8 显著上调了多种趋化因子相关基因的表达（例如，FibC6 中的 Ccl7 和 FibC8 中的 Cxcl12，补充图 4c ），这些是已知能够招募免疫细胞的分子。成纤维细胞簇动态变化表明，在食管鳞状细胞癌（ESCC）发生过程中，免疫反应被调节。

Turndown of adaptive anticancer immune during tumorigenesis
肿瘤发生期间适应性抗肿瘤免疫的下调

We then explored the T cell status during ESCC tumorigenesis and identified 4 clusters designated as CD8⁺ T cell, CD4⁺ T cell, CD4⁻CD8⁻ T cell 1 and CD4⁻CD8⁻ T cell 2, respectively (Fig. 5a). We further classified CD8⁺ T cells into 7 clusters (Fig. 5b) and CD4⁺ T cells into 7 clusters (Fig. 5c), respectively. During 4NQO-induced tumorigenesis, all CD8⁺ T cells expressed naive and memory marker genes (e.g., Tcf7, Il7r, and Ccr7; Supplementary Fig. 5a)^39–41 rather than cytotoxic markers. In stage NOR, the major type was naive cells (CD8-TN1) while in stage INF, the major cell type was replaced by memory T cells, mainly CD8⁺ tissue resident memory T (CD8-T_RM) cells and effector memory T (CD8-T_EM) cells (Fig. 5d), suggesting that the immune responses were active in these two stages. Furthermore, we observed a high expression of the exhaustion markers (e.g., Pdcd1 and Lag3)^42,43 in addition to the effective markers (e.g., Gzmk, Gzmb, and Prf1)^44,45 in these memory T cells (Fig. 5e); but the proportion of these memory T cells were declined along with tumorigenic processing after stage INF, indicating a non-effective CD8⁺ T cell dominant microenvironment throughout precancerous stages.
我们随后探索了食管鳞状细胞癌（ESCC）肿瘤发生过程中 T 细胞的状态，并识别出 4 个簇，分别命名为 CD8 ⁺ T 细胞、CD4 ⁺ T 细胞、CD4 ⁻ CD8 ⁻ T 细胞 1 和 CD4 ⁻ CD8 ⁻ T 细胞 2（图 5a ）。我们进一步将 CD8 ⁺ T 细胞分为 7 个簇（图 5b ），将 CD4 ⁺ T 细胞分为 7 个簇（图 5c ）。在 4NQO 诱导的肿瘤发生过程中，所有 CD8 ⁺ T 细胞均表达幼稚和记忆标记基因（例如，Tcf7、Il7r 和 Ccr7；补充图 5a ）而非细胞毒性标记基因。在 NOR 阶段，主要类型是幼稚细胞（CD8-TN1），而在 INF 阶段，主要细胞类型被记忆 T 细胞取代，主要是 CD8 ⁺ 组织驻留记忆 T 细胞（CD8-T _RM ）和效应记忆 T 细胞（CD8-T _EM ）（图 5d ），表明这两个阶段的免疫反应是活跃的。此外，我们观察到这些记忆 T 细胞中除了有效标记基因（例如，Gzmk、Gzmb 和 Prf1） ^44,45 的高表达外，还存在耗竭标记基因（例如，Pdcd1 和 Lag3） ^42,43 的高表达（图 5e ）；但在肿瘤发生阶段 INF 之后，这些记忆 T 细胞的比例下降，表明整个癌前阶段存在非有效的 CD8 ⁺ T 细胞主导微环境。

Considering CD4⁺ T cells (Fig. 5f; Supplementary Fig. 5b), we noticed an imbalance among CD4-Th1-like, CD4-Th2, and CD4-Th17 cells, involving in three major kinds of cell-mediated effector immunity, respectively⁴⁶. Beginning at stage HYP, the proportion of CD4-Th1-like cells was decreased and replaced by CD4-Th2 and CD4-Th17 cells (Fig. 5f, g). CD4-Th1-like cells expressing high levels of Tbx21, Gzmb, and Ifng (Fig. 5h; Supplementary Fig. 5c) had been shown to have anti-tumor capability through activating monocytes and responding to interferon-γ⁴⁶. Differential gene expression analysis also revealed that some genes relative to immune activation (e.g., Ifitm1, Ifitm2, and Rora)^47,48 were significantly higher in CD4-Th1-like cells than in CD4-Th2 cells (Supplementary Fig. 5d), a type of immune cells that their role in tumorigenesis was uncertain. The activation of CD4-Th17 produced high levels of Il22, Il17a, and Rorc at stage HYP (Fig. 5h) that may activate monocytes and recruit neutrophils to respond to epithelial injuries⁴⁶. The proportion of CD4-Th17 cells increased sharply after stage INF (Fig. 5f). The reduced anti-tumor role of CD4⁺ T cells and increased inflammatory response suggest an immune response Type 1 to Type 3 transition⁴⁶ in the esophageal microenvironment during tumorigenesis.
考虑到 CD4 ⁺ T 细胞（图 1# 补充图 2#），我们注意到 CD4-Th1 样、CD4-Th2 和 CD4-Th17 细胞之间存在不平衡，分别涉及三种主要的细胞介导效应免疫类型。从 HYP 阶段开始，CD4-Th1 样细胞的比例下降，并被 CD4-Th2 和 CD4-Th17 细胞所取代（图 4#）。表达高水平 Tbx21、Gzmb 和 Ifng 的 CD4-Th1 样细胞（图 5#；补充图 6#）已被证明具有抗肿瘤能力，通过激活单核细胞并响应干扰素-γ发挥作用 ⁴⁶ 。差异基因表达分析还显示，一些与免疫激活相关的基因（如 Ifitm1、Ifitm2 和 Rora） ^47,48 在 CD4-Th1 样细胞中的表达显著高于 CD4-Th2 细胞（补充图 9#），而 CD4-Th2 细胞在肿瘤发生中的作用尚不确定。CD4-Th17 的激活在 HYP 阶段产生了高水平的 Il22、Il17a 和 Rorc（图 10#），可能激活单核细胞并招募中性粒细胞以响应上皮损伤 ⁴⁶ 。CD4-Th17 细胞的比例在 INF 阶段后急剧增加（图 12#）。 CD4+ T 细胞的抗肿瘤作用减弱和炎症反应增加表明在肿瘤发生过程中食管微环境中免疫反应从类型 1 向类型 3 过渡。

Activated myeloid cells create inflammatory microenvironment
激活的髓系细胞创建炎症微环境

We found that, among the 11 clusters of myeloid cells (Fig. 6a; Supplementary Fig. 6a), the proportions of monocyte/macrophage (Mo/Mφ)-C 1, Mo/Mφ-C 3, cDC2-C 1, tolerogenic dendritic cell (tDC), and neutrophils changed substantially along the stages of tumorigenesis (Fig. 6b, c). After stage INF, Mo/Mφ-C 1 cells were replaced by Mo/Mφ-C 3 cells, featured by decreased expression levels of the genes (e.g., H2-Ab1 and C1qa) relative to antigen presentation and complement (Supplementary Fig. 6b)^49,50. Meanwhile, the expression level of genes involved in immunosuppression (e.g., Il10, Arg1, and Adora2b, Fig. 6d)^51–53 elevated along with the increase in Mo/Mφ-C 3 proportion. Furthermore, the proportion of cDC2-C 1 sharply decreased since stage NOR with the downregulation of Icosl, a gene representing the immune stimulation activity⁴³ (Fig. 6d). In contrast, the immune suppressive tDC cells with high expressions of Cd274, Pdcd1lg2, and Ido1 were persisted throughout all the tumorigenic stages (Fig. 6d). In addition, we found that the number of neutrophils continuously increased, and these cells expressed high levels of Mmp9, Prok2, Vegfa, and Nos2, but not tumor suppressors (Supplementary Fig. 6c). These results suggested that myeloid cells were activated during the processes of ESCC tumorigenesis, resulting in an inflammatory microenvironment and CD8⁺ T-associated adaptive immune suppression, which might allow initiated ESCC cells to survive and proliferate. Transcriptomic analysis also revealed that the immune suppression microenvironment might be created by the complex chemotaxis among various types of cells in this microenvironment. For example, neutrophils recruitment could be mediated through its high expression of Cxcr1/2 receptors by their ligands Cxcl1/2/5 from other cell types or by their brother neutrophils through the simultaneous expression of Ccl3 and Ccr1. Similarly, tDC, Th17, and CD4⁻CD8⁻ T cells could be recruited by chemokines secreted by different cell types (Fig. 6e).
我们发现，在 11 个髓系细胞簇（图 6a; 补充图 6a ）中，肿瘤发生不同阶段中单核细胞/巨噬细胞（Mo/Mφ）-C1、Mo/Mφ-C3、cDC2-C1、耐受性树突状细胞（tDC）和中性粒细胞的比例发生了显著变化（图 6b, c ）。在 INF 阶段之后，Mo/Mφ-C1 细胞被 Mo/Mφ-C3 细胞取代，后者表现出与抗原呈递和补体相关的基因（如 H2-Ab1 和 C1qa）表达水平降低（补充图 6b ） ^49,50 。同时，与免疫抑制相关的基因（如 Il10、Arg1 和 Adora2b，图 6d ） ^51–53 的表达水平随着 Mo/Mφ-C3 比例的增加而升高。此外，cDC2-C1 的比例自 NOR 阶段以来急剧下降，伴随 Icosl 基因表达下调，该基因代表了免疫刺激活性（图 6d ） ⁴³ 。相比之下，高表达 Cd274、Pdcd1lg2 和 Ido1 的免疫抑制性 tDC 细胞在整个肿瘤发生阶段持续存在（图 6d ）。 此外，我们发现中性粒细胞的数量持续增加，这些细胞表达了高水平的 Mmp9、Prok2、Vegfa 和 Nos2，但不表达肿瘤抑制基因（补充图 6c ）。这些结果表明，在食管鳞状细胞癌（ESCC）肿瘤发生过程中，髓系细胞被激活，导致炎症微环境和 CD8 ⁺ T 细胞相关的适应性免疫抑制，这可能使启动的 ESCC 细胞得以存活和增殖。转录组分析还揭示，这种免疫抑制微环境可能是由微环境中各种细胞复杂的趋化作用共同创造的。例如，中性粒细胞的招募可能是通过其高表达的 Cxcr1/2 受体由其他细胞类型的 Cxcl1/2/5 配体或通过其兄弟中性粒细胞同时表达 Ccl3 和 Ccr1 而介导的。同样，tDC、Th17 和 CD4 ⁻ CD8 ⁻ T 细胞可以通过不同细胞类型分泌的趋化因子被招募（图 6e ）。

Besides the roles of chemokines, cells in the esophageal microenvironment were also stimulated by each other through various signaling pathways. We found that from stage HYP, the expression levels of Il17a in CD4⁺ T cells and Il17ra in myeloid cells were elevated (Fig. 6f), which were well-known to have ability to activate macrophages and neutrophils⁵⁴. On the other hand, we found significant correlations between Il1b and Il23a levels in macrophages and neutrophils and their receptors Il1r1 and Il23r in CD4⁺ T cells, implying close-loop activation among these cells. The activated immune cells secreted high levels of cytokines (Il1b, Il6, Il18, Tnf) and formed an inflammatory microenvironment starting from stage HYP (Fig. 6f). Analyses of interaction indicated increased interactions between malignant epithelial cells (EpiC 5 and EpiC 6) and inflammatory immune cells, Th17 and neutrophils (Supplementary Fig. 6d, e). The inflammatory microenvironment might further promote carcinogen-initiated epithelial cells to transit along the carcinogenic roadmap as illustrated in Figs. 2b and 4b.
除了趋化因子的作用外，食管微环境中的细胞还通过各种信号通路相互刺激。我们发现从 HYP 阶段开始，CD4 T 细胞中的 Il17a 表达和髓系细胞中的 Il17ra 表达升高（图 1），这两种因子已知能够激活巨噬细胞和中性粒细胞。另一方面，我们发现巨噬细胞和中性粒细胞中的 Il1b 和 Il23a 水平与其受体 Il1r1 和 Il23r 在 CD4 T 细胞中的表达之间存在显著相关性，表明这些细胞之间存在紧密的激活循环。激活的免疫细胞分泌高水平的细胞因子（Il1b、Il6、Il18、Tnf），从 HYP 阶段开始形成炎症微环境（图 4）。相互作用分析表明，恶性上皮细胞（EpiC 5 和 EpiC 6）与炎症免疫细胞（Th17 和中性粒细胞）之间的相互作用增加（补充图 5）。炎症微环境可能进一步促进致癌物引发的上皮细胞沿着致癌路径过渡，如图 6b 和图 7 所示。

Validation of the expression profiles in human ESCC samples
人类 ESCQ 样本中表达谱的验证

We last performed a validation examination of the identified expression profiles with mouse model in tissue samples obtained from four ESCC patients having both precancerous and cancerous lesions in the esophagus. RNA-seq was conducted with mini-bulk tissues of invasive carcinoma (ICA, n = 9), carcinoma in situ (CIS, n = 6), dysplasia (DYS, n = 17) and tumor-adjacent normal epithelia (taNOR, n = 8) obtained by LCM (Supplementary Fig. 7a). PCA showed a clearly different expression profile between cancer tissues (CIS and ICA) and non-cancer tissues (taNOR and DYS); however, an extensive overlap was present between taNOR and DYS (Fig. 7a; Supplementary Fig. 7b). In addition, we found that the six representative genes identified in mouse different epithelial cell clusters were also expressed in the precancerous lesions from ESCC patients (Fig. 7b) and the expression patterns were in line with those in mice. ALDH3A1 and S100A8 were highly expressed in the early precancerous lesions while MMP14 and ITGA6 were mainly present in invasive ESCC, which was concordant with their biological functions. As LCM captured mainly epithelial cells and the samples were not suitable for analyzing the status of other cell types, we performed bulk RNA-seq using biopsy samples. The results also showed a trend of increased neutrophils in the microenvironment, although this trend along the stage progress were marginally significant (trend test, P = 0.066; Fig. 7c). These results suggested that strongly supported our findings in mouse model and provided new clues in understanding ESCC initiation.
我们最后使用来自 4 例同时具有癌前病变和癌变病变的食管鳞状细胞癌（ESCC）患者组织样本的 mini-bulk 组织进行了已鉴定表达谱的验证检查。进行了 RNA-seq 分析，样本包括浸润性癌（ICA，n=9）、原位癌（CIS，n=6）、异型增生（DYS，n=17）和肿瘤相邻正常上皮（taNOR，n=8），这些样本通过激光捕获显微切割（LCM）获得（补充图 7a ）。主成分分析（PCA）显示癌组织（CIS 和 ICA）和非癌组织（taNOR 和 DYS）之间有明显不同的表达谱；然而，taNOR 和 DYS 之间存在广泛的重叠（图 7a; 补充图 7b ）。此外，我们发现，在小鼠不同上皮细胞簇中鉴定出的 6 个代表性基因也在 ESCC 患者的癌前病变中表达（图 7b ），其表达模式与小鼠中的模式一致。ALDH3A1 和 S100A8 在早期癌前病变中高度表达，而 MMP14 和 ITGA6 主要存在于侵袭性 ESCC 中，这与其生物学功能一致。 LCM 主要捕获了上皮细胞，且样本不适合分析其他细胞类型的状态，因此我们使用活检样本进行了 bulk RNA-seq。结果也显示微环境中中性粒细胞的趋势增加，尽管这种趋势沿阶段进展仅边缘显著（趋势检验，P=0.066；图 7c ）。这些结果支持了我们在小鼠模型中的发现，并提供了理解食管癌起始的新线索。

Human biospecimen collection
人类生物样本采集

ESCC tumor, dysplasia lesions and tumor-adjacent (>5 cm) normal tissues of the same patients (n = 4) used for LCM and esophageal lesions of different pathological stages (n = 45) used for bulk RNA sequencing were collected during surgery or endoscopy in Linzhou Esophageal Cancer Hospital (Henan Province, China) from 2018 to 2019. The various lesions were diagnosed independently by at least two pathologists according to the American Joint Committee on Cancer Eighth edition. No patient had received chemotherapy or radiotherapy before biopsy or surgery. This study was approved by the Institutional Review Boards of Cancer Hospital, Chinese Academy of Medical Sciences and informed consent was obtained from each patient. Clinical information was collected from patients’ medical records.
2018 年至 2019 年，在河南省林州市食管癌医院（中国）进行手术或内镜检查时，收集了 4 名患者的 ESCC 肿瘤、异型增生病变和距肿瘤>5 厘米的正常组织（用于激光捕获显微切割，LCM），以及 45 名不同病理阶段的食管病变组织（用于批量 RNA 测序）。各种病变由至少两位病理学家根据美国癌症联合委员会第八版独立诊断。在活检或手术前，患者未接受化疗或放疗。本研究经中国医学科学院肿瘤医院机构审查委员会批准，并从每位患者处获得了知情同意。临床信息从患者的医疗记录中收集。

Induction of multi-staged ESCC development and sample preparations
诱导多阶段食管鳞状细胞癌的发展及样本准备

Animal experiments in this study were conducted in compliance with approved protocols and guidelines from the Institutional Animal Care and Use Committee of the Chinese Academy of Medical Sciences. Eight-week-old female C57BL/6 mice, purchased from the Beijing Huafukang bioscience company in China, were maintained in local housing facility of a controlled condition (23 ± 1 °C, 50 ± 10% humidity and 12–12 h light-dark cycle). Mice were treated with 4NQO (Sigma-Aldrich) in drinking water (100 μg/ml) for 16 weeks to induce multi-staged ESCC carcinogenesis. Drinking water containing the carcinogen was replaced once a week with freshly prepared one and mice were allowed to access drinking water ad libitum during treatment. After 16 weeks of carcinogen treatment, 4NQO drinking water was replaced by sterile water until the mice were killed.
本研究中的动物实验符合中国医学科学院实验动物护理和使用委员会批准的协议和指南。使用中国北京华阜康生物公司提供的 8 周龄雌性 C57BL/6 小鼠，在受控条件下（温度 23±1℃，湿度 50±10%，12-12 小时光照-黑暗周期）饲养。小鼠通过饮用含 4NQO（Sigma-Aldrich，100μg/ml）的水诱导多阶段食管鳞状细胞癌（ESCC）发生，持续 16 周。每周更换一次含有致癌物的饮用水，并允许小鼠自由饮用。16 周致癌物处理后，用无菌水替代 4NQO 饮用水，直至处死小鼠。

The esophageal lesions were identified by two independent pathologists based on the histopathological criteria described previously⁵⁸ (Fig. 1a). Briefly, stage NOR was well-oriented stratified epithelium consists of basal zone and superficial zone. Stage INF was normal epithelium with focal aggregates of epithelial lymphocytes. Stage DYS was defined as loss of polarity in the epithelial cells, nuclear pleomorphism, hyperchromatic, and increased or abnormal mitoses. In stage HYP, these abnormalities were confined to the lower third of the epithelium while in DYS they present in lower two thirds of the epithelium. Lesions with such abnormal changes involving the entire thickness of epithelium were considered as carcinoma in situ (stage CIS). Stage ICA was defined as a lesion with invasion into the sub-epithelial tissues. Esophageal samples of 4NQO-induced mice were subjected to single-cell RNA sequencing at six different time points and indicated different pathological stages: stage NOR (at week 0), stage INF (at week 12), stage HYP at (week 20), stage DYS (at week 22), stage CIS (at week 24), and stage ICA (at week 26). A group of control mice treated without 4NQO were killed at month 1, 2, 8, and 25 (n = 2, respectively). The esophagus was removed immediately when the animal was killed. Cross-sections of the esophagus were cut and stored at –80 °C and sections of the frozen tissues were stained with hematoxylin–eosin (H&E) for histopathological examination and microdissection for bulk RNA sequencing.
食管病变根据先前描述的组织病理学标准由两位独立的病理科医师识别 ⁵⁸ (图 1a )。简而言之，NOR 阶段为有基底区和表层区的定向良好分层上皮。INF 阶段为正常上皮，伴有局部上皮淋巴细胞聚集。DYS 阶段定义为上皮细胞失去极性，核多形性，染色深，以及增多或异常的有丝分裂。在 HYP 阶段，这些异常局限于上皮的下三分之一，而在 DYS 阶段则出现在下三分之二。涉及整个上皮厚度的异常变化被认为是原位癌（CIS 阶段）。ICA 阶段定义为侵入上皮下组织的病变。诱导食管癌的 4NQO 小鼠在六个不同时间点进行了单细胞 RNA 测序，并表明不同的病理阶段：NOR 阶段（第 0 周），INF 阶段（第 12 周），HYP 阶段（第 20 周），DYS 阶段（第 22 周），CIS 阶段（第 24 周），ICA 阶段（第 26 周）。 一组未接受 4NQO 处理的对照小鼠在第 1、2、8 和 25 个月时被处死（每组分别有 2 只）。动物被处死后立即取出食管。食管横截面被切成薄片并保存在-80 °C 下，冷冻组织的切片用苏木精-伊红（HE）染色进行病理学检查和 bulk RNA 测序的宏分离。

Bulk RNA sequencing and data analysis
批量 RNA 测序和数据分析

Tissue samples from human or mice were cut into 5–10 consecutive sections (8 μm) and the epithelial layer contained 30–50 cells on each section was micro-dissected with a Leica LMD7000 laser-capture microdissection system. RNA was isolated from the mini-bulk samples and cDNA was prepared for sequencing based on the Geo-seq protocol⁵⁹. Sequencing libraries were built using the TruePrep DNA Library Prep Kit V2 for Illumina (Vazyme), and evaluated by Bioanalyzer (DNA HS kit, Agilent). RNA-seq data were mapped to GRCh38 human genome and GRCm38 murine genome by HISAT2 (version 2.1.0)⁶⁰ with default parameter for human and murine samples, respectively. The gene expression matrix of raw reads counts after annotation by HTSeq (version 0.6.1p1) was processed using the DESeq2 (version 1.22.2)⁶¹ and visualized by showing the first 3 dimensions calculated by plotPCA function. We used TRIZOL to extract bulk RNA from patient biopsy (n = 45), and constructed sequencing library using NEBNext Ultra II RNA Library Prep Kit for Illumina. Sequencing data were processed by HISAT2 and HTseq as described above. The normalized expression from bulk samples of human precancerous lesions was used to estimated neutrophil fraction with CIBERSORT (version 1.06)⁶².
人或小鼠的组织样本被切成 5–10 个连续切片（8 μm），每个切片中的上皮层包含 30–50 个细胞，使用 Leica LMD7000 激光捕获显微切割系统进行微切割。从 mini-bulk 样本中提取 RNA，并根据 Geo-seq 协议制备 cDNA 进行测序 ⁵⁹ 。使用 Vazyme 的 TruePrep DNA Library Prep Kit V2 for Illumina 构建测序文库，并通过 Bioanalyzer（DNA HS kit，Agilent）进行评估。RNA-seq 数据使用 HISAT2（版本 2.1.0） ⁶⁰ 分别对人和小鼠样本使用默认参数进行比对，映射到 GRCh38 人基因组和 GRCm38 小鼠基因组。使用 HTSeq（版本 0.6.1p1）对注释后的原始读数计数生成基因表达矩阵，并使用 DESeq2（版本 1.22.2） ⁶¹ 处理，通过 plotPCA 函数计算的前 3 个维度进行可视化。我们使用 TRIZOL 从患者活检样本（n = 45）中提取 bulk RNA，并使用 NEBNext Ultra II RNA Library Prep Kit for Illumina 构建测序文库。测序数据的处理方法如上所述。 从 bulk 样本中的人癌前病变的归一化表达量使用 CIBERSORT（版本 1.06）估算中性粒细胞比例 ⁶² 。

Single-cell RNA sequencing (scRNA-seq) and data analysis
单细胞 RNA 测序（scRNA-seq）和数据分析

For mice at stage NOR and INF, the whole esophagi were taken and for mice at other stages, the dysplastic/malignant lesions were taken immediately after killed. The tissue samples of ESCC and various precursor lesions were gently minced into small pieces and digested for with in RPMI-1640 medium (Invitrogen) containing collagenase IV (Gibco) and hyaluronidase (Sigma-Aldrich). CD45-FITC (553080, BD Biosciences, dilution 1:20) antibody staining was performed for fluorescence activated cell sorting (FACS) on a FACSAria sorter (BD Biosciences). Single cells with or without GFP signal, representing immune or non-immune cells, were sorted and captured respectively in nanoliter droplets using Chromium (10× Genomics). scRNA-seq libraries were prepared using Chromium Single Cell 5′ Reagent Kits (10× Genomics) and sequencing was accomplished with an Illumina HiSeq x10 System.
对于 NOR 和 INF 阶段的小鼠，整个食管被取下；而对于其他阶段的小鼠，则在处死后立即取下异常增生/恶性病变组织。ESCC 及其各种前驱病变组织样本被轻轻剪碎成小块，并用含胶原酶 IV（Gibco）和透明质酸酶（Sigma-Aldrich）的 RPMI-1640 培养基（Invitrogen）消化。使用 BD Biosciences 的 CD45-FITC（553080，稀释度 1:20）抗体进行荧光激活细胞排序（FACS）染色，使用 BD Biosciences 的 FACSAria 分选器。带有或不带有 GFP 信号的单个细胞，代表免疫或非免疫细胞，分别被分选并用 Chromium（10× Genomics）捕获在纳升液滴中。使用 10× Genomics 单细胞 5′试剂盒制备 scRNA-seq 文库，并使用 Illumina HiSeq x10 系统进行测序。

Raw gene expression matrices obtained per sample using CellRanger (version 2.1.0, 10× genomics) were combined using the Seurat R package (version 2.3.4)⁶³. Genes detected in <0.1% of all cells were filtered. We further excluded cells with gene counts <500 and cells that had >10% of mitochondrial gene expressions. After quality control, 66,089 cells were further analyzed for their gene expression profiles. The genes with normalized expression between 0.0125 and 3, and dispersion >0.5 were selected as highly variable genes. The resultants were first summarized by principle component analysis (PCA) and then first several PCs were selected for tSNE dimensional reduction using the default settings of the RunTSNE function. The numbers of resulting highly variable genes and the select PCs are shown in Supplementary Table 2. Cell clusters in the resulting two-dimensional representation were annotated as known biological cell types using canonical marker genes.
每个样本使用 CellRanger（版本 2.1.0，10× genomics）获得原始基因表达矩阵，并使用 Seurat R 包（版本 2.3.4）进行合并 ⁶³ 。过滤掉在所有细胞中检测不到<0.1%的基因。进一步排除基因计数<500 的细胞以及线粒体基因表达>10%的细胞。经过质量控制后，66,089 个细胞进一步分析其基因表达谱。选择归一化表达在 0.0125 到 3 之间且离散度>0.5 的基因作为高度可变基因。首先通过主成分分析（PCA）汇总结果，然后使用 RunTSNE 函数的默认设置进行 tSNE 维度降维，选择前几个主成分。高度可变基因的数量和选择的主成分数量显示在补充表 2 中。在结果的二维表示中，细胞簇被注释为已知的生物细胞类型，使用经典标记基因。

Major cell-type clustering and marker gene identification
主要细胞类型聚类和标志基因识别

We reanalyzed epithelial cells and stromal cells separately to identify their sub-clusters by using highly variable gene identification and dimensional reduction as described above (Supplementary Table 2). Cells with mix features were removed from further analysis (e.g., cells with both Cd3d and Cd19 expression indicating T cell and B cell Multiplet). Clusters were identified using FindClusters function, and the specific gene markers for each cluster were determined using the FindAllMarkers function implanted in Seurat package.
我们分别重新分析了上皮细胞和间质细胞，通过上述方法（见补充表 2 ）识别高度可变基因并进行维度减少，以识别其亚簇。具有混合特征的细胞被从进一步分析中移除（例如，同时表达 Cd3d 和 Cd19 的细胞表明是 T 细胞和 B 细胞的 Multiplet）。使用 FindClusters 函数识别簇，并使用 Seurat 包内置的 FindAllMarkers 函数确定每个簇的特定基因标记。

Gene set variation analysis (GSVA)
基因集变异分析（GSVA）

Pathway analyses were predominantly performed on the 50 hallmark pathways described in the molecular signature database, exported using the MSigDB database (version 6.2)⁶⁴. We also assessed biological process activities using a described biological process of Gene Ontology (GO) dataset. To assign pathway activity estimates to individual cells, we applied the GSVA using standard settings, as implemented in the GSVA package (version 1.30.0)⁶⁵. To assess differential activities of pathways between sub-cluster of cells, we contrasted the activity scores for each cell using Limma package (version 3.38.3)⁶⁶. Differential activities of pathways were calculated for each identified cluster. T-values of the results of some significant differential pathways (P < 0.05) in top 10 were visualized using heatmaps with average pathway activity scores of each cluster.
通路分析主要针对分子签名数据库（MSigDB 数据库，版本 6.2）中描述的 50 个标志性通路进行。我们还使用基因本体（GO）数据集描述的生物学过程评估了生物学过程活性。为了将通路活性估计值分配给单个细胞，我们使用 GSVA 包（版本 1.30.0）的标准设置进行了 GSVA 分析 ⁶⁵ 。为了评估细胞亚簇之间通路活性的差异，我们使用 Limma 包（版本 3.38.3）对比了每个细胞的活性评分 ⁶⁶ 。对于每个识别出的簇，计算了通路活性的差异。对于一些显著差异的通路（P < 0.05），在前 10 个中，使用热图可视化每个簇的平均通路活性评分。

Analysis of transcription factor expression
转录因子表达分析

SCENIC (version 1.1.0)⁶⁷ was used to assess the transcriptional activity of epithelial cells with high quality (UMI > 5500). The analysis used the motifs database for RcisTarget and GRNboost (corresponding to GENIE3 1.4.3, AUCell 1.4.1 and RcisTarget 1.2.1; with mm10__refseq-r80__10kb_up_and_down_tss.mc9nr). The input matrix was read counts.
SCENIC（版本 1.1.0） ⁶⁷ 用于评估具有高质量（UMI > 5500）的上皮细胞的转录活性。分析使用了 RcisTarget 和 GRNboost 的 motifs 数据库（对应于 GENIE3 1.4.3、AUCell 1.4.1 和 RcisTarget 1.2.1；使用 mm10__refseq-r80__10kb_up_and_down_tss.mc9nr）。输入矩阵为读取计数。

Cell transition trajectory and diffusion map analysis
细胞转换轨迹和扩散图分析

Monocle 2 (version 2.10.1)⁶⁸ was used for the trajectory analysis on high quality epithelial cells (UMI > 5500). All the top 100 markers of each cluster were used for the cell ordering. Dimensionality reduction and trajectory construction were performed on the selected genes with default methods and parameters. We calculated diffusion components using the RunDiffusion function as implemented in Seurat package with default parameters. The first three dimensions were used to draw diffusion maps. Mean coordinates of all of each cluster’s cells were considered as the center of the cluster. The farthest cell of stage NOR from the total distance of other stages was the start point.
Monocle 2（版本 2.10.1） ⁶⁸ 用于高质量上皮细胞（UMI > 5500）的轨迹分析。每个簇的前 100 个标记基因用于细胞排序。在选定基因上进行降维和轨迹构建，默认方法和参数。我们使用 Seurat 包中的 RunDiffusion 函数计算扩散成分，默认参数。使用前三个维度绘制扩散图。每个簇所有细胞的平均坐标被视为簇的中心。从总距离来看，阶段 NOR 中最远的细胞作为起点。

Analysis of interaction between cell types
细胞类型间的相互作用分析

We used CellPhoneDB (version 2.0.6)⁶⁹ with default arguments to reveal interaction between cell types. Each cluster analyzed was downsampled to 100 cells since the low cell number of some epithelial cluster. Interactions of epithelial cells with immune cells and immune cells with epithelial cells were demonstrated respectively.
我们使用了 CellPhoneDB（版本 2.0.6） ⁶⁹ 并采用默认参数来揭示不同细胞类型之间的相互作用。由于某些上皮细胞簇的细胞数量较低，每个分析的簇被下采样至 100 个细胞。上皮细胞与免疫细胞之间的相互作用以及免疫细胞与上皮细胞之间的相互作用分别得到了展示。

Immunohistochemistry and immunofluorescent detection
免疫组化和免疫荧光检测

Formalin-fixed paraffin-embedded (FFPE) sections of esophageal precursor lesions were collected from 30 patients between 2016 and 2018 in Linzhou Cancer Hospital, including INF (n = 5), HYP (n = 13), DYS (n = 8) and CIS (n = 3), to validate the results obtained in mice. The protein expression levels of the marker genes were detected by IHC staining for mice tissues and immunofluorescence for human specimens with antibodies (Abcam) shown in Supplementary Table 3. The samples were incubated with antibody against Ki67 (1:50 for IHC, ab16667), Top2a (1:8000 for IHC, 1:10,000 for IF, ab52934), Aldh3a1 (1:200 for IHC, 1:600 for IF, ab76976), Atf3 (1:200 for IHC, 1:600 for IF, ab216569), S100a8 (1:500 for IHC, 1:1500 for IF, ab92331), Mmp14 (1:2000 for IHC, 1:6000 for IF, ab51074), or Itga6 (1:250 for IHC, 1:750 for IF, ab181551). Opal multiplex staining was performed according to the Opal 5-Color Manual IHC Kit (Perkin Elmer). Opal DAPI, Opal 520, Opal 570, Opal 620, and Opal 690 were used to generate different signals. Slides were counterstained with DAPI (1:2000) for nuclei visualization, and subsequently coverslipped using a VectaShield Hardset mounting media. The slides were imaged using Vectra Polaris Automated Quantitative Pathology Imaging System (Perkin Elmer). We used inForm software (Perkin Elmer) to unmix and remove autofluorescence and to analyze the multispectral images.
食管前驱病变的甲醛固定石蜡包埋（FFPE）切片来自 2016 年至 2018 年林州癌症医院的 30 名患者，包括 INF（n=5）、HYP（n=13）、DYS（n=8）和 CIS（n=3），用于验证在小鼠中获得的结果。通过 IHC 染色检测小鼠组织中标志基因的蛋白表达水平，并通过免疫荧光检测人类标本，使用抗体（Abcam），详见补充表 3 。样本用抗体孵育，针对 Ki67（IHC 1:50，ab16667），Top2a（IHC 1:8000，IF 1:10,000，ab52934），Aldh3a1（IHC 1:200，IF 1:600，ab76976），Atf3（IHC 1:200，IF 1:600，ab216569），S100a8（IHC 1:500，IF 1:1500，ab92331），Mmp14（IHC 1:2000，IF 1:6000，ab51074），或 Itga6（IHC 1:250，IF 1:750，ab181551）。根据 Perkin Elmer 的 Opal 5-Color Manual IHC Kit 进行 Opal 多色染色。使用 Opal DAPI、Opal 520、Opal 570、Opal 620 和 Opal 690 生成不同信号。切片用 DAPI（1:2000）复染以观察核，随后使用 VectaShield Hardset 封片介质封片。 切片使用 Perkin Elmer 的 Vectra Polaris 自动定量病理成像系统成像。我们使用 Perkin Elmer 的 inForm 软件进行去混和去除自发荧光，并分析多光谱图像。

Statistical analysis  统计分析

Statistical analyses were conducted by using R v3.5.1⁷⁰ and Prism 7 (Graphpad Software). Pearson’s correlation was calculated with the R function cor() and the significance was determined using two-sided unpaired Wilcoxon rank-sum test. P < 0.05 was considered statistically significant.
统计分析使用了 R v3.5.1 和 Prism 7（Graphpad Software）。使用 R 函数 cor()计算了皮尔森相关系数，并使用双侧非配对 Wilcoxon 秩和检验确定了显著性。P < 0.05 被认为具有统计学意义。

Reporting summary  报告摘要

Further information on research design is available in the Nature Research Reporting Summary linked to this article.
进一步的研究设计信息请参见本文链接的 Nature Research Reporting Summary 。