Between January 2013 to December 2016, 18 non-metastatic NSCLC (Nome_LC) patients with tumor size over 3 cm, and 18 T1 NSCLC patients with lymph node metastases (Meta_LC) from Jiangsu Cancer Hospital were enrolled for this study. A total of 59 NSCLC tissues including primary lesion and metastatic lymph nodes were extracted from the 36 patients, fixed in formalin and embedded in paraffin (FFPE) for somatic mutations analyses. The pathological FFPE findings were validated using hematoxylin and eosin (HE) staining. The age, gender, TNM stage, tumor histology and smoking history of the patient were all recorded.

An NGS approach was performed on genomic DNA isolated from FFPE surgically resected NSCLC tissue samples. Tissues smaller than 1 mm or with less than 20% tumor cells were all excluded. A targeted pan‐cancer panel with 1021 cancer related genes was conducted and was performed on the Hiseq NGS platforms (Illumina Inc., San Diego, CA, United States) using of base substitutions, indel, copy number variations (CNV) and DNA arrangement analysis. Targeted capture sequencing revealed a median effective depth of coverage of 671× in the tissue samples. The detailed genes for the 1021 cancer-gene panel were provided in Supplementary Table S1.

The tumor mutation burden (TMB) of a tumor a tumor was calculated based on the number of non-synonymous somatic mutations (single nucleotide variants and small insertions/deletions) per mega-base in coding regions. All germline variants were filtered using paired blood controls.

Phylogenetic analysis of tumor genes sequences was analyzed using the ape package. For each patient, an optimal phylogenetic tree was inferred using the maximum parsimony method. Candidate driver mutations were mapped to the trunk and branches of each phylogenetic tree to illustrate the molecular processes underlying NSCLC metastasis.

Subclones generated from SNVs were clustered sing PyClone. The optimal tree was generated using the clustering outputs. The analysis was performed using iterative version of citup. The fish plot and phylogenetic tree of each patient were constructed using timescape.

Categorical variables were expressed using number and proportions. Somatic mutation data of NSCLC patients from TCGA database were downloaded from cBioPortal [8]. The relationships between groups were analyzed using Pearson’s chi-square or Fisher’s exact test, based on the Yates continuity correction coefficient, with p < 0.05 considered statistically significant. Data was analyzed using SPSS22.0 software (SPSS, Inc., Chicago, IL, United States).

The baseline characteristics of the 36 NSCLC patients included in this study are summarized in Table 1. Overall, 59 NSCLC tissues were extracted and analyzed. Of the 36, 27 patients (75%) were male, with majority (83.33%) above 65 years old. The patients were divided into two groups based on the absence or presence of lymph node metastasis. We found 18 patients with one or two lymph nodes metastases (Meta_LC group). The median tumor size in these patients was 2.5 cm (range 1.0–3.0 cm). The rest 18 patients (50%), had no lymph node metastasis (Nome_LC group). Nine tissues were selected from each group of lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) for further analyses. The combined median tumor size of Meta_LC group was 2.5 cm (range 1.0–3.0), almost three times the size of the smallest tumor in the Nome_LC group (7.1 cm, range from 5.0 to 12.0). There was no distant metastasis in any of the cases (M stage: M0).

Next-generation sequencing of the 59 samples uncovered a pan-cancer panel of 1021 tumorigenesis related genes. The average sequencing depth of DNA extracted from the aforementioned tissues was 671×[9]. The somatic mutations investigated included small nucleotide variants (SNVs), copy number variants (CNVs), short insertions or deletions (indels) and structure variations (SVs). Overall, we identified 486 somatic mutations, with the majority in TP53 (58.3%), EGFR (21.7%), MLL3 (13.3%), LRP1B (13.3%) and TERT (13.3%) genes (Figure 1A). We downloaded sequence data for 68 LUAD and 70 LUSC patients from The Cancer Genome Atlas (TCGA) for further analyses (Supplementary Figure S1A). The top five most common alterations were not quite the same as genes identified in our study. For Nome_LC group, mutations were detected in 126 genes. Comparatively, there were 51 gene variations in the Meta_LC group, 26 genes of which were also detected in the LN group. Notably, mutations in 11 genes were unique to the LN group alone (Figure 1B). Mutation frequency of common driver genes was quite different between Nome_LC and Meta_LC groups, although it was not statistically significant (Figure 1C).

Intriguingly, there was a clear distinction in the frequency and type of mutant genes between Nome_LC and Meta_LC groups. For instance, for Nome_LC group, the majority of mutations occurred in the TP53 (66.7%), PIK3CA (27.8%), LRP1B (22.2%), CTNNB1 (22.2%), and EGFR (16.7%) genes. However, TP53 (72.2%), EGFR (33.3%), LRP1B (16.7%), MLL3 (16.7%) and TERT (16.7%) were the top most mutated genes in the Meta_LC group (Figure 1D). EGFR (33.3% vs. 15.8%) and TERT (16.7% vs. 5.2%) were the most mutated genes in the Meta_LC group (Supplementary Table S2). In LN group, TP53 (43.5%), EGFR (17.4%), TERT (17.4%), ERBB2 (13.0%) and FAT1 (13.0%) were the top most mutated genes (Figure 1D). Notable, the difference in the most common mutations between Meta_LC and LN group was very small (Supplementary Table S3).

We further compared the distribution of the most frequently altered genes in Meta_LC, Nome_LC and LN groups. The highest gene mutations were observed in the Meta-ca group. Interestingly, there was a significant difference in the distribution of altered genes between Nome_LC and LN group (Supplementary Figure S1C).

On each chromosome, the gene was marked when its mutation frequency was greater than 50% and we defined it as having a high mutation frequency. With regard to high-level frequency of mutated genes among the three subgroups, Meta_LC group exhibited two high-level frequency mutations (EGFR, CDKN2A), whereas in LN group, high-level frequency gene mutation was only observed in one gene (NCOR1). However, high-level frequency mutations were observed in ten genes (EPHA3, CDKN2A, ATM, MAPK3, TP53, NOTCH3, STK11, CIC, FAM123B, RBM10) in the Nome_LC group (Figure 1E). Overall, patients with large primary tumor size and those without lymphatic metastasis displayed more high-level gene mutations. Even so, precise mechanism underlying the phenomenon needs further investigation.

Further analyses reveled that patients with higher Tumor Mutation Burden (TMB/Tumor Mutation Load (TML)) responded better to tumor immunotherapy [10]. Compared to patients in Meta_LC and LN groups, those in Nome_LC group presented with higher mutation burden. However, no significance difference in median TMB was seen between primary and metastatic tumors (Supplementary Figure S1B).

Sequence data for 562 LUAD and 470 LUSC patients (1032) was downloaded from the TCGA database. Additional validation data for 5597 patients was downloaded from the Memorial Sloan Kettering Cancer Center (MSK) database. As shown in Supplementary Figure S2, about 90% of high frequency mutated genes both in LUAD (Supplementary Figure S2A) and LUSC (Supplementary Figure S2B) detected in this study were indentified in cancer tissues of patients in the TCGA and MSK database, which means the results obtained by panel sequencing were reliable.

Further comparisons in genetic aberrations between LUAD and LUSC of our data were made. We identified 173 somatic mutations in 32 lung adenocarcinoma tissues (Figure 2A). The top most altered genes included TP53, EGFR, MLL3, APC, ERBB2, LRP1B, BAP1, BRAF, BRCA2 and CTNNB1. Nome_LC and Meta_LC groups only shared eight mutated genes (Figure 2B). In addition, there was no significance difference in mutation frequency between the two groups (Supplementary Table S4). In 27 lung squamous carcinoma tissues, we identified 302 somatic mutations (Figure 2C). The top five most commonly altered genes included TP53, FAT1, PIK3CA, TERT and EGFR (Supplementary Table S5). Meta_LC and Nome_LC groups shared 21 mutated genes, conflicting findings of the TCGA analyses (Supplementary Figures S2C,D). Interestingly, the rate of mutations in squamous cell carcinoma tissues was significantly high, relative to adenocarcinoma tissues. The average number of mutations per LUSC and LUAD tissue was 8.1 and 5.4, respectively (Figures 2A,C). These findings demonstrate the marked mutational distinction between LUSC and LUAD tissues. However though, LUSC exhibits greater gene mutations.

Given that small primary tumor (T1) may present with lymph node metastasis, and relatively large primary tumors (>T2) can present without lymph node metastasis, we sought to identify the genetic mutations associated with lymph node metastasis or proliferation of primary tumor cells.

The differences of mutated genes were analyzed in Meta_LC and Nome_LC group. We found that mutations in CTNNB1 was significantly (p < 0.05) different and only presented in Nome_LC group in LUAD (Figure 3A). Further analysis of 27 LUSC tissues only identified mutations in ARID1A genes (Figure 3B) among patients in the Nome_LC group. Analysis of data downloaded in the TCGA database identified significant mutations in KCNH7, UNC9, ADAMTLS3, HCN1, AHNAK2, EPHA7, MXRA5, CTNNA2 and SEMA5A (p < 0.05) in the two groups of LUAD (Figure 3C). and 17 genes (E4F1, EPYC, FGF12, etc.) were detected in TCGA data (Figure 3D). In general, findings of TCGA data analyses were substantially different from those generated from analysis of our data. This discrepancy may have arisen from numerous sources.

To further understand the evolution of lymph node metastasis, we compared genetic aberrations in primary tumors and metastatic lymph node in the Meta_LC group. The trunks of phylogenetic trees represented genes with mutations both in primary tumors and all lymph nodes. Branches represent mutated genes with lymph nodes only (Figure 4). In some cases, we observed high similarity in genetic variation between metastatic lymph nodes (Figure 4M) and their paired primary tumors. This suggest that metastatic lymph nodes “drifted” from the primary tumor relatively late in the evolutionary process. However, other genetic alterations were highly distinct between groups (Figure 4N), demonstrating early “drift” between evolutionary groups. Overall, we identified two distinct evolutionary models both in LUAD and LUSC. These results suggest that metastatic lymph nodes followed the linear progression model and parallel progression model. In the cases of the primary tumor with two metastatic lymph nodes, there were also different evolutionary patterns (Figures 4A,G,I–K). In one pattern, two metastatic lymph nodes shared exactly the same branch (Figure 4A). In another pattern, two metastatic lymph nodes shared a trunk but have distinct branches (Figures 4I,J). And in the other pattern, similarity was observed between one metastatic lymph node and its paired primary tumor, the other metastatic lymph node however, was quite different between them (Figures 4G,K).

Clonal analysis of altered genes in 14 NSCLC lymphatic metastases revealed the hypothetical subclonal development process and cancer evolution. The mutant gene with the maximum load was regarded as having a cell proportion of 1 (100%), and the other genes decreased in turn. The mutation load was negatively proportional to the occurrence time. The evolution patterns of subclones differed from patient to patient, characteristic of intratumor heterogeneity (Figure 5). It’s worth noting that several LUAD tissues displayed fewer or slightly more mutated (Figures 5A,B,E,I). This phenomenon corresponded with low level TMB. However, there were greater mutations in most LUSC, which corresponded with high TMB (Figures 5H,L,M).