Four GEJ cancer patients were enrolled from Tianjin Medical University Cancer Institute and Hospital (from May 2020 to October 2020). Inclusion criteria: (1) treatment naive patients; (2) patients undergoing radical surgery. The information of the GEJ cancer patients was summarized in Table S1 . All the four patients’ surgical primary tumor tissues and adjacent normal tissues were collected. Three patients were diagnosed with positive lymph nodes. Positive lymph nodes of two patients were also collected (The lymph nodes in one patient were not suited for single cell sequencing). This study was approved by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (NO.bc2020180).

The acquired GEJ tumor tissues and positive lymph nodes were mechanically dissociated, part of which was used to extract genomic DNA and part was used to isolate single cells. The normal tissues were only used to extract genomic DNA. The genomic DNAs were extracted with QIAamp^® DNA Micro Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s protocol.

For single cell isolation, the tissues were digested with tumor dissociation kit (Miltenyi Biotec, Bergisch Gladbach, Germany). Erythrocytes were depleted by red blood cell lysis buffer (Solarbio, Beijing, China). Leukocytes were removed with Dynabeads™ CD45 (Invitrogen, Waltham, MA, USA). The single cells were isolated with mouth pipetting and transferred to tubes containing the lysis buffer (0.02 M Tris-HCl, 2 mM EDTA, 0.02 M KCl, 0.015 M DTT, 0.25 μM 5N3G Primer, 5 μg QIAGEN protease). The cells were lysed and stored at −80°C until amplification.

Single cell amplification was carried out according to the previously reported MALBAC method (24). This whole genome amplification method introduced quasi-linear preamplification, which could reduce the amplification bias. In brief, nine cycle preamplification and subsequent exponential amplification were performed. The amplification products were purified with AMPure XP beads (Beckman Coulter, Brea, CA, USA) and quantified with Qubit 4.0 Fluorometer (Invitrogen, Waltham, MA, USA). Quantitative PCR (qPCR) was used to check the integrity of the genome (the sequences of qPCR primers were shown in Table S2 ). Cells having no less than six out of eight randomly selected loci with Ct value<30 could be used for library preparation and next generation sequencing.

For single cells, we constructed the whole genome library for each cell. The selected single cell DNAs were sonicated with Covaris S220 to produce short DNA fragments. Then, libraries were prepared with KAPA Hyper Prep Kit (KAPA, Boston, MA, USA) and NEBNext Multiplex Oligos for Illumina (NEB, Ipswich, MA, USA) according to manufacturer’s protocol. Briefly, main steps included end repair and A-tailing, adapter ligation, post-ligation cleanup and library amplification. qPCR and fragment analysis were used for library quality control. The selected libraries were sequenced with Illumina NovaSeq 6000.

The exome libraries were constructed with Agilent SureSelect Human All Exon V6 kit (Agilent Technologies, Santa Clara, CA, USA) according to manufacturer’s protocol. Briefly, the genomic DNA was sonicated with Covaris S220 to produce short DNA fragments. After end repair with exonuclease/polymerase, A-tailing in 3’ ends of DNA fragments and adapter ligation, the DNA fragments with adapter at both ends were enriched by PCR. The biotin labeled probe was used to hybridize with DNA libraries and the exons were captured with magnetic beads and streptomycin. Then, the captured libraries were added index tags and enriched by PCR reaction. qPCR and fragment analysis were used for library quality control. The selected libraries were sequenced with Illumina NovaSeq 6000. The sequencing data has been submitted to NCBI SRA under BioProject accession number PRJNA718709.

First, the raw sequencing reads were trimmed to remove adaptors and low quality bases using Trimmomatic (25). Pair-End mode with ‘HEADCROP:35’ was used to remove the low quality bases from the start of the read and other default parameters. Then, the clean reads were mapped to the human reference genome hg38 with Burrows–Wheeler Alignment tool (BWA) using ‘MEM’ command in Paired-End mode and unmapped reads were realigned to the same reference genome in Single-End mode (26). The PCR duplicates were removed using SAMtools with ‘markdup’ (27). After alignment, the single cell copy number values were calculated with Control-FREEC at a resolution of 2M genomic window size (28). The raw copy number matrices obtained from Control-FREEC were processed in R software. First, copy number values less than one were converted to one, and the copy number values greater than four were converted to four. Then, the copy number values were log2 transformed and they were centered to zero by subtracting one for each cell.

Frequencies of CNVs in primary GEJ cancer were calculated based on the copy number value matrix after preliminary transformed. With a window size of 2M, the values greater than 0.3 were regarded as copy number gains while less than -0.3 were regarded as copy number losses. The number of cells with gains or losses in each window were divided by the total primary tumor cell number respective to each patient to eliminate effect of the unbalanced cell numbers among the four patients. Then the four weighted values of each window were summed together and divided by four.

To identify the number of clusters, Euclidean distances between cells were calculated and ‘hclust’ was used to perform hierarchical clustering. To delineate the clonal evolution during invasion, the genomic lineages were inferred, and the data were plotted using R package Timescape (29). The clones were defined according to the clustering results and the proportions of different origins of cells that belonged to each clone were calculated.

Quality control of the raw sequencing data was performed by discarding reads containing adapter contamination, low-quality nucleotides, and unrecognizable nucleotides. Then the clean reads were mapped to the UCSC human reference genome by BWA-MEM algorithm (26). If a read was mapped to multiple positions, BWA allowed to choose the most likely placement (if the multiple placements were most likely, the choice was random). SAMtools was used to sort and merge the mapped reads and Picard was used for duplicate removal. After removing the duplicate reads, recalibrating the base quality scores and local realignment (30), Samtools mpileup and bcftools were used to do variant calling and the somatic SNVs were detected by muTect, which was based on a Bayesian classifier (31). The variants were determined by analyzing the LOD score and the systematic false positive was decreased with the use of filters. The functional annotation of variants, including information in gene transcript annotation databases, such as Consensus CDS, RefSeq, Ensembl and UCSC, and other related databases, such as dbSNP and 1000 Genome, was given using ANNOVAR (32).

The phylogenetic trees were constructed based on the distribution of nonsynonymous SNVs by MEGA5 software. The max-mini branch-and-bound algorithm was used to infer the maximum-parsimony trees (33). We used the average pathway method (34), to calculate branch lengths. The normal tissues were set as the outgroup to gain the consensus tree. Then, the phylogenetic trees were redrawn in Adobe Illustrator software. The lengths of trunk and branch were proportional to the acquired average number of nonsynonymous mutations, with the terminal branches colored in red, the internal branches colored in yellow, and the trunk colored in blue. The angles between branches were chosen just for convenient display. Representative driver SNVs were marked next to the corresponding branches.

In this study, four patients were enrolled. We obtained the primary tumors from all the patients and also two lymph nodes from two of the patients (Pt.3, Pt.4) respectively. For Pt.3 and Pt.4, L1 represented the lymph node anatomically closer to the primary tumor compared to L2. The sample information and clinical characteristics of the four GEJ cancer patients were shown in Table S1 . In order to better understand the genomic heterogeneity and phylogeny of GEJ cancers, we first performed whole exome sequencing (WES) on primary tumor and metastatic lymph node samples to analyze the mutations. The overview of project was shown in Figure 1A .

The somatic SNVs of the primary tumor and lymph nodes were called using the matched normal tissue samples as control. The mutation numbers of the four GEJ cancer patients were shown in Figure 1B . For the primary sites, Pt.1 had the largest number of both nonsynonymous SNVs and total SNVs. As was known, mutations accumulated as age increased and Pt.1 was the oldest among the four patients. The number of SNVs in Pt.2 was the least. The primary tumor in Pt.2 was large, but no clinical metastasis was diagnosed, which might be related to less tumor mutation burden. We also analyzed the mutation spectra, showing similar mutation types among the four patients ( Figure 1C ). Whether in the primary or lymph node sites, the C > T mutation was significantly enriched, consistent with the results of a relatively large set of Chinese GEJ cancer patients (11).

Although there were similar mutation spectra among patients, the numbers of mutations varied among patients or different sampled sites in individuals. Figures 1D, E showed the common SNVs and the primary or lymph node specific SNVs in Pt.3 and Pt.4. The proportion of primary or lymph node specific mutation (the number of primary or lymph node mutations divided by the common mutations) were distinct between Pt.3 (0.40, 2.49) and Pt.4 (0.94, 1.40). Also, numbers of primary tumor and lymph node specific mutations were distinct within the same patient. In Pt.3, the lymph node specific SNVs were significantly more than the primary site (6.31-fold). In Pt.4, the lymph node specific SNVs versus primary site was 1.49-fold.

Next, the recurrent mutated genes among the four patients were analyzed. Figure 1F showed that the most frequently mutated gene was TTN, which occurred in three patients. Eight mutation sites were identified, six of which occurred in Pt.4. The mutation sites of TTN gene were shown in Figure S1 . Most of the mutation sites of TTN occurred in the immunoglobulin I-set domain. Other frequent mutations were in TP53, TG, LRP1B, AAR2, CCDC108, COL7A1, DNAH17, DYNC1H1, PLXNB2, etc., which happened in two patients.

The mutation allele frequency was used to identify the clonal or subclonal mutations between different sampled sites (35). For Pt.1 and Pt.2 that only had primary tumors, the mutation frequencies of nonsynonymous SNVs were shown in Figure S2 . Almost all the SNVs were subclonal, indicated the intra-tumor heterogeneity. Next, we analyzed the mutation allele frequency between primary tumor and matched lymph nodes ( Figure S3 ). Figures S3A, B showed that the mutation subclonal architecture was relatively more conserved between the primary tumor and L1 than primary tumor and L2 in Pt.3. While the subclonal structures between primary tumor and L1 were similar to that between primary tumor and L2 in Pt.4 ( Figures S3C, D ).

For Pt.3 and Pt.4 who had two matched lymph nodes, the mutation frequencies were analyzed between the lymph nodes. Figures 2A, B showed that the mutation frequencies between the two lymph nodes in Pt.4 were more similar than that in Pt.3, suggesting that the subclonal structures of SNVs were relatively more conserved between the two lymph nodes in Pt.4.

The above results revealed the mutation profiles and subclonal architectures of the four GEJ cancer patients. However, the lymph node metastasis patterns of different patients still needed further analysis. The SNV based phylogenetic trees in Pt.3 and Pt.4 with lymph node metastasis were constructed. Figure 2C showed that the two lymph nodes of Pt.3 were in different branches, while one of the lymph nodes and primary tumor were on the same branch, indicating that the two lymph nodes metastasized independently. The distant lymph node (L2) branched from the trunk earlier, indicating that L2 metastasized earlier than L1. However, the phylogeny pattern of Pt.4 showed that the primary site was branched from the trunk, and the two lymph nodes were separated from another branch ( Figure 2D ). Figure 2D indicated that the two lymph nodes had the same genomic ancestor descended from the primary site and metastasized following the lymphatic drainage.

The genes included in COSMIC (Catalogue Of Somatic Mutations In Cancer) database and KEGG pathways in cancer were defined as the inferred driver mutations and marked in the phylogenetic trees. Mutation of TP53 was trunk mutation in both Pt.3 and Pt.4, indicating that TP53 mutations occurred earlier and played an important role in tumor progression and lymph node metastasis. For Pt.3, PIK3CA, JAK2 etc. mutations occurred only in L1, not in L2. This might be related to the phenomenon that the two lymph nodes were separated different branches from the trunk. While in Pt.4, the two lymph nodes were separated from the same branches, simultaneously showing the mutation of POTEF, OSBPL7 and CHRFAM7A.

The SNV-phylogenetic analysis revealed distinct phylogeny patterns of lymph node metastasis in the two patients. The lymph nodes could separate independently from the primary tumor at different time points, and also metastasize station by station following the lymphatic drainage.

In the above results, we found that SNVs were heterogeneous between primary tumor and matched metastatic lymph nodes. Next, we analyzed the heterogeneity of CNVs at single cell level, which could not be found by bulk sequencing. 210 single cells from the four patients passed quality control and the single cell information was shown in Table S3 . First, we analyzed the heterogeneity of Pt.1 and Pt.2, who only had the primary tumors collected. The heatmap in Figure 3A showed the global copy number profiles in Pt.1. The data identified the CNV events to be clonal, due to the similar CNV patterns of all the single cells. Figure 3B showed the normalized CNVs, including the amplification in Chr8, 13, 20 and deletion in Chr14, 15. The heatmap in Figure 3C showed the global CNV profiles in Pt.2. The data identified two major clones in the primary tumor. Figure 3D showed the CNV patterns of the two subclones. The shared CNVs included the amplification in Chr7, 13, 20 and deletion in Chr5q. There were also distinct CNV regions between the two subclones, such as a greater degree of deletion on Chr4 and 14 in subclone 2 compared to subclone 1.

Next, we analyzed the intra-patient heterogeneity in Pt.3 and Pt.4, who had the single cell information of the primary tumor and lymph nodes. Figures 4A , 5A showed the global CNV profiles of the individual cells in Pt.3 and Pt.4 respectively. The data identified four major clones in each patient. Figure 4B showed the CNV patterns of the four subclones in Pt.3. The shared CNVs included the amplification in Chr7p,20 and deletion in Chr5q,21. The subclone specific CNV events were also observed. Figure 5B showed the CNV patterns of the four subclones in Pt.4. The differences of CNV profiles between subclones revealed obvious intra-patient heterogeneity. This heterogeneity not only existed in the primary tumor or lymph nodes, but also existed between the primary tumor and lymph nodes.

It was reported that CNV events in primary tumor often also occurred in single or multiple matched lymph nodes, but a few CNVs in the lymph nodes were not found in primary tumor (36). This indicated that the degree of heterogeneity between primary tumor and lymph nodes were different. Figure 4A revealed that in Pt.3 the primary tumor consisted of three subclones, while the L1 and L2 lymph nodes consisted of one or two subclones. Figure 5A also revealed that in Pt.4 the primary tumor consisted of three subclones, while the L1 and L2 lymph nodes consisted of one or two subclones. The results indicated that the heterogeneity of primary tumor was higher than that of lymph nodes in the same patient.

Unlike the evolutionary model of mutation, which evolved gradually, resulting in extensive clonal diversity, the CNVs occurred in the early stage of tumor evolution and were highly stable with the clonal expansion of tumor (13). Our above mutation data showed that the phylogeny patterns of lymph node metastasis in Pt.3 and Pt.4 were distinctive ( Figure 2 ). Next, we analyzed the modes of clonal evolution using CNV data. The clonal frequencies could be analyzed to visualize the clonal changes with tumor metastasis (29). Based on previous SNV data, the two lymph nodes in Pt.3 originated from the primary tumor independently, so we analyzed the proportion changes of each subclone from primary to L1 and L2 lymph nodes separately.

In Pt.3, there were three subclones in primary tumor, including clone 2, clone 3 and clone 4 with percentages of 16%, 58% and 26% ( Figure 4C ). In the lymph nodes, the clone 3, which was the main proportion of primary tumors, and clone 2 disappeared. However, a new subclone (clone 1) appeared and accounted for 100% in L1 lymph node. Two subclones (clone 1 and clone 4) were found in L2 lymph node and clone 1 accounted for 81%. The main differences of CNVs between clone 3 and clone 1 were clone 1-specific deletion in Chr1p and 18. These CNVs might be related to lymph node metastasis in Pt.3.

In Pt.4, there were three subclones in the primary tumor, including clone 1, clone 2 and clone 3 with percentages of 17%, 65% and 17% ( Figure 5C ). The clone 1 and clone 2 disappeared in the lymph nodes. However, a new subclone (clone 4) appeared in the lymph nodes and accounted for 100% in L1 lymph node and 94% in L2 lymph node. The main differences of CNVs between clone 4 and other subclones were that the obvious deletion of Chr4 and 6. These CNVs might be related to lymph node metastasis in Pt.4.

The clonal evolution analysis in Pt.3 and Pt.4 showed that the proportion of different subclones changed significantly with the metastasis of GEJ cancer. The findings further revealed the molecular mechanism of tumor metastasis in GEJ cancer.