Multiregional Sequencing Reveals Genomic Alterations and Clonal Dynamics in Primary Malignant Melanoma of the Esophagus

Mucosal melanoma is most often the result of metastases from primary cutaneous melanoma (PCM), whereas primary mucosal melanoma (PMM) is an exceedingly rare disease that behaves more aggressively and has a poorer prognosis (1). The molecular and pathological characteristics of cutaneous melanoma have been well studied, and more than 1,000 PCM exomes have been sequenced (2, 3). Key cellular pathways, such as CDKN2A/CDK4/CCND1/RB1, MAPK and PI3k/AKT, and significant genes, such as BRAF, NRAS, TERT, TP53, NF1, and RB1, are considered to play critical roles in PCM tumorigenesis (4–7). Recently, two studies with very limited sample sizes have claimed that distinct mechanisms drive PMM and PCM (8, 9). Moreover, melanoma subtypes affecting different body regions exhibit different genetic features, even those of the same histological type (2, 10, 11). Despite these findings, there is a paucity of data elucidating the etiopathogenesis and predictive factors, as well as the effectiveness of therapeutic options, for these diseases.

Primary malignant melanoma of the esophagus (PMME) is a rare gastrointestinal (GI) mucosal melanoma, representing 0.1% to 0.2% of all esophageal malignancies, with only 339 cases reported by 2016 (12–14). It is reported that PMME has a tendency for vertical growth within the esophageal wall (15). Melanosis is another important feature favoring PMME over a metastatic melanoma to the esophagus (15). However, the molecular background of PMME is largely unknown because of its rarity. Recent reports have included only small series or case studies with analyses of a single or few molecular aberrations. The BRAF, NRAS, KRAS, and c-Kit genes have been reported to be the most commonly altered genes in PMME (16, 17). Although evidence indicates that GI mucosal melanoma arises primarily at certain areas within the GI system, it is still difficult to distinguish PMM from metastatic melanoma in the GI tract that may be derived from an unknown or regressed cutaneous primary tumor. Tumorigenesis and metastatic mechanisms, as well as other clinical performance of PMME, are lack of scientific interpretation due to the vacancy of studies.

Here, we applied exome sequencing and chromosomal aberration analysis to study multiple spatially separated resected samples from a PMME patient, including samples of the primary tumor, metastatic lymph nodes, normal esophageal mucosal tissues, and nonneoplastic melanosis lesion. The presence of a precursor lesion or melanosis could be a marker for GI mucosal melanoma (18); however, this notion is still under debate, and more evidence is needed. In general, this PMME exhibited distinct genetic features, with both similarities to and differences from those of cutaneous melanoma. Significant driver genes in PCM, such as BRAF and CDKN2A, were determined to be potential driver genes in this PMME, whereas infrequently altered genes in PCM, such as KRAS, also played a role in PMME tumorigenesis. By characterizing the genetic basis and intratumor heterogeneity (ITH), we found that the genome of this PMME fit a branching evolution model, as ancestor clones evolved along different trajectories and diversified into primary tumor subclones and early and late metastatic subclones. Interestingly, the mutational signature switched from a C>T to a C>A transition with tumor progression, indicating that different inducers promoted each trunk and branch mutation during the mutational processes.

Comprehensively resolving the genetic architecture and temporal dynamics of the mutagenic processes in PMME could help us better understand the genetic basis of this rare malignancy. With more PMME cancer genomes being deciphered, advanced diagnostics and therapeutics could be developed in the near future.

A 59-year-old Chinese female was admitted to the Thorax Surgical II Department of Beijing Cancer Hospital in October 2013 with a 3-month history of dysphagia. Physical examination and CT scan showed no evidence of skin melanoma lesions. Esophagogastroduodenoscopy revealed a 3.0 × 2.5-cm polypoidal tumor with black pigmentation; 3 cm below the bottom edge of the main tumor, a flat, 0.5 × 1.0-cm mucosal melanosis lesion was observed. The patient did not receive any treatment before surgery. The patient was informed and signed an informed consent form before surgery. This study was approved by the Medical Ethics Committee of Peking University Cancer Hospital (#2013KT31) and conducted in accordance with the Declaration of Helsinki.

The resected esophageal specimen was photographed and marked with a vertical ruler. Then, spatially separated samples, including four primary tumor samples (T1–T4), five normal mucosal samples (N1–N5), and one mucosal melanosis sample (ME1), were collected discretely. Lymph node samples (LN1–LN4) were isolated from surgical specimen tissue blocks, and a peripheral blood sample was acquired after surgery (Fig. 1A and Supplementary Fig. S1). Each tumor and normal mucosal sample (10–20 mg) was cut into thirds; one piece was used for pathological analysis, the second was used for genomic DNA extraction, and the third was kept in a storage solution (RNAlater, Thermo Fisher Scientific) as an extra piece. Experienced local pathologists performed a histopathological review of all primary tumors, including IHC for detection of HMB-45, Melan-A, and S-100 protein expressions (Supplementary Fig. S2). The ME1 sample had a jelly-like texture, which may have led to its failure in the formalin-fixed, paraffin-embedded (FFPE) process; thus, no further pathological examination was performed on this sample. Genomic DNA was isolated using a DNeasy Blood and Tissue Kit (QIAGEN) or a QIAamp DNA FFPE Tissue Kit (QIAGEN).

Exome capture was performed using 3 to 6 μg of genomic DNA per sample with an Agilent SureSelect Human All Exon V4 Kit according to the manufacturer's instructions. Briefly, genomic DNA was fragmented with an ultrasonicator (Covaris) to obtain fragments that were primarily between 250 and 300 bp in length. These fragments were amplified using an NEBNext Ultra DNA Library Prep Kit (NEB) and then hybridized with an Illumina TruSeq Capture Kit for enrichment. Each captured library was loaded onto an Illumina HiSeq2000 platform, and paired-end sequencing was performed to the desired median sequencing depth (>200 ×). Raw image files were processed using HCS1.4.8 for base calling with the default parameters.

The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive in BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under accession number CRA000346, that are publicly accessible at http://bigd.big.ac.cn/gsa.

The sequencing reads were aligned to a human reference genome (hg19) using BWA (19). Mutect (v1.16) was used for somatic mutation calling. The following filtering criteria were applied to the mutation profiles produced with Mutect to obtain highly probable somatic nucleotide variants (SNVs): (i) only SNVs marked with a “keep” status were considered; (ii) minimum depths of 10 × and 14 × were required for the normal and tumor samples, respectively; and (iii) a minimum of four reads were required to support the mutant allele. Small indels were identified using Pindel version 0.2.4 (20) in paired tumor–normal mode. Highly probable indels were identified using an in-house pipeline as follows: (i) all identified indels were required to have a minimum coverage of 20 reads in both the tumor and normal samples; and (ii) a minimum of 5 independent reads was required to support the indel event in the tumor sample, with no evidence of such an event in the paired normal sample. The qualified variants were annotated with oncotator (21).

To identify possible driver events, we classified the altered genes in this PMME according to reports in the literature, as well as databases. Mutated genes were classified into one of six tiers: (i) Tier A: cutaneous melanoma-related gene (2); (ii) Tier B: probability score for oncogene (OG) or tumor suppressor gene (TSG) > 0.75, as indicated by TUSON Explore (22); (iii) Tier C: suggested TSG or OG by Vogelstein and colleagues (23); (iv) Tier D: suggested significant cancer gene by Lawrence and colleagues (24); (v) Tier E: cancer-related gene in the Cancer Gene census catalog (http://cancer.sanger.ac.uk/census); and (vi) Tier F: recurrent gene mutation in mucosal melanoma (8). SIFT (25), PolyPhen-2 (26), and MutationAssessor (27) were used to predict the potential biological function changes of each missense mutation. If any mutation received a rating of "Damaging" in SIFT, "Damage"/"Possible Damage" in Polyphen-2, or "Medium"/"High" in MutationAssessor, it was considered to have a high functional impact (HiFI). Nonsense and frameshift mutations, as well as splicing mutations, were also regarded as HiFI alterations. Non-HiFI mutations were not considered to contribute to PMME tumorigenesis or invasion.

Whole-genome SNP array for each samples was performed on Illumina Human OmniZhongHua BeadChips (Illumina). Raw intensity values were processed to obtain a normalized B allele frequency (BAF) and a Log R ratio (LRR) for each probe using the Genome Studio Software V2011.1. LRR values were segmented by the ASCAT algorithm (28). Tumor and matched normal samples were used in pairs to obtain somatic alterations.

Somatic copy-number alterations (SCNA) were predicted with ADTEx (29) using both depth of coverage ratios and B allele frequencies (BAF) calculated from the WES data. Briefly, ADTEx consists of two hidden Markov models (HMM) to predict SCNAs and genotypes in WES data. One HMM is used to predict CNAs, which, in combination with BAFs, can be used to estimate tumor ploidy and predict absolute copy numbers. The second HMM is used to predict the zygosity or genotype of each CNA segment. In this study, ADTEx v2.0 was run with the BAM and BAF inputs using the default parameters.

ABSOLUTE (v1.0) was used to predict the tumor content, ploidy, and absolute copy number based on the copy-number profile and somatic SNVs for each tumor. The cancer cell fraction (CCF) of each mutation was subsequently estimated according to the mutant VAFs after adjustment for the tumor purities and local absolute copy number. Private, branch, trunk mutations were determined by the extent of shared mutations among the different samples.

In principle, a phylogenetic tree was generated according to the numbers of shared mutations among the samples. Of the sequenced samples, three were considered to be normal tissue samples (N1, N5, and LN4) with few mutations, and another three (LN3, ME1, and N3) had a low tumor content (36%, 11%, and 7%, respectively). The four primary tumor samples (T1–T4) possessed a high tumor content of over 80%, whereas the other two lymph node samples (LN1 and LN2) had a tumor cell proportion of >50%. Utilizing low-tumor-content samples can introduce bias into the phylogenetic tree. Therefore, we selected SNVs from six high-tumor-content samples to create the topological structure that reflected the main evolutionary trajectory of this tumor. Then, we added the remaining samples to the prototype tree according to the numbers of SNVs they shared with the existing taxon. Because ME1 and N3 had a relatively low tumor content, we were not able to estimate the precise divergence time when they derived from the founding clone; thus, we suggested that they should be located at the beginning of the trunk of the tree (Fig. 1B). During tree construction, we also considered tumor content, copy-number changes, and intermixed clones that might have led to the misinterpretation of our findings. By calculating branch/trunk lengths inferred from the mutation counts, the final trees were drawn manually.

Whole exon-capture multiregion sequencing (WES) was performed on DNA from 12 samples derived from esophagectomy specimens, including four samples obtained from primary tumor dissection (T1–T4), three normal mucosal samples (N1, N3, and N5), one melanosis sample (ME1), and four metastatic lymph node samples (LN1–LN4). Germline control DNA was extracted from the peripheral blood. From 13 sequenced samples, high-coverage sequencing (mean, 200 ×) data were obtained. Overall, 386 nucleotide substitutions and 134 indels (insertions and deletions) were detected in at least one sample (Supplementary Tables S1 and S2). Fifty single nucleotide variants (SNV) were randomly selected and validated in four primary tumor samples (T1–T4) by nucleotide mass spectrum (Sequenom). The genotypes were successfully validated for 188 (94%) of the 200 (50 mutations per sample) tested loci. Tumor purity of the primary tumor was estimated to be nearly 80% using ABSOLUTE software according to the exome sequencing data (Table 1), indicating that our data were sufficiently robust to detect subclonal mutations with low variant allele frequencies (VAF).

To compare the differences in the mutation spectra between PMME and PCM, we acquired genomic data for PCM from TCGA database (http://www.cbioportal.org/). A total of 281,339 SNVs from 341 PCMs were used in this analysis. We also utilized WES data from 81 squamous esophageal carcinomas (derived from unpublished data obtained in our laboratory) to compare similarities in mutational signatures among ESCC, PCM, and PMME. As reported previously, the PCM genome exhibited a specific predominance of C>T transitions at TpC dinucleotides, which were probably associated with exposure to ultraviolet light. The ESCC genome was characterized by a mutational process associated with an APOBEC-mediated and aging signature; thus, the mutations tended to be C>T or C>G alterations. In this PMME genome, we found that 46.1% (178/386) of the somatic mutations were C>A transitions, 41% (73/178) of which occurred in the CpCpA trinucleotide context (Fig. 2A and B). However, no similar signature was identified in the COSMIC data, suggesting that the CpCpA signature found in this PMME genome may be caused by an unknown mutational process that needs to be further investigated. In summary, although these three types of tumors shared similar anatomical or pathological characteristics, they underwent distinct mutational processes, indicating that different molecular mechanisms drive the tumorigenesis and progression of each tumor type.

SCNAs were investigated in this PMME. At the arm level, we observed copy-number loss of 1p, 3p, 4p, 11p, 11q, and 14q and copy-number gain of 1q, 6p, 8p, and 8q (Fig. 3A). To identify focal amplifications or deletions in these samples, we applied GISITC to determine copy numbers for these 13 regional samples and found 23 amplification and 21 deletion peaks across the genome, some of which were associated with potential oncogenic genes, such as EIF3E, KRAS, and SOX5, and the tumor suppressor genes PTEN and CDKN2A. Interestingly, we observed characteristics of chromothripsis in chromosome 3p in the primary tumor samples (T1–T4) and two metastasis lymph nodes (LN1 and LN2; Fig. 3B). These observations were based on the following principles: (i) only one chromosome was affected; (ii) distinctive oscillation between two copy-number states arose across the affected regions; and (iii) within the regions with “higher” copy-number states, retention of heterozygosity was observed.

To assess ITH, we examined 235 nonsynonymous point mutations from the primary tumor, including the T1–T4 regions, to determine the regional mutation distribution. Only 20% (47/235) of the mutations were detectable across all four sampled tumor regions, whereas 54.5% (128/235) of the mutations existed in a single region. Among all of the sequenced spatial samples, N1, N5, and LN4 were confirmed to be in the normal mucosal region due to their infrequent mutations, and N3 and ME1 showed pretumor clonogenic potential because they both shared pathogenic mutations with other samples of high tumor content, suggesting that they were clonally related. As the tumor content was relatively low in these two samples, the number of SNVs was not as high as that in the primary tumor regions. T1–T4 and LN1–LN3 were considered to be tumor regions, and we found an average of 136 mutations per sample, with the number of SNVs ranging from 98 to 188, in regions with a tumor content of over 50% (Table 1).

To characterize the ITH and evolutionary trajectories in the PMEE genome, a phylogenetic tree was constructed based on the extent of shared somatic mutations identified in each of the tumor cell–infected regions (T1–T4, LN1–LN3, ME1, and N3), and this tree was used to determine the clonal relatedness of the tumor subclones that coexisted in this PMME. The results showed that this PMME tumor had a typical branched rather than linear tumor evolution. Subclones of the primary tumor and other regions shared genetic mutations inherited from a common ancestor. Allopatric clones evolved and followed distinct evolutionary trajectories. ME1 and N3 remained at the original phylogenetic position and slowly developed alterations, whereas the T1–T4 sample region acquired the necessary genetic alterations to form the primary tumor clone. The original tumor clone continued to progress and diversified into two major branches: clade1 and clade2 (Fig. 1B). Interestingly, we observed that clade1 seemed to undergo a trunk mutation accumulation stage, during which time dozens of mutations emerged in this branch, which then shaped the main clonal constituents in clade1 that evolved into new subclones with metastatic ability. Notably, we found that the tumor sample T4 shared mutations with both clades, suggesting that the subclones had originated from clade1 and clade2 and that they had broken the spatial barriers of the solid tumor and intermixed in the T4 region. We further divided T4 into T4a and T4b according to their respective subclones, and we then compared the VAFs of the shared mutations between the spatial samples in clade1 and clade2. Our results showed that nearly 20% and 60% of the tumor subclones of T4 came from clade1 and clade2, respectively. Additionally, the VAFs of the mutations in T4a and T4b were lower than those of the mutations in the other regions in each clade (Supplementary Fig. S3A), which further highlights the fact that T4 was derived from two progenitors.

At the copy-number level, this PMME also exhibited ITH. Copy-number loss of 1p, 3p, 11p, 11q, and 14q and gain of 1q, 8p, and 8q were ubiquitous and conserved across all tumor clones, and these SCNAs were considered to have occurred as early evolutionary events (Fig. 3A and C). Focal amplification of SOX5 and EIF3E, as well as homogeneous deletion of CDKN2A, was in the trunk clones. Telomere-bounded 1q and 2q amplification, internal 4q deletion, neutral LOH of 7p, chr13 arm-level deletion, and 12q gain were specific occurred in clade1; telomere-bounded 4q deletion, 7p amplification, chr10 arm-level deletion (Supplementary Fig. S3B), and chr18 arm-level deletion were identified only in clade2 (Fig. 3A and Supplementary Fig. S4A–S4C). Furthermore, we found that clade-specific SCNAs intermixed in T4, such as deletions involving 4q, telomere-bounded amplification of 2q, chromosome 10 and 18 deletion, further confirming that T4 is intermixed clones from two different clades. For example, 10q deletion involving PTEN (chr10q23.31) was detected in clade2 but absence in clade1, whereas PTEN deletion of T4 showed a lower VAF (Supplementary Fig. S3B). The phylogenetic relationships between each of the spatial regions reflected by the SCNAs were highly consistent with the mutation tree.

Potential parallel evolution was observed in this PMME at copy-number level, although no evidence of parallel evolution was found at the mutational level. A mirrored subclonal allelic imbalance was observed in chromosome 4. If the maternal allele is gained or lost in a subclone in one region, yet the paternal allele is gained or lost in a different subclone in another region, it will result in a mirrored subclonal allelic imbalance profile (30). Different SCNAs of chr4 deletion that affected maternal allele of chr4:82–143M and paternal allele of chr4:93M–telomere were observed in clade1 and clade2, respectively (Fig. 4), suggesting parallel evolution in branched stage. Moreover, even inside the independent clades, such chromosomal instability persisted. In clade1, 27% subclones of T1 obtained a privately chr4 arm-level deletion, which coexisted with the clade1 common chr4 aberration (chr4:82–143M); in clade2, 46% tumor subclones of T2 shared clade-specific chr4 loss (93M–telomere) with T3 in paternal allele, whereas the remaining 54% subclones of T2 carried chr4:103M–telomere loss in the maternal allele (Fig. 4).

To further explore the dynamics of the mutational processes shaping this PMME genome over time, the temporal spectra of the somatic mutations in the different evolutionary stages were analyzed. The PMME genome exhibited a significant increase in the proportion of C>A transversions in branch mutations compared with trunk mutations (P = 0.006; Fig. 2B), which was accompanied by a decrease in C>T transitions. The private mutations in ME1 and N3 presented with similar signatures as the founder mutation (Fig. 2B).

To determine whether specific driver genes carried predominantly trunk or branch mutations, we identified potential driver mutations according to the known cancer gene catalogs. A total of 22 altered genes were present in the catalogs, and each of them carried a single mutation, except CDH13 (Fig. 4). Among the 26 potential driver mutations identified, 17 were predicted to be pathogenic by SIFT/Polyphen2/MutationAssessor. Potential functional trunk mutations in BRAF, KRAS, DNAH7, and CDH13 were detected in all of the primary tumor regions; therefore, these mutations were considered to have occurred during the early stages of tumorigenesis. In addition, branch mutations in JAK, ARHGEF4, and UMOD were observed. The JAK mutations belonged to clade2, and the mutations in the other two genes were specific to clade1. Mutations in PIK3CA, GPRASP1, and ARHGAP26 occurred exclusively in T1 and lymph node metastatic regions, whereas their cancer cell fraction (CCF) values were increased in the metastatic regions, which suggested that the tumor cell populations that harbored these mutations might have experienced subclonal expansion after exposure to the metastatic lymph node environment. Activating mutations in TRPS1, MAPK3K13, AHNAK, DNAH3, PCDHGA3, PRUNE2, and TAS1R2 were private and distributed throughout different spatial regions (Fig. 1B). However, their CCFs were much lower than those of the founding and branch mutations, as they appeared to be intraregional subclonal mutations (Fig. 5).

The genomic landscape of cutaneous melanoma has been well described through the use of next-generation sequencing (NGS). New diagnostic and prognostic biomarkers, as well as therapeutic approaches, are being developed based on known germline and somatic genetic alterations in PCM. These advances will improve the treatment of PCM patients. However, the genetic nature of PMM still has not been elucidated; thus, patients with PMM, which is an important subtype of melanoma, may not benefit from the advances facilitated by the new technologies. Therefore, sequencing efforts must be continued, and particular attention should be paid to rare melanoma subtypes to comprehensively characterize the somatic mutational landscape of melanoma. In this study, we sequenced a rare primary esophageal melanoma using a multiregion exome sequencing strategy to characterize the genetic profile, ITH, and clonal evolution of this PMME.

Our data revealed that the four primary tumor samples harbored a total of 110 coding region nsSNVs, whereas each tumor site carried an average of 55 nsSNVs (T1: 67, T2: 48, T3: 38, and T4: 67). The real mutational burden was underestimated by one-fold when only one tumor sample was sequenced, which suggests that high ITH existed in this PMME genome. The high somatic mutational burden in melanoma has been previously linked to favorable outcomes in immune-based therapeutic approaches (31); thus, it is necessary to consider the influence of ITH in the PMME genome when evaluating mutational burden, particularly in clinical practice.

To better interpret the genomic ITH and tumor evolutionary history in this PMME, we selected 23 important mutations from among all of the altered genes as representative genetic evolutionary events. The candidate driver genes for mucosal melanoma are currently unclear; thus, we used pan-cancer and cutaneous melanoma catalogs, as well as data from nine reported mucosal melanoma genomes, to predict putative driver genes. Genomic ITH of this PMME was observed for both somatic mutations and chromosomal imbalances. Among the 22 selected putative driver mutations, only four were clonal events; the other mutations, including a canonical hotspot mutation, PIK3CA^E545A, were restricted to tumor subclones. In copy-number level, ITH was also significant. Some SCNVs involving important cancer genes such as homogeneous deletion of CDKN2A, and high-level amplification of KRAS and SOX5 (Fig. 3C) were ubiquitous across the clade 1 and clade2. Each evolution branch had its own specific SCNAs, such as PTEN focal deletion that exclusively occurred in clade2 and chr13 arm-level deletion that only appeared in clade1. The ITH status could have important implications for targeted therapeutic strategies. Strikingly, we found that mutations in genes involved in the PI3K–AKT–mTOR pathway, such as PIK3CA, tended to be present in subclones, and some reports have suggested that PIK3CA mutations often lead to subclonal expansion, whereas mutations in genes associated with RAS–MEK or cyclin-dependent kinase often result in founding clones. These results are consistent with previous research (32). Therefore, our findings suggest that the current methods utilized to identify cancer genes are likely biased to detect clonal drivers that occur at a higher frequency in single samples; thus, it is necessary to use a more efficient strategy to obtain a full list of clinically actionable genes, such as collecting spatial samples from the same individual and pooling DNA for further analysis.

A phylogenetic tree was constructed for this PMME by analyzing the genetic variants identified in tissues from each of the geographically separated regions. A branching evolution model was proposed, as two of the tumor clades exhibited subclonal differentiation. Genetic diversity was greater in the T4 region than at the other sites; it seemed that T4 comprised two distinct subclones (T4a and T4b), which belonged to clade1 and clade2, respectively. By comparison of VAFs of T4 mutations and SCNAs with other primary tumor samples, it was very obvious that branch-specific mutations or SCNAs showed lower VAF value in T4 than others (Supplementary Fig. S3A and S3B). These results indicated that the subclones of clade1 and clade2 converged in T4 but that competitive subclones with selective advantages had not yet adapted in that region. Recently, a common parallel evolution pattern that involved chromosome 4 instability was suggested by Jamal-Hanjani and colleagues in lung cancer (30). In this study, we also observed chromosome 4 instability, which performed diverse forms in evolution branches (Fig. 4). Ongoing dynamic chromosomal instability was associated with intratumor heterogeneity and resulted in parallel evolution of driver somatic copy-number alterations, which may possibly be associated with poor survival in cancer patients (30).

Early (truncal) mutations likely occur either before or during tumor initiation or during early development, whereas late (branch) mutations reflect mutational processes that shape the genome during tumor maintenance and progression (33). By focusing on putative driver events, we found that critical trunk mutations tended to be clonal, whereas branch mutations tended to be present in subclones (Fig. 4 and Supplementary Fig. S5). Among the four trunk mutations, BRAF and KRAS, particularly BRAF, are notorious cancer-promoting mutations in various cancers, and they account for 25% of melanoma cases in the Chinese population (34). We did not observe a hotspot mutation in BRAF^V600E in this case, but another pathogenic mutation, BRAF^H574Q, was detected, which may also play a significant role in tumorigenesis. KRAS gene mutations are rarely observed in cutaneous melanoma. However, the hotspot mutation KRAS^K12S has been reported in one PMME (17). Notably, trunk mutations in CDH13 and DNAH7 exhibited similar evolutionary patterns as those in BRAF and KRAS, indicating that they might also function during tumor initiation. CDH13 encodes a member of the cadherin superfamily, and it is hypermethylated in many types of cancers (35–37). DNAH7 is a component of the inner dynein arm of ciliary axonemes and drives the sliding of microtubules, which is associated with papillary craniopharyngioma (38). Thus, we suggest that CDH13 and DNAH7 might also be critical during initial tumorigenesis and that more attention should be paid to these genes. However, more biological evidence is needed.

From the phylogenetic topology, we were also able to infer the possible evolutionary trajectories of the tumor metastases. Among the four tested lymph nodes, three were found to be infiltrated by metastatic tumor cells. All of the metastatic tumor subclones were nested in clade1, and T1 was their most recent common ancestor. Metastatic subclones containing pathogenic mutations in ARHGAP26 and PIK3CA, as well as Chr13 arm-level deletion experienced subclonal expansion in metastatic environments (Supplementary Fig. S6). The subclone that carried ARHGEF4 and UMOD gene mutations was detected in T1, T4, LN1, and LN2, and it maintained subclonal dominance in LN1 and LN2. We used TCGA data to test whether above PMME metastatic molecular changes are also metastatic features in other melanoma (Supplementary Table S3 and Fig. S7). It is noticeable that chr13 arm-level deletion is significantly occurred in metastasis TCGA samples (Supplementary Fig. S7), which suggested that this chromosome aberration could be a common metastatic marker in melanoma regardless the sample type. The microenvironment in the lymph nodes might have facilitated positive selection for metastatic tumor cells with such mutations, providing them with growth advantages. Spatially separated lymph nodes had similar subclones and belonged to the same monophyletic group, although most VAFs in the LN3 region were very low. Our results highlight the fact that spatially segregated metastatic subclones can follow very similar evolutionary trajectories.

Another interesting finding was that two of the “normal” mucosal samples were related with the primary tumor clone and may be derived in the very early stage of tumorigenesis. ME1 and N3 shared a certain number of SNVs with the primary tumors. A total of 41.3% of the SNVs (19/46) in N3 were shared with the tumor clones, and 36.2% (29/80) of those in ME1 were also present in the primary tumor (Supplementary Fig. S8). Nevertheless, no significant gene mutations were identified in ME1 or N3. Intensive SCNA analysis could further confirm this relationships. Because ME1 (CCF = 11%) and N3 (CCF = 7%) had extremely low tumor cell content, only SCNAs with BAF > 0.535 based on ASCAT algorithm (Supplementary Fig. S9) could be detected. Two chromosome aberrations that were located in chr7q and chr8 (Fig. 3A) were identified in ME1, and chr8 amplification could also be observed in N3. Chr7q and Chr8q amplification were trunk SCNAs and across both clade1 and clade2 samples, which indicated that N3 and ME1 had a common ancestor with the primary tumor. Moreover, chr15 (60M-telomere) LOH was absent in ME1 or N3, which was another ubiquitous SCNA and had a high detection threshold (BAF > 0.535). These results indicated that ME1 and N3 diverged from other tumor samples in the very early stage. It is reported that PMME has a tendency for vertical growth within the esophageal wall and has a tendency to disseminate by the submucosal lymphatic route (15). In this case, ME1, N3, and the primary tumor distributed vertically along the esophageal wall (Supplementary Fig. S1), so it was most probably that early metastatic clade3 (ME1 and N3) was formed by dissemination through the submucosal lymphatic route. Findings of this study illustrated the mechanisms of early intramural metastasis, which might be a common cause for bad survival in cancer patients (30). This finding prompts the clinicians pay more attention to improve the surgical strategy of PMME to decrease the risk of recurrence or metastasis. Technically, it also warned that using “normal” tissues as germline control is risky for missing some important molecular aberrations.

In this PMME, we observed significant shifts in the mutational signatures during different tumor progression stages. The decreased proportion of C>T mutations was accompanied by an increase in C>A mutations over time. In skin melanoma, C>T is the predominant nucleotide base substitution that may be caused by the misrepair of ultraviolet-induced covalent bonds between adjacent pyrimidines (39). C>T mutations may also be caused by another inductive event, such as the spontaneous deamination of 5-methyl-cytosine, which predominantly occurs at NpCpG trinucleotides (40). We observed a relatively high proportion of C>T transitions at CpG sites in the trunk of the phylogenetic tree; however, within the branches, C>A changes became the most common alteration (Table 1 and Fig. 2B). Some studies have suggested that smoking may be a risk factor for C>A substitutions (24, 41). However, this patient was not a smoker, so it is unlikely that tobacco caused damage to the genomic DNA. We presume that specific processes promoted the generation of genetic diversity in this PMME and that the spontaneous deamination mutational process might have played an important role in founding clone formation; subsequently, some unknown mutational processes resulted in C>A alterations and eventually promoted subsequent subclonal expansions. Patterns of mutation spectrum could also confirmed our presumption that metastatic event of ME1 and N3 started at the early stage, because the mutation spectrum of ME1 and N3 were as similar as the founding clone.

As in any case study, the present study has notable limitations. It involves only one patient, and thus further studies are needed to determine whether the principles discovered here apply to other patients. Despite such limitations, this case provides evidence for tumor evolution models and tumor cell disseminate route at the molecular level, which to our knowledge has not yet been demonstrated in patients with PMME.

In conclusion, by sequencing multiregional samples obtained from a PMME patient, we determined the evolutionary history of this tumor. The presence of ITH and regionally separated driver events, in addition to the relentless and dynamic nature of genomic instability processes, were observed in this case. These findings could help us better understand the processes of PMME tumorigenesis, maintenance, and metastasis. Although large-scale sequencing efforts must be continued for mucosal melanoma and PMME to comprehensively characterize their genomic landscapes, our results have highlighted the importance of viewing tumor development as a dynamic evolutionary process. We have revealed valuable genetic information that contributes to the understanding of the biological mechanism of PMME and that also provides insights that may facilitate the development of biomarkers and therapeutic strategies.

No potential conflicts of interest were disclosed.

Conception and design: J. Li, H. Cai, N. Wu, Y. Ke

Development of methodology: J. Li, Y. Zhou

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S. Yan, Z. Liu, Y. Pan, W. Yuan, N. Wu

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Li, Z. Liu, Y. Zhou, W. Yuan, G. Tian, B. Dong

Writing, review, and/or revision of the manuscript: J. Li, S. Yan, Y. Zhou, N. Wu

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S. Yan, Z. Liu, Y. Pan, M. Liu, Q. Tan, N. Wu

Study supervision: Y. Ke

The authors thank the patient and her family to participate this study.This work was supported by the UMHS-PUHSC Joint Institute for Translational and Clinical Research (grant number BMU20140483 to Y. Ke), Beijing Municipal Science and Technology Commission (grant number Z141100002114046 to Y. Ke), and the Charity Project of National Ministry of Health (grant number 201202014 to Y. Ke).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.