Trajectory of Lymph Node Metastasis Reflects the Molecular Features of Esophageal Squamous Cell Carcinoma

Jiho Park; Seong Yong Park; Ha Eun Kim; Se-Young Jo; Jeongsoo Won; Dae Joon Kim; Sangwoo Kim

doi:10.4143/crt.2025.257

Articles

Page Path: HOME > J Korean Cancer Assoc > Ahead-of print articles > Article

Original Article Trajectory of Lymph Node Metastasis Reflects the Molecular Features of Esophageal Squamous Cell Carcinoma: Jiho Park¹, Seong Yong Park^2,a), Ha Eun Kim², Se-Young Jo¹, Jeongsoo Won¹, Dae Joon Kim², Sangwoo Kim¹; DOI: https://doi.org/10.4143/crt.2025.257
Published online: October 15, 2025

¹Department of Biomedical Systems Informatics, Yonsei University College of Medicine, Seoul, Korea

²Department of Thoracic and Cardiovascular Surgery, Yonsei University College of Medicine, Seoul, Korea

Correspondence: Dae Joon Kim, Department of Thoracic and Cardiovascular Surgery, Yonsei University College of Medicine, 50-1 Yonsei-ro, Seodaemun-gu, Seoul 03722, Korea
Tel: 82-2-3410-1696 E-mail: kdjcool@yuhs.ac

Co-correspondence: Sangwoo Kim, Department of Biomedical Systems Informatics, Yonsei University College of Medicine, 50-1 Yonsei-ro, Seodaemun-gu, Seoul 03722, Korea
Tel: 82-2-2228-2589 E-mail: swkim@yuhs.ac

^*Jiho Park and Seong Yong Park contributed equally to this work.
^a)Present address: Department of Thoracic and Cardiovascular Surgery, Samsung Medical Center, Sungkyunkwan University School of Medicine, Suwon, Korea

• Received: March 5, 2025 • Accepted: October 14, 2025

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

2,371 Views
91 Download

Full Article

Download PDF

Abstract

Purpose
Esophageal squamous cell carcinoma (ESCC) is frequently accompanied by lymph node metastasis (LNM) to the neck, chest, and abdomen. Despite its significance as a key prognostic factor, the genomic trajectory of LNM remains poorly understood. This study aimed to characterize the underlying patterns and genomic characteristics of LNM.
Materials and Methods
Whole-exome sequencing and transcriptome sequencing were performed on 45 multiregional samples (10 primary tumors, 10 normal esophageal tissues, and 25 lymph node tumors) from 10 ESCC patients who underwent esophagectomy with three-field lymphadenectomy. The temporal trajectory of metastasis was reconstructed through phylogenetic analysis, leveraging somatic mutations identified in the primary tumor and lymph nodes.
Results
Somatic mutations preceding metastasis included major driver mutations, such as TP53 and KMT2D, and displayed a mutational process associated with alcohol consumption (SBS16), emphasizing its influence on early tumorigenesis. In contrast, post–lymph node metastatic mutations were sporadic. Lymph nodes seeded later acquired mutations at a faster rate, suggesting increased genomic instability. In three of nine patients (33.3%), nodal skip metastasis (NSM) was observed, including two cases detected exclusively via genomic analysis, highlighting the necessity of phylogenetic assessment to avoid misclassification. Transcriptome analysis revealed activation of epithelial-mesenchymal transition and KRAS signaling pathways in NSM tumors, indicative of poor prognostic outcomes.
Conclusion
Our study provides a molecular understanding of LNM, emphasizes the potential importance of node-skipping patterns in ESCC, and underscores the utility of genomic analysis in elucidating the connection between LNM.
Key words: Phylogeny, Esophageal squamous cell carcinoma, Nodal skip metastasis, Lymphatic metastasis

Introduction

Esophageal squamous cell carcinoma (ESCC) is the predominant subtype of esophageal cancer, accounting for approximately 90% of cases in East Asia [1]. The 5-year survival rate is low, typically below 25% [2], primarily due to its frequent late-stage diagnosis and the high incidence of lymph node metastasis (LNM). The complex lymphatic system of the esophagus, encompassing longitudinal lymphatic vessels within the submucosa and lamina propria, complicates the interpretation and management of LNM. The multidirectional propagation of tumor cells further contributes to the poor prognosis of ESCC.

Nodal skip metastasis (NSM) is a phenomenon in which disseminated tumor cells bypass proximal lymph nodes to metastasize directly to more distal nodes. NSM has been reported not only in ESCC but also in other cancers, including breast cancer [3], oral cavity squamous cell carcinoma [4], and non–small cell lung cancer. Reported frequencies of NSM in ESCC vary considerably, from 18.2% to 64.0%, and its prognostic significance remains controversial [5-7]. Conventional NSM classification, based solely on histopathological assessment, does not adequately reflect the intricate genetic relationships among lymph nodes and between lymph nodes and the primary tumor.

Phylogenetic analysis provides a powerful approach for reconstructing the genetic evolution of cancer [8]. By examining somatic mutation patterns across primary tumors and metastatic sites, it enables precise delineation of tumor progression and the sequence of metastatic events. This framework facilitates a deeper understanding of the genetic relationships between primary tumors and metastatic lymph nodes, offering critical insights into LNM. Despite its potential, phylogenetic analysis has rarely been applied to ESCC, where it could substantially enhance our understanding of metastatic dynamics.

In this study, we conducted phylogenetic analysis on 45 multiregional tissue samples, comprising paired primary tumors and metastatic lymph nodes, from 10 ESCC patients to elucidate the genomic and transcriptomic landscape of LNM. Reconstruction of metastatic trajectories across multiple lymph nodes revealed temporal and site-specific molecular features, such as mutational rates, etiological factors, NSM occurrence, and activation of specific oncogenic pathways. This study aims to provide novel perspectives on LNM patterns and improve our understanding of the molecular characteristics of metastatic progression in ESCC.

Materials and Methods

1. Patients and sample preparation

Ten Korean patients diagnosed with ESCC who underwent esophagectomy with three-field lymphadenectomy were enrolled. Eligibility criteria included: (1) no prior neoadjuvant therapy, (2) presence of multiple LNM in the final pathology report, and (3) no history of other cancers.

Across the cohort, a total of 26 to 92 lymph nodes (mean, 69.9) were resected per patient. Of these, 2 to 6 lymph nodes (mean, 4.6) were histologically confirmed to harbor metastases, while the remaining nodes were negative for tumor involvement. All tumor-positive nodes meeting technical criteria were sequenced; the SureSelect Human Exon V6 kit requires ≥ 0.200 μg double-stranded genomic DNA, and lesser amounts were excluded as insufficient. Detailed information on the number and anatomical distribution of metastatic lymph nodes is provided in S1 Table. S2 Table lists the technical reason sequencing was not performed for each lymph node.

Pathological staging was determined according to the American Joint Committee on Cancer/Union for International Cancer Control (AJCC/UICC) 7th edition [9], and lymph node stations were classified using the Japanese Classification of Esophageal Cancer (JES), 11th edition [10].

2. Whole-exome sequencing

The Agilent SureSelect Target Enrichment protocol for the Illumina paired-end sequencing library (version C2, December 2018) was used for whole-exome sequencing (WES) with 1 μg of input genomic DNA and 200 ng of input formalin-fixed paraffin-embedded DNA. The SureSelect Human All Exon V6 probe set was used for target enrichment. PicoGreen and agarose gel electrophoresis were employed to measure the quantity and quality of DNA. We used 1 μg of each tissue’s genomic DNA and 200 ng of formalin-fixed paraffin-embedded DNA diluted in EB buffer and sheared to a target peak size of 150-200 bp using the Covaris LE220 focused-ultrasonicator (Covaris) according to the manufacturer’s protocols. The captured DNA was then washed and amplified. The final purified product was quantified using qPCR following the qPCR Quantification Protocol Guide (KAPA Library Quantification Kits for Illumina Sequencing Platforms) and qualified using the TapeStation DNA screen tape D1000 (Agilent Technologies). Sequencing was performed using the HiSeq 2500 platform (Illumina).

3. Somatic variant calling and filtration

Sequence reads were aligned to the human reference genome (hg38) using BWA-MEM (v. 0.7.17) with default parameters (https://github.com/lh3/bwa). Post-alignment processing, including MarkDuplicates and FixMateInformation, was performed with GATK (v. 4.2.3.0) (https://gatk.broadinstitute.org/hc/en-us/). Somatic single-nucleotide variants (sSNVs) were identified using GATK Mutect2 (v. 4.2.3.0) for each tumor-normal pair, followed by variant filtering with FilterMutectCalls (v. 4.2.3.0). Small insertions and deletions (indels) were detected using Strelka2 (v. 2.9.10) (https://github.com/Illumina/strelka/). After sSNV and indel calling, two lymph nodes (101L in P3 and 1 in P8) were excluded due to the absence of shared variants with their corresponding primary tumors.

A post-variant rescue method (“force calling”) [11] was applied, in which sites initially classified as non-variants were re-evaluated and reclassified as variants if detected in other samples from the same patient and meeting all of the following criteria: (1) read depth ≥ 10, (2) number of mutant alleles ≥ 3 in the tumor sample, (3) base quality ≥ 20, (4) mapping quality ≥ 30, and (5) number of mutant alleles ≤ 2 in the paired normal sample. Indels were manually inspected using Integrative Genomics Viewer (https://github.com/igvteam/igv/). For indels detected across multiple tumor samples from a patient, those sharing the same indel form and absent in the paired normal sample were rescued if previously filtered out. All sSNVs and indels, including rescued variants, were annotated with ANNOVAR (v. 20191024) (http://annovar.openbioinformatics.org/). Tumor mutation burden (TMB) was defined as the number of nonsynonymous mutations per megabase (Mb) of the examined genome.

4. RNA sequencing and data processing

Transcriptome sequencing was performed on seven primary tumors (excluding patients P2, P6, and P7) and their matched normal mucosal tissues. Total RNA concentration was measured using Quant-IT RiboGreen (#R11490, Invitrogen). To assess RNA integrity, samples were analyzed on the TapeStation RNA ScreenTape (#5067-5576, Agilent Technologies), and those with RNA integrity number scores below 7.0 were excluded from RNA library construction. Libraries were independently prepared with 1 μg of total RNA per sample using the Illumina TruSeq Stranded mRNA Sample Prep Kit (#RS-122-2101, Illumina). The products were purified and enriched by PCR to create the final cDNA library. The libraries were quantified using the KAPA Library Quantification Kits for Illumina Sequencing Platforms, following the qPCR Quantification Protocol Guide (#KK4854, KAPA Biosystems) and qualified using the TapeStation D1000 ScreenTape (#5067-5582, Agilent Technologies). Indexed libraries were then subjected to paired-end (2×100 bp) sequencing on the Illumina NovaSeq platform (Illumina) at Macrogen.

Sequence reads were aligned to a genome index constructed from the human reference genome (hg38) and GENCODE annotation (v38) using STAR (v. 2.7.3a) (https://github.com/alexdobin/STAR/). After alignment, the primary tumor from P4 and normal mucosa samples from P8 and P10 were excluded because their expression profiles did not match the expected ESCC-specific or normal esophageal reference profiles from GEPIA2 (log-fold change > 2; p < 0.05) (S3 Fig.) (http://gepia2.cancer-pku.cn/). For gene expression quantification, raw counts were extracted using featureCounts (v. 2.0.1) (https://subread.sourceforge.net/) and normalized via variance-stabilizing transformation (VST) with DESeq2 (v. 1.30.1) (https://www.bioconductor.org/packages/release/bioc/html/DESeq2.html/).

5. Judging the functional impact of mutations

The functional impact of somatic variants was evaluated based on two criteria: (1) the presence of mutations in cancer driver genes and (2) assessment of their functional effects on protein function. A reference set of 621 cancer driver genes was compiled by integrating three databases: 299 The Cancer Genome Atlas (TCGA) pan-cancer driver genes [12], 483 “tier 1” COSMIC cancer gene census (https://cancer.sanger.ac.uk/census/), and 40 esophageal cancer-specific driver genes from IntoGen (https://www.intogen.org/search/).

The functional impact of sSNVs was predicted using SnpEff (v. 5.0) (https://pcingola.github.io/SnpEff/snpeff/introduction/) and categorized as low, moderate, and high. High-impact variants were considered functionally damaging. Moderate-impact variants were considered as functionally damaging if at least two of the following criteria were met: predicted as “deleterious” by SIFT [13], “deleterious” by PROVEAN [14], or “functional” by REVEL [15]. Variants with low impact were excluded. All indels were classified as functionally damaging. Functionally damaging variants were considered as potential drivers of tumorigenesis, and the genes harboring these variants were designated as cancer driver genes.

6. Phylogenetic tree construction and analysis

Phylogenetic relationships among multiple samples (normal, primary, and LNM) from each patient were inferred from sSNVs, including synonymous variants. Phylogenetic trees were reconstructed using the maximum parsimony method implemented in the PHYLogeny Inference Package v.3.698 (PHYLIP) (https://phylipweb.github.io/phylip/), with the matched normal sample designated as the outgroup root. Trees were manually visualized to ensure that all branch lengths were proportional to the number of somatic mutations. Somatic mutations were then classified by temporal occurrence: those acquired between initial tumorigenesis (t₀) and the first metastatic event (t₁) were defined as pre–lymph node metastatic (pre-LN-metastatic), while those acquired after t₁ were classified as post–lymph node metastatic (post-LN-metastatic).

For each LNM (leaf node), two measures were derived from branch lengths: the timing of metastasis (ToM) and mutation accumulation rate after metastasis (MARM). ToM represents the relative timing of LNM during tumor progression and was calculated as the ratio of the primary tumor branch length from initiation (t₀) to the time of the metastatic event over the total primary tumor branch length from initiation to resection. MARM represents the rate of de novo somatic mutation accumulation in the lymph nodes after metastasis and was calculated as the ratio of the private branch length of the lymph node to the length of the primary tumor branch from the time of the metastatic event to resection. An overview is provided in S4A Fig. The association between ToM and MARM was evaluated using Spearman’s correlation coefficient.

7. Mutational signature analysis

Mutational signatures were identified using DeconstructSigs (v. 1.9.0) (https://github.com/raerose01/deconstructSigs/). All SNVs, including synonymous mutations, were included as input. COSMIC mutational signature data (signatures.exome.cosmic.v3.may2019; v. 3.0) was applied to estimate the relative contribution of each signature and reconstruct the mutational profiles. Patient-level comparisons in signature enrichment were assessed using the Mann-Whitney U test (Wilcoxon rank-sum test).

8. Identification of NSM and gene set enrichment analysis

All lymph nodes were classified into three groups (N1, N2, and N3) according to the JES classification [10]. Patients with sequential metastasis from N1 to N2 to N3 were categorized as non-NSM. Patients deviating from this order, metastasizing to N2/N3 before N1 or bypassing N1 entirely, were labeled as NSM. One patient with a single LNM (P2) was excluded from the analysis.

Gene set enrichment analysis (GSEA v. 4.2.3) (https://www.gsea-msigdb.org/) was performed using the Hallmark gene set database, with permutation type set to “gene set” and the ranking metric set to “Diff of Classes.” Three pairwise comparisons were tested (NSM vs. non-NSM, NSM vs. normal, and non-NSM vs. normal). Gene sets with a false discovery rate (FDR) < 0.005 were considered significantly enriched.

9. Comparative analysis of epithelial-mesenchymal transition markers

A list of the top 10 epithelial-mesenchymal transition (EMT) markers was obtained from the EMTome resource (http://www.emtome.org/), comprising one epithelial marker (CDH1) and nine mesenchymal markers (CDH2, FN1, MMP2, SNAI1, SNAI2, SPARC, VIM, ZEB1, ZEB2). EMTome ranks these genes based on their prevalence across 445 publications related to EMT and mesenchymal-epithelial transition. The VST values were normalized for each gene using min-max scaling, and unsupervised hierarchical clustering was performed on the 10 EMT markers across 11 samples (five normal, two NSM, and four non-NSM).

10. Statistical analysis

Statistical analyses were performed in R (v. 4.2.0) (https://cran.r-project.org/bin/windows/base/). The Mann-Whitney U test was applied to compare two independent groups. Spearman’s rank correlation was used to assess associations between variables.

Results

1. Patient cohort and characteristics

Ten patients (P1-P10) with ESCC and multiple LNM were enrolled in this study. From each patient, tissues from the primary tumor, adjacent normal mucosa, and one to three LNM were obtained for WES. In total, 45 multiregional samples were collected: 10 primary tumors, 10 normal mucosae, and 25 metastatic lymph nodes (Fig. 1A). Clinical characteristics, including age and sex, were evenly distributed among patients (Table 1). LNM distribution varied by primary tumor location, consistent with a recent systematic review highlighting differences in lymph node metastatic spread depending on the primary tumor site. Patients with primary tumors in the upper thoracic region (P1, P5, and P8) predominantly exhibited LNM in the upper thoracic esophagus, whereas those with middle thoracic tumors (P2, P3, P4, P6, P7, P9, and P10) showed LNM distributed throughout the esophagus, rather than confined to a specific region.

2. Somatic mutations in primary tumors and LNM

WES of 35 tumor tissues (10 primary tumors and 25 LNM) and matched controls (10 normal mucosae) identified 5,580 somatic SNVs and 267 indels. The median TMB was 1.55 mutations/Mb, within the previously reported range (0.840-5.60 mutations/Mb) observed in ESCC [16,17]. Somatic mutation counts were comparable between primary tumors (99-258; mean, 165.40) and metastases (46-324; mean, 167.72).

Somatic mutation profiling revealed recurrent alterations in cancer driver genes, underscoring their contribution in ESCC progression (Fig. 1B). All primary tumors (10/10, 100%) and matching LNM (25/25, 100%) harbored damaging or truncating mutations in TP53, indicating its essential role in tumorigenesis. In addition, we identified three other frequently mutated genes (DNAH5, KMT2D, and SYNE1) that may act as tumor suppressors in ESCC. DNAH5 mutations have previously been associated with poor overall survival in ESCC [18]. KMT2D has been identified as a tumor suppressor in multiple tumor types, including ESCC [19]. Hypermethylation of the promoter of SYNE1 in gastric cancer [20] and somatic mutations in renal cell carcinoma [21] have been associated with tumor aggressiveness and poorer outcomes. Notably, the de novo acquisition of these mutations in LNM (SYNE1 in P9 and DNAH5 in P1 and P7) suggests a selective advantage during tumor progression.

3. Trajectory of LNM

Cancer development can be modeled as an evolutionary process [8], wherein progressive events, such as initial transformation (at time t₀) and subsequent spread to N lymph nodes (at time t₁,…, t_N > t₀), are mapped to an ordered pseudo-time (Fig. 2A, left). The endpoint (t_E > t₀,…, t_N) denotes the time of observation (resection), corresponding to tumor status at resection. Throughout the trajectory, somatic mutations accumulate continuously, driving genomic heterogeneity within and across tumors. This information can be leveraged to reconstruct phylogenetic trees that represent the order of metastatic events (Fig. 2A, right).

For classification, mutations arising between initial tumorigenesis (t₀) and the first metastatic event (t₁) were defined as pre-LN-metastatic, whereas those acquired after t₁ were defined as post-LN-metastatic.

Phylogenetic trees of the 10 patients demonstrated diverse histories of cancer progression (Fig. 2B). Most LNM (21/25, 84%) branched directly from primary tumors, except for two lymph node pairs, 104R/106recR in P1 and 104/106recR in P5, which metastasized concurrently, suggesting possible node-to-node spread. The interval between initial tumorigenesis (t₀) and the first metastatic event (t₁) (Fig. 2B, dark blue), measured by the number of somatic mutations at the pre-LN-metastatic stage (t₀-t₁), varied across patients (shortest in P4, P8, and P10, and longest in P9).

Classifying mutations by timing and location clarified their genomic roles in ESCC progression (Fig. 2C). All patients harbored pre-LN-metastatic TP53 mutations, underscoring its crucial role in tumor initiation and early progression. Similarly, three truncating mutations in KMT2D (Gly133fs in P1, Ser2807fs in P5, and Leu2957fs in P6) were acquired prior to LNM, consistent with its established tumor suppressor role in ESCC [19]. Additional recurrent pre-LN-metastatic mutations were identified in NFE2L2 (Ile28Thr in P5 and Glu82Asp in P6) and SPTA1 (Pro2115Thr in P1 and Gln413Lys in P7), suggesting potential involvement of the Nrf2 signaling pathway in tumorigenesis [22]. By contrast, post-LN-metastatic mutations were sparse and heterogeneous, reflecting the underlying genomic diversity of ESCC.

4. Characteristics of somatic mutations in pre-LN and post-LN metastasis

We performed genomic analyses to identify time-specific mutational patterns along the metastasis trajectory. The overall spectrum revealed predominant C > T at NpCpG, C > G at TpCp(A/G), and C > A at (G/T)Cp(A/G), corresponding primarily to SBS1 (clock-like), SBS5 (clock-like), and SBS18 (damage by reactive oxygen species [ROS], respectively (Fig. 3A top, S4B Fig. left), all of which have been previously reported as major mutational signatures in ESCC [23].

Comparison between pre- and post-LN-metastatic mutations revealed stage-specific mutational processes. T > C transitions were more common in pre-LN-metastatic mutations, whereas C > A transversions predominated post-LN-metastatic mutations (Fig. 3A bottom, S4B Fig.). Patient-level analysis identified SBS16 enrichment in pre-LN-metastatic mutations (Mann-Whitney U test, p=0.006). By contrast, SBS5 (p=0.04), SBS6 (defective DNA mismatch repair, p=0.02), and SBS18 (p=0.07) were more prominent in post-LN-metastatic mutations (Fig. 3B).

We further examined whether the occurrence of de novo mutations in metastatic lymph nodes was influenced by the ToM. Later-arising metastases showed a higher MARM, reflecting increased acquisition of de novo mutations (Spearman’s coefficient=0.79, p=3×10^–6) (Fig. 3C, S4C Fig.). Given the progressive decline of genome integrity and repair mechanisms in primary tumors [24], these findings suggest that later-arising metastases may derive from cells with higher genomic instability, potentially contributing to aggressive tumor progression.

5. NSM and molecular characteristics

Of the nine patients, only P3 (1/9, 11.1%) satisfied NSM criteria based on conventional histopathological assessment (Fig. 4A). Interestingly, phylogenetic analysis uncovered two additional patients (P6, P8) with N2 metastasis preceding N1 involvement, a pattern undetected by routine assessment. This raised the NSM prevalence to 33.3% (3/9), underscoring the value of genomic analysis in revealing hidden metastatic trajectories.

Transcriptomic analysis revealed distinct molecular pathway activations between NSM and non-NSM groups (Fig. 4B). In NSM, EMT was significantly enriched (FDR < 2×10^–16), along with angiogenesis (FDR < 2×10^–16), apical junction (FDR < 2×10^–16), and KRAS signaling pathways (FDR=0.001), all of which are directly or indirectly linked to EMT [25]. Gene expression profiles of 10 EMT markers also indicated elevated EMT activity in NSM compared with non-NSM and normal groups (Fig. 4C). Unsupervised hierarchical clustering separated NSM from the other two groups, marked by decreased expression of the epithelial marker (CDH1) and increased expression of nine mesenchymal markers (CDH2, FN1, MMP2, SNAI1, SNAI2, SPARC, VIM, ZEB1, and ZEB2) [26].

In contrast, non-NSM tumors exhibited activation of interferon alpha (FDR < 2×10^–16) and gamma responses (FDR < 2×10^–16) and increased cell proliferation. Enrichment of KRAS signaling down (FDR < 2×10^–16) indicated suppressed KRAS activity. As KRAS activation is known to promote EMT [27], this finding is consistent with the lower EMT activation observed in non-NSM group. Both groups were enriched for general cancer hallmarks relative to normal tissue, including E2F targets, G2M checkpoint, MYC targets, TNFα-signaling via NF-KB, and mTORC1 signaling (S5 Table).

Discussion

This study analyzed the genomic and transcriptomic landscape of ESCC with LNM, highlighting both its molecular characteristics and progression patterns. By employing genomic trajectory reconstruction, we delineated facets of LNM progression beyond what traditional approaches capture. Somatic mutations were classified as pre-LN-metastatic or post-LN-metastatic, with each group exhibiting distinct mutation profiles in frequently mutated genes, functional implications, and mutational processes. We also robustly identified NSM, enabling precise transcriptomic comparisons between NSM and non-NSM groups. Overall, this study underscores the value of multiregional sequencing-based phylogenetic analysis in tracing the sequential progression of LNM and enhancing our understanding of molecular mechanisms in ESCC metastasis.

We demonstrated that tumorigenesis and metastasis are characterized by distinct mutational profiles, indicating divergent genetic underpinnings. The esophagus is particularly vulnerable to carcinogens such as alcohol and tobacco, both well-known risk factors for ESCC [28]. These exposures are closely linked to TP53 mutations, explaining the notably high TP53 mutation frequency in this cancer type [29]. The strong enrichment of SBS16, a mutational signature associated with alcohol consumption, in the pre-LN-metastatic stage supports its role in early tumorigenesis. In contrast, defective mismatch repair and ROS-induced damage, together with the accumulation of sporadic somatic mutations, predominated in the post-LN-metastatic stage. This reflects genomic instability, uncontrolled proliferation, and increased tumor heterogeneity, features that may give rise to the high recurrence rates observed in ESCC.

Phylogenetic analysis enabled identification of sentinel lymph nodes (SLNs) in ESCC. SLNs, defined as the first nodes to receive lymphatic drainage from the primary tumor, are well established in breast cancer [30] and melanoma [31]. In our cohort, the initial LNMs were most commonly detected in the bilateral 106rec and paracardial nodes. Notably, however, phylogenetic evidence revealed that nodes 107 (P6) and 108 (P8) can also serve as initial metastatic sites, broadening the conventional SLN spectrum. This result supports the genetic validity of the SLN concept while emphasizing the importance of radical, three-field lymphadenectomy to account for such exceptions.

Reported frequencies and prognostic significance of NSM remain inconsistent across studies, which we believe it largely because conventional histopathological assessment provides only partial insights into metastatic patterns. This method fails to capture the dynamic nature of metastasis or account for irregular LNM progression (e.g., spread from primary tumors to N2 before N1 nodes) [5-7]. Genetic analyses suggest that the prevalence of NSM has likely been underestimated, contributing to variability in reported frequencies and confounding its association with prognosis. To address this, we applied robust genetic analysis for NSM identification and evaluated its prognostic relevance. In our cohort, survival analysis showed a negative trend toward poorer outcomes in patients with NSM (log-rank test, p=0.4). Although the comparison did not reach statistical significance, likely due to the small sample size and short follow-up period, the data indicate a potential survival difference between the groups.

Between the two ESCC staging systems, the JES classification was more effective in detecting NSM and revealing EMT activation in NSM compared with non-NSM. Given the established link between EMT and poor prognosis, this observation suggests that patients with NSM may have worse outcomes, supporting a potential prognostic role for NSM. In comparison, the AJCC/UICC system, which defines the N stage solely by the number of metastatic lymph nodes, classified all nine patients as N2, obscuring critical genetic distinctions. These findings underscore the need to integrate genetic validation with JES classification to more accurately identify NSM and assess its prognostic relevance [32]. Moreover, although sequencing of both primary tumors and lymph nodes is informative, it is often challenging and impractical in routine clinical practice. Based on the unsupervised hierarchical clustering of EMT markers, we propose that their expression levels in primary tumors may serve as a surrogate indicator of NSM when multiregional genome sequencing is not feasible.

This study has several limitations. First, the rarity of patients with multiple metastatic lymph nodes spanning different JES-defined groups limited cohort recruitment and led to a relatively small sample size. Second, some pathology-positive lymph nodes could not be sequenced because of limited tissue or poor DNA quality, which may have reduced the completeness of phylogenetic reconstruction and the resolution of inferred dissemination patterns; a ‘phantom skip’ cannot be excluded. We therefore regard NSM as a genetically inferred pattern constrained by node availability and advise caution in patient-level interpretation. Nonetheless, unsequenced nodes were typically small deposits with limited tumor content, making it unlikely that major clonal lineages were missed. Third, the choice of nodal classification system may influence interpretations, thus findings should be considered within the framework applied.

While these factors may limit generalizability, the rarity of such cases highlights the value of this cohort in providing novel insights into LNM trajectories and the potential clinical relevance of NSM. Further studies with larger cohorts would be needed to validate our observations, enhance prognostic accuracy, and optimize treatment strategies for ESCC. If confirmed, this approach could offer a cost-effective alternative to multiregional genome sequencing, enabling earlier and more accessible NSM detection in clinical practice.

In summary, this study provides an in-depth understanding of the genomic and transcriptomic landscape of ESCC, with particular focus on LNM. By integrating advanced genomic and phylogenetic analyses, we identified key features such as node-skipping patterns and EMT activation in NSM, proposing their potential as prognostic significance. Multiregional sequencing further uncovered genetic and molecular relationships among metastatic nodes, offering novel perspectives on metastatic behavior. Together, we anticipate refining diagnostic precision and improving patient outcomes in ESCC.

Electronic Supplementary Material

Supplementary materials are available at Cancer Research and Treatment website (https://www.e-crt.org).

crt-2025-257_S1_Table.pdf

crt-2025-257_S2_Table.pdf

crt-2025-257_S3_Fig.pdf

crt-2025-257_S4_Fig.pdf

crt-2025-257_S5_Table.pdf

NOTES

Ethical Statement

All procedures were conducted in accordance with the Declaration of Helsinki, and informed consent was obtained from each participant. This study was approved by the Institutional Review Board of Yonsei University Health System (IRB approval no. 4-2018-1210).

Author Contributions

Conceived and designed the analysis: Kim DJ, Kim S.

Collected the data: Park SY, Kim HE.

Contributed data or analysis tools: Park J, Jo SY, Won J.

Performed the analysis: Park J.

Wrote the paper: Park J, Park SY.

Project administration: Park SY, Kim DJ, Kim S.

Conflicts of Interest

Sangwoo Kim is cofounder of AIMA Inc., which seeks to develop techniques for early cancer diagnosis based on circulating tumor DNA. The other authors declare no competing interests.

Funding

This research was supported by the Bio&Medical Technology Development Program of the National Research Foundation funded by the Ministry of Science and ICT, Republic of Korea [RS-2023-00261820] and by a National Research Foundation of Korea grant funded by the Korean government (Ministry of Science and ICT) [2022R1A2C209310611].

Data Availability

Raw sequence reads have been deposited at SRA as PRJNA1224845. They are available upon request if access is granted. To request access, please contact Sangwoo Kim. All original code has been deposited at ESCC-LNM and is publicly available at https://github.com/jennyp76/ESCC-LNM as of the date of publication.

Fig. 1.

Mutational spectrum of multiregional whole-exome sequencing (WES) in 10 esophageal squamous cell carcinoma patients. (A) Locations of the primary tumors and lymph node metastasis sites harvested for WES. Names of the lymph nodes are presented according to the Japanese Classification of Esophageal Cancer. (B) Tumor mutation burden (TMB) (top) and distribution of somatic mutations (n=35 samples, bottom).

Fig. 2.

Temporal order of lymph node metastasis in 10 esophageal squamous cell carcinoma patients generated using phylogenetic analysis. (A) Schematic representation of the clonal evolution used to infer the temporal order of tumor development. The primary tumor acquires somatic mutations, leading to selective sweeps and the emergence of new subclones. Lymph node metastasis can originate from different subclones at multiple time points (t₁, t₂, and t₃) and may undergo additional evolution within the lymph nodes. Phylogenetic trees were constructed for each patient by assessing the genetic similarity between samples, using a paired normal sample as the outgroup root. (B) Phylogenetic trees were constructed using the Wagner parsimony method with the PHYLIP software package. Colors indicate the mutation’s distribution: dark blue if detected in all regions, orange if in not all but some regions, and beige if in a single region. The branch lengths are proportional to the number of acquired mutations. Known cancer driver genes are annotated according to their predicted time of acquisition. The anatomic chart provides the spatial locations of the primary tumor and lymph nodes. (C) Distribution of prelymph node metastatic (pre-LN-metastatic) somatic mutations. Black circles indicate cancer driver genes. Post-LN-metastatic, post–lymph node metastatic.

Fig. 3.

Mutational characteristics in lymph node metastasis trajectory. (A) Overall spectra of base-substitution contributions in total mutations (top), pre–lymph node metastatic (pre-LN-metastatic) mutations (middle), and post–lymph node metastatic (post-LN-metastatic) mutations (bottom) across all patients. (B) Mutation signature analysis at the patient level using the COSMIC Mutational Signature Framework version 3. A dot in the boxplot represents a single patient. SBS16, unknown; SBS5, unknown (clock-like); SBS6, defective DNA mismatch repair; SBS1, spontaneous deamination of 5-methylcytosine (clock-like); SBS18, damage by reactive oxygen species. p-values were calculated using the Mann-Whitney U test. Boxplots are represented by a centerline, median; box, 25th to 75th percentile (interquartile range [IQR]); and whiskers, data within 1.5 times the IQR. ^*p < 0.05, ^**p < 0.01. (C) Positive correlation between the timing of metastasis (ToM) and mutation accumulation rate after metastasis (MARM). ToM and MARM were calculated per patient per lymph node, using the lengths of the tree (S4A Fig.). Spearman’s correlation and p-values are shown. The positive tendency within patients is shown in S4C Fig.

Fig. 4.

Identification of nodal skip metastasis (NSM) and its characteristics in esophageal squamous cell carcinoma tumors. (A) Identification of non-NSM and NSM. A total of three patients (3/9, 33.3%) were classified as NSM, including two patients (P6 and P8) who were additionally included through phylogenetic analysis. Lymph nodes are arranged by the order of metastasis, and the darkness of the circle represents the grade of a lymph node classified by the Japanese Classification of Esophageal Cancer. (B) Gene set enrichment analysis of NSM vs. non-NSM. Positive values of normalized enrichment score (NES) indicate enrichment in NSM, and negative values imply enrichment in non-NSM. (C) Unsupervised hierarchical clustering using gene expression profiles of 10 epithelial-mesenchymal transition (EMT) markers. These EMT markers are commonly used in immunohistochemistry staining to assess differences in EMT activation levels. Patients were clustered as NSM, non-NSM, and normal, corresponding to the metastatic patterns observed in the phylogenetic tree analysis, and markers were classified as mesenchymal and epithelial, consistent with the acknowledged information.

Table 1.

Clinical characteristics of the ESCC cohort

Patient	Sex/Age (yr)	Location	pTN stage (AJCC)	No. dissected lymph nodes	Lymph node metastasis
P1	M/71	Upper	T3N2	91	104R (N2), 106recR (N1)
P2	F/67	Middle	T1bN1	76	104R (N2)
P3	M/73	Middle	T3N2	56	107 (N2), 106tbL (N3), 101L (N2)^a)
P4	M/63	Middle	T3N2	85	7 (N2), 106recL (N1), 112aoA (N2)
P5	F/51	Upper	T3N2	68	104R (N2), 106recR (N1), 106recL (N1)
P6	M/65	Middle	T3N2	92	1 (N1), 107 (N2), 106recL (N1)
P7	M/69	Middle	T1bN2	26	1 (N1), 110 (N2), 106recL (N1)
P8	M/72	Upper	T3N2	75	106recR (N1), 108 (N2), 1 (N3)^a)
P9	M/62	Middle	T2N2	62	7 (N2), 105 (N2), 106recR (N1)
P10	M/76	Middle	T4N2	68	1 (N1), 3a (N1), 106tbL (N3)

7, lymph nodes along the left gastric artery; 106recL, left recurrent nerve lymph nodes; 108, middle thoracic paraesophageal lymph nodes; 1, right paracardial lymph node; 3a, lesser curvature lymph nodes along the branches of the left gastric artery; 106tbL, left tracheobronchial lymph nodes; 104R, right supraclavicular lymph nodes; 106recR, right recurrent nerve lymph nodes; 107, subcarinal lymph nodes; 105, upper thoracic paraesophageal lymph nodes; 110, lower thoracic paraesophageal lymph nodes; 112aoA, anterior thoracic paraaortic lymph nodes; 101L, left cervical paraesophageal lymph nodes; AJCC, American Joint Committee on Cancer; ESCC, esophageal squamous cell carcinoma.

^a) Lymph nodes filtered out for low quality (see Materials and Methods).

REFERENCES

1. Park SY, Kim DJ. Esophageal cancer in Korea: epidemiology and treatment patterns. J Chest Surg. 2021;54:454–9. Article PubMed PMC
2. Pennathur A, Gibson MK, Jobe BA, Luketich JD. Oesophageal carcinoma. Lancet. 2013;381:400–12. Article PubMed
3. Chung HL, Sun J, Leung JWT. Breast cancer skip metastases: frequency, associated tumor characteristics, and role of staging nodal ultrasound in detection. AJR Am J Roentgenol. 2021;217:835–44. Article PubMed
4. Warshavsky A, Rosen R, Nard-Carmel N, Abu-Ghanem S, Oestreicher-Kedem Y, Abergel A, et al. Assessment of the rate of skip metastasis to neck level IV in patients with clinically node-negative neck oral cavity squamous cell carcinoma: a systematic review and meta-analysis. JAMA Otolaryngol Head Neck Surg. 2019;145:542–8. Article PubMed PMC
5. Wang F, Zheng Y, Wang Z, Zheng Q, Huang Q, Liu S. Nodal skip metastasis in esophageal squamous cell carcinoma patients undergoing three-field lymphadenectomy. Ann Thorac Surg. 2017;104:1187–93. Article PubMed
6. Zhu Z, Yu W, Li H, Zhao K, Zhao W, Zhang Y, et al. Nodal skip metastasis is not a predictor of survival in thoracic esophageal squamous cell carcinoma. Ann Surg Oncol. 2013;20:3052–8. Article PubMed PDF
7. Cavallin F, Alfieri R, Scarpa M, Cagol M, Ruol A, Fassan M, et al. Nodal skip metastasis in thoracic esophageal squamous cell carcinoma: a cohort study. BMC Surg. 2017;17:49.Article PubMed PMC PDF
8. Um SW, Joung JG, Lee H, Kim H, Kim KT, Park J, et al. Molecular evolution patterns in metastatic lymph nodes reflect the differential treatment response of advanced primary lung cancer. Cancer Res. 2016;76:6568–76. Article PubMed PDF
9. Rice TW, Blackstone EH, Rusch VW. 7th edition of the AJCC Cancer Staging Manual: esophagus and esophagogastric junction. Ann Surg Oncol. 2010;17:1721–4. Article PubMed PDF
10. Japan Esophageal Society. Japanese classification of esophageal cancer, 11th edition: part I. Esophagus. 2017;14:1–36. Article PubMed PDF
11. Stachler MD, Taylor-Weiner A, Peng S, McKenna A, Agoston AT, Odze RD, et al. Paired exome analysis of Barrett’s esophagus and adenocarcinoma. Nat Genet. 2015;47:1047–55. Article PubMed PMC PDF
12. Bailey MH, Tokheim C, Porta-Pardo E, Sengupta S, Bertrand D, Weerasinghe A, et al. Comprehensive characterization of cancer driver genes and mutations. Cell. 2018;173:371–85. PubMed PMC
13. Ng PC, Henikoff S. SIFT: predicting amino acid changes that affect protein function. Nucleic Acids Res. 2003;31:3812–4. Article PubMed PMC
14. Choi Y, Sims GE, Murphy S, Miller JR, Chan AP. Predicting the functional effect of amino acid substitutions and indels. PLoS One. 2012;7:e46688Article PubMed PMC
15. Ioannidis NM, Rothstein JH, Pejaver V, Middha S, McDonnell SK, Baheti S, et al. REVEL: an Ensemble method for predicting the pathogenicity of rare missense variants. Am J Hum Genet. 2016;99:877–85. Article PubMed PMC
16. Oshima K, Kato K, Ito Y, Daiko H, Nozaki I, Nakagawa S, et al. Prognostic biomarker study in patients with clinical stage I esophageal squamous cell carcinoma: JCOG0502-A1. Cancer Sci. 2022;113:1018–27. Article PubMed PMC PDF
17. Wang L, Jia YM, Zuo J, Wang YD, Fan ZS, Feng L, et al. Gene mutations of esophageal squamous cell carcinoma based on next-generation sequencing. Chin Med J (Engl). 2021;134:708–15. Article PubMed PMC
18. Qing T, Zhu S, Suo C, Zhang L, Zheng Y, Shi L. Somatic mutations in ZFHX4 gene are associated with poor overall survival of Chinese esophageal squamous cell carcinoma patients. Sci Rep. 2017;7:4951.Article PubMed PMC PDF
19. Dhar SS, Lee MG. Cancer-epigenetic function of the histone methyltransferase KMT2D and therapeutic opportunities for the treatment of KMT2D-deficient tumors. Oncotarget. 2021;12:1296–308. Article PubMed PMC
20. Qu Y, Gao N, Wu T. Expression and clinical significance of SYNE1 and MAGI2 gene promoter methylation in gastric cancer. Medicine (Baltimore). 2021;100:e23788Article PubMed PMC
21. Li P, Xiao J, Zhou B, Wei J, Luo J, Chen W. SYNE1 mutation may enhance the response to immune checkpoint blockade therapy in clear cell renal cell carcinoma patients. Aging (Albany NY). 2020;12:19316–24. Article PubMed PMC
22. Wang K, Li Z, Xuan Y, Zhao Y, Deng C, Wang M, et al. Pan-cancer analysis of NFE2L2 mutations identifies a subset of lung cancers with distinct genomic and improved immunotherapy outcomes. Cancer Cell Int. 2023;23:229.Article PubMed PMC PDF
23. Moody S, Senkin S, Islam SMA, Wang J, Nasrollahzadeh D, Cortez Cardoso Penha R, et al. Mutational signatures in esophageal squamous cell carcinoma from eight countries with varying incidence. Nat Genet. 2021;53:1553–63. Article PubMed PDF
24. Chen M, Linstra R, van Vugt M. Genomic instability, inflammatory signaling and response to cancer immunotherapy. Biochim Biophys Acta Rev Cancer. 2022;1877:188661.Article PubMed
25. Takahashi H, Oshi M, Yan L, Endo I, Takabe K. Gastric cancer with enhanced apical junction pathway has increased metastatic potential and worse clinical outcomes. Am J Cancer Res. 2022;12:2146–59. PubMed PMC
26. Taube JH, Herschkowitz JI, Komurov K, Zhou AY, Gupta S, Yang J, et al. Core epithelial-to-mesenchymal transition interactome gene-expression signature is associated with claudin-low and metaplastic breast cancer subtypes. Proc Natl Acad Sci U S A. 2010;107:15449–54. Article PubMed PMC
27. Kim RK, Suh Y, Yoo KC, Cui YH, Kim H, Kim MJ, et al. Activation of KRAS promotes the mesenchymal features of basal-type breast cancer. Exp Mol Med. 2015;47:e137Article PubMed PMC PDF
28. Tajiri A, Ishihara R, Sakurai H, Nakamura T, Tani Y, Inoue T, et al. Clinical features of superficial esophagus squamous cell carcinoma according to alcohol-degrading enzyme ADH1B and ALDH2 genotypes. J Gastroenterol. 2022;57:630–9. Article PubMed PDF
29. Li XC, Wang MY, Yang M, Dai HJ, Zhang BF, Wang W, et al. A mutational signature associated with alcohol consumption and prognostically significantly mutated driver genes in esophageal squamous cell carcinoma. Ann Oncol. 2018;29:938–44. Article PubMed PMC
30. Byeon J, Lim C, Kang E, Jung JJ, Kim HK, Lee HB, et al. Comparison of long-term oncological outcome of sentinel lymph node mapping methods (dye-only versus dye and radioisotope) in breast cancer patients following neoadjuvant chemotherapy. Cancer Res Treat. 2026;58:175–81. Article PubMed PDF
31. Cascinelli N, Belli F, Santinami M, Fait V, Testori A, Ruka W, et al. Sentinel lymph node biopsy in cutaneous melanoma: the WHO Melanoma Program experience. Ann Surg Oncol. 2000;7:469–74. Article PubMed PDF
32. Shang QX, Yang YS, Xu LY, Yang H, Li Y, Li Y, et al. Prognostic role of nodal skip metastasis in thoracic esophageal squamous cell carcinoma: a large-scale multicenter study. Ann Surg Oncol. 2021;28:6341–52. Article PubMed PDF

Figure & Data

REFERENCES

Citations

Citations to this article as recorded by

PubReader
ePub Link

Cite

CITE: export Copy Download Format; Close

Download Citation

Download a citation file in RIS format that can be imported by all major citation management software, including EndNote, ProCite, RefWorks, and Reference Manager.

Format:

RIS — For EndNote, ProCite, RefWorks, and most other reference management software
BibTeX — For JabRef, BibDesk, and other BibTeX-specific software

Include:

Citation for the content below
Citation and abstract for the content below

Trajectory of Lymph Node Metastasis Reflects the Molecular Features of Esophageal Squamous Cell Carcinoma

DOI: https://doi.org/10.4143/crt.2025.257

XML Download

Abstract
Introduction
Materials and Methods
Results
Discussion
Electronic Supplementary Material
NOTES
REFERENCES

Trajectory of Lymph Node Metastasis Reflects the Molecular Features of Esophageal Squamous Cell Carcinoma

Fig. 1. Mutational spectrum of multiregional whole-exome sequencing (WES) in 10 esophageal squamous cell carcinoma patients. (A) Locations of the primary tumors and lymph node metastasis sites harvested for WES. Names of the lymph nodes are presented according to the Japanese Classification of Esophageal Cancer. (B) Tumor mutation burden (TMB) (top) and distribution of somatic mutations (n=35 samples, bottom).

Fig. 2. Temporal order of lymph node metastasis in 10 esophageal squamous cell carcinoma patients generated using phylogenetic analysis. (A) Schematic representation of the clonal evolution used to infer the temporal order of tumor development. The primary tumor acquires somatic mutations, leading to selective sweeps and the emergence of new subclones. Lymph node metastasis can originate from different subclones at multiple time points (t1, t2, and t3) and may undergo additional evolution within the lymph nodes. Phylogenetic trees were constructed for each patient by assessing the genetic similarity between samples, using a paired normal sample as the outgroup root. (B) Phylogenetic trees were constructed using the Wagner parsimony method with the PHYLIP software package. Colors indicate the mutation’s distribution: dark blue if detected in all regions, orange if in not all but some regions, and beige if in a single region. The branch lengths are proportional to the number of acquired mutations. Known cancer driver genes are annotated according to their predicted time of acquisition. The anatomic chart provides the spatial locations of the primary tumor and lymph nodes. (C) Distribution of prelymph node metastatic (pre-LN-metastatic) somatic mutations. Black circles indicate cancer driver genes. Post-LN-metastatic, post–lymph node metastatic.

Fig. 3. Mutational characteristics in lymph node metastasis trajectory. (A) Overall spectra of base-substitution contributions in total mutations (top), pre–lymph node metastatic (pre-LN-metastatic) mutations (middle), and post–lymph node metastatic (post-LN-metastatic) mutations (bottom) across all patients. (B) Mutation signature analysis at the patient level using the COSMIC Mutational Signature Framework version 3. A dot in the boxplot represents a single patient. SBS16, unknown; SBS5, unknown (clock-like); SBS6, defective DNA mismatch repair; SBS1, spontaneous deamination of 5-methylcytosine (clock-like); SBS18, damage by reactive oxygen species. p-values were calculated using the Mann-Whitney U test. Boxplots are represented by a centerline, median; box, 25th to 75th percentile (interquartile range [IQR]); and whiskers, data within 1.5 times the IQR. *p < 0.05, **p < 0.01. (C) Positive correlation between the timing of metastasis (ToM) and mutation accumulation rate after metastasis (MARM). ToM and MARM were calculated per patient per lymph node, using the lengths of the tree (S4A Fig.). Spearman’s correlation and p-values are shown. The positive tendency within patients is shown in S4C Fig.

Fig. 4. Identification of nodal skip metastasis (NSM) and its characteristics in esophageal squamous cell carcinoma tumors. (A) Identification of non-NSM and NSM. A total of three patients (3/9, 33.3%) were classified as NSM, including two patients (P6 and P8) who were additionally included through phylogenetic analysis. Lymph nodes are arranged by the order of metastasis, and the darkness of the circle represents the grade of a lymph node classified by the Japanese Classification of Esophageal Cancer. (B) Gene set enrichment analysis of NSM vs. non-NSM. Positive values of normalized enrichment score (NES) indicate enrichment in NSM, and negative values imply enrichment in non-NSM. (C) Unsupervised hierarchical clustering using gene expression profiles of 10 epithelial-mesenchymal transition (EMT) markers. These EMT markers are commonly used in immunohistochemistry staining to assess differences in EMT activation levels. Patients were clustered as NSM, non-NSM, and normal, corresponding to the metastatic patterns observed in the phylogenetic tree analysis, and markers were classified as mesenchymal and epithelial, consistent with the acknowledged information.

Fig. 1.

Fig. 2.

Fig. 3.

Fig. 4.

Trajectory of Lymph Node Metastasis Reflects the Molecular Features of Esophageal Squamous Cell Carcinoma

Patient	Sex/Age (yr)	Location	pTN stage (AJCC)	No. dissected lymph nodes	Lymph node metastasis
P1	M/71	Upper	T3N2	91	104R (N2), 106recR (N1)
P2	F/67	Middle	T1bN1	76	104R (N2)
P3	M/73	Middle	T3N2	56	107 (N2), 106tbL (N3), 101L (N2)^a)
P4	M/63	Middle	T3N2	85	7 (N2), 106recL (N1), 112aoA (N2)
P5	F/51	Upper	T3N2	68	104R (N2), 106recR (N1), 106recL (N1)
P6	M/65	Middle	T3N2	92	1 (N1), 107 (N2), 106recL (N1)
P7	M/69	Middle	T1bN2	26	1 (N1), 110 (N2), 106recL (N1)
P8	M/72	Upper	T3N2	75	106recR (N1), 108 (N2), 1 (N3)^a)
P9	M/62	Middle	T2N2	62	7 (N2), 105 (N2), 106recR (N1)
P10	M/76	Middle	T4N2	68	1 (N1), 3a (N1), 106tbL (N3)

Table 1. Clinical characteristics of the ESCC cohort

Lymph nodes filtered out for low quality (see Materials and Methods).