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

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.

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.

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

Table 1.

Clinical characteristics of the ESCC cohort

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

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.

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

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.

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.

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.

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