A total of 13 ESCC patients from Zhejiang Province, China, after esophagectomy were enrolled, in which 8 cases had two primary lesions and 5 cases had a single lesion. Fresh tumor tissues and PNM were obtained during the surgical operation. Totally 21 tumor and 13 PNM samples of ESCC were collected. We resected the PNM samples from the upper or lower part of the esophageal mucosa when the primary lesions were present in the lower or upper esophagus. Furthermore, the distance between each PNM and its closest tumor was more than 5 cm. Besides, we collected esophageal PNM samples from five patients with gastric cancer (GC) undergoing total gastrectomy. All their PNM samples were collected in the upper esophagus from > 10 cm far away from the GC lesions, and no lesion in the esophagus was found in these five GC patients. All cases were operated on in the Cancer Hospital of the University of Chinese Academy of Sciences (Zhejiang Cancer Hospital) from April 2019 to February 2021. Diagnosis, validity, and tumor purity of the specimens were confirmed by two independent pathologists of Zhejiang Cancer Hospital. In addition, 5 ml peripheral blood was collected from each patient and placed into EDTA-coated tubes (BD Biosciences). White blood cells (WBCs) were extracted as the control to determine germline variations, and our study was focused on somatic alterations. This study was approved by the ethics committee of Zhejiang Cancer Hospital (Approval No. IRB-2022-154). All participants were informed and consent to sample collection, intended research, and publication usage. Written consent was collected according to the ethical regulations of Zhejiang Cancer Hospital. The next-generation sequencing (NGS) was performed in a Clinical Laboratory Improvement Amendments (CLIA)- certified and College of American Pathologists (CAP)- accredited clinical testing laboratory (Nanjing Geneseeq Technology Inc., China). All samples were shipped to the clinical testing laboratory following the required conditions.

Genomic DNA from fresh tumor tissue, PNM, and WBCs was extracted using DNeasy Blood & Tissue kit (Qiagen) following the manufacturer’s instruction. Purified genomic DNA was qualified by Nanodrop2000 for A260/280 and A260/A230 ratios (Thermo Fisher Scientific). All DNA samples were quantified by Qubit 3.0 using the dsDNA HS Assay Kit (Life Technologies) according to the manufacturer’s recommendations.

For WES library construction, we fragmented 2μg DNA using Covaris M220 sonication system (Covaris), followed by end-repairing, A-tailing, and adaptor ligation and purification by KAPA Hyper Prep Kit (KAPA Biosystems). The resultant libraries were amplified and purified before exome capture using the xGen Exome Research Panel v1.0 (Integrated DNA Technologies). The enriched libraries were then sequenced using the Illumina HiSeq 4000 platform with 2×150 bp pair-end reads. The mean raw coverage depth was ~60× for the WBC samples and ~200× for the tumor and PNM samples.

Paired-end sequencing data were aligned to the reference human genome (build hg19) with the Burrows-Wheeler Aligner (bwa-mem) (7). Alignment results (BAM files) were further processed for de-duplication, base quality recalibration, and indel realignment using the Picard suite (http://picard.sourceforge.net/) and the Genome Analysis Toolkit (GATK) (8). MuTect with default parameters was applied to the paired PNM and tumor BAM files to identify somatic single nucleotide variants (SNVs) (9). SNVs in the 1000 Genomes project and dbSNP with frequency >1% were excluded. Small insertions and deletions (indels) were detected using SCALPEL (10). SNV and indel annotation was performed by ANNOVAR using the hg19 reference genome and 2014 versions of standard databases and functional prediction programs (11). Gene-level copy number ratios were calculated by CNVKit using the CNVKit algorithm, relative copy-ratios for each exon were calculated by correcting for imbalanced library size, GC bias, sequence repeats, and target density. The log2 ratio values of 2.0 and 0.6 were used as the cut-off for copy number gain and copy number loss of tissue samples, respectively. Chromosome arm-level somatic copy number variations (CNVs) were analyzed by FACETS with a 0.2 drift cut-off for unstable joint segments. Chromosome instability score (CIS) was defined as the proportion of the genome with aberrant (purity-adjusted segment-level copy number >=3 or <=1) segmented copy number. Mutational signature enrichment weights were calculated using the sigminer R package (12). Treeomics (13) was used to reconstruct the phylogenetic relationships with maximum likelihood. Different tissues that were grouped together into the same clades were determined as convergent, otherwise, they were considered as divergent. The mutant-allele tumor heterogeneity (MATH) analysis was performed using inferHeterogeneity in the Maftools package (14).

Quantitative data were displayed as the median value (range) or the number of patients (percentage). Comparisons of proportion between two groups were done using Fisher’s exact test. Wilcoxon rank-sum test was performed to compare the mutation number, CIS, and MATH between different groups. Differences in mutation number among PNM of ESCC, PNM of GC, and tumor of ESCC were analyzed using the Kruskal-Wallis test. A two-sided P value of less than 0.05 was considered significant for all tests unless indicated otherwise. All statistical analyses were done in R (v.3.6.0).

As demonstrated in Figure S1 , we enrolled 13 ESCC patients (double-primary: 8; single-primary: 5) whose tumor and PNM samples were collected to perform WES analysis. Five GC patients were also enrolled whose PNM samples were collected and subject to WES. Participants’ demographics and clinical characteristics are listed in Tables 1 and S1 . The median age of the 13 ESCC patients was 66 years. 12 of them were male. 7, 2, 10, and 11 of the patients had hypertension history, ESCC family history, smoking, and drinking history, respectively. The study included Stage I, II, and III patients (1, 4, and 8 cases, respectively). For tumor differentiation level, 3 of the ESCC patients were defined as high differentiation, 4 were moderate, and the other 6 were low ( Table 1 ). All the GC participants were male, with a median age of 71 years ( Table S1 ).

To understand the molecular characteristics underlying ESCC, we set out mutational profiling using tumor and normal mucosa tissues from ESCC patients. The WES results in Figure 1 outlined the key differences in genomic alterations between tumor and normal mucosa. We did not observe significant enrichment of any somatically mutated gene in the tumor sample. Conversely, the mutation frequencies of NOTCH1 and NOTCH2 are significantly higher in PNM ( Figure 1 , top panel). The frequency of copy number changes was also examined, and we found a significant increase associated with CCND1 in tumors ( Figure 1 , middle panel). In all the chromosomal arms, the frequency of copy number variation (CNV) in tumors was remarkably higher than that in PNM ( Figure 1 , bottom panel).

We further investigated the distribution of somatic gene mutations in ESCC tumors and PNM. In total, we found 854 mutations in the esophageal PNM and 1571 mutations in the tumors ( Figure 2A ). The median number of mutations in tumors (71, range 23-166) is higher than that in PNM (56, range 31-128) but lacks statistical significance (p=0.290, Wilcoxon test) ( Figure 3A ). Ten of the somatic mutations are found to present in both tumor and PNM, while 8 of them are actually from Patient P6 ( Figures 2A, B ). Detailed information about these tumor-PNM shared mutations was summarized in Table S2 . In addition, we scrutinized the somatic mutations for their base mutation patterns, mutation types, and functional outcomes (inactivation vs. non-inactivation) and found high-degree similarity between tumor and PNM from ESCC ( Figure S2 ). The somatic mutation pattern was also analyzed for the NOTCH genes, which are more frequently mutated in PNM ( Figures 1 and S3 ). There were many code-shifting mutations of NOTCH1 in PNM from the mutation type, but no significant difference. The mutation frequency of NOTCH2 gene in PNM was significantly higher than that in the tumor. However, there was no significant difference between the mutation types.

The genome stability conditions varied between the tumor and PNM samples. We calculated the chromosome instability score (CIS) of tumor and PNM, and the scores in the tumors are significantly higher than those in PNM ( Figure 3B ). The CNV analysis revealed significantly higher incidences of copy number amplification in the tumor than in PNM at both gene and chromosomal arm levels (P= 0.003 and 0.009, respectively, Fisher’s exact test, Figures 3C, D ). As for the incidence of copy number deletion, tumor and PNM samples are similar at the gene level ( Figure 3C ). The incidence of chromosomal arm level copy number deletion in the tumor is moderately higher than that in PNM (P=0.139, Fisher’s exact test, Figure 3D ).

For signaling pathway analysis, we classified key cancer-associated genes into ten canonical mitogenic signaling pathways (15). Then, we counted the number of samples that have at least one gene mutated in each pathway and computed the proportion of total samples altered in each pathway. Compared with normal mucosa, ESCC tumors exhibited significantly higher mutation frequencies in several cancer-related pathways, including the HIPPO (P=0.038, Fisher’s exact test), WNT (P=0.013, Fisher’s exact test), and NRF2 (P=0.029, Fisher’s exact test) pathways ( Figure 4A ). The NOTCH pathway is highly mutated in PNM relative to the tumor but lacks statistical significance ( Figure 4A ). Additionally, no somatic mutation was detected in genes of the NRF2 or WNT pathway in PNM ( Figure 4B ), but these signaling pathways are frequently mutated in the tumor. We also performed mutation signature analysis ( Figure S4 ) and found that the Age signature (SBS1) is highly prevalent in both PNM and tumor (16). The enrichment of APOBEC signature (SBS2) in the tumor is significantly higher than that in PNM (P=0.021, Wilcoxon test). The signatures of Ultraviolet (SBS7, P=0.081, Wilcoxon test), POLE (SBS10a/b, P=0.092, Wilcoxon test), MMRdeficiency (SBS15, P=0.550, Wilcoxon test), and BRCA (SBS3, P=0.630, Wilcoxon test) in the tumor are also more enriched but lack statistical significance.

We performed phylogenetic analyses of somatic SNVs on specimens from 8 ESCC patients with multiple primary tumors to examine their interfocal heterogeneity, as well as the relationships between tumor and normal mucosa. Based on the genomic similarity among the two primary tumors and PNM of the same case, the genetic divergence patterns of the 8 patients could be categorized into three modes. In Mode I, the PNM and two tumor samples are mutually divergent and distribute at different clades. In Mode II, one of the two tumor samples shows convergence with PNM, and they cluster in the same clade. In Mode III, the two tumor samples but no PNM have convergence and cluster in the same clade. We found 5 patients in Mode I, 1 patient in Mode II, and 2 patients in Mode III ( Figure 5A ).

For every tumor and PNM sample, the MATH score was calculated, grouped into the aforementioned three modes, and compared to evaluate the heterogeneity of the different modes. We observed a gradual downward trend of the MATH scores from Mode I to Modes II and III ( Figure 5B ), consistent with the patterns shown in the hierarchical clustering dendrograms ( Figure 5A ).

To explore the unique genomic alterations in ESCC, we incorporated genomic analysis on esophageal PNM of GC patients and made the comparison with ESCC tissue samples. Remarkably, there is a significant difference in the prevalence of mutations in the TP53 gene ( Figures 1 and S5 ), as this gene is commonly mutated in both the tumor and PNM samples of ESCC but not mutated in the GC samples (p=0.0003, P=0.002). The detailed TP53 somatic mutation information is listed in Table S3 . We observed that the most frequent TP53 single base substitution (SBS) is the G:C>T:A transversion (26.2%) related to smoking in ESCC, followed by the G:C>A:T transition (19.0%) probably associated with chronic inflammation and the aristolochic acid-related A:T>T:A transversion (19.0%) (17). The genomic profiling of GC PNM also identified alterations in other cancer-related genes, including NOTCH1, NOTCH2, NF1, and STK11 ( Figure S5A ). The median number of mutations per sample was 15 (5–62) in PNM of GC, 56 (31–128) in PNM of ESCC patients, and 71 (23–166) in the tumor of ESCC ( Figure S5B ). The number of mutations in PNM of GC was significantly lower than that in ESCC tumor (P=0.021, Wilcoxon test) and was also lower than that in ESCC PNM but not significant (P=0.100, Wilcoxon test). Based on the signaling pathway analysis comparing the ESCC and GC samples, we noticed that tumor and PNM of ESCC patients had higher frequencies of mutations in cancer-related signaling pathways in general ( Figure S5C ). Specifically, mutations are significantly enriched in the tumor and PNM of ESCC patients than that of GC PNM (P=0.002 and 0.022, respectively, Fisher’s exact test) for the TP53 signaling pathway. The comparison of mutation signatures is shown in Figure S6 , showing no significant difference between ESCC and GC patients in the prevalent mutation signatures.