Clonality and Evolutionary History of Rhabdomyosarcoma

Summary

To decipher the dynamic mutational process and identify the causative genomic events in rhabdomyosarcoma (RMS), we developed a systematic method that incorporates multiple types of genomic information to estimate normal cell contamination, tumor clonality, and a timeline of somatic events that occurred prior to the tumor presentation. Our results demonstrate two distinct evolutionary paths resulting in PAX-fusion-negative-rhabdomyosarcoma (PFN-RMS) and PAX-fusion-positive-rhabdomyosarcoma (PFP-RMS): (1) In PFN-RMS, genomic loss of heterozygosity on chromosome 11p15.5 and non-synonymous mutations in RAS pathway and cell cycle genes (FGFR4, KRAS, NRAS, HRAS and CCDN1), as well as several genes not previously known to be drivers of RMS, including PKN1, CUL2, and TTK, occurs early in the evolutionary history of tumor; (2) In contrast, the PAX gene fusion event in PFP-RMS tumors is an early detectable event which consistently occurs prior to a whole genome duplication event. These findings provide new insights into the biology and molecular events that initiate RMS tumorigenesis and may help identify actionable drivers for targeted therapies.

Introduction

Cancer development is driven by dynamic mutational processes and selective pressures which allow a tumor to adapt over time from the initiating oncogenic lesion towards clinical presentation [1–4]. High coverage, next-generation sequencing technologies have provided an unprecedented view of the mutational landscape of whole cancer genomes and demonstrated that cancer genomes have typically acquired thousands of somatic alterations by the time they are clinically detected [5–8]. Although a majority of these alterations do not have clear biological consequences, some alterations are recurrently found; implicating them as critical events in that tumor’s evolution. Subsequently, external factors such as the tumor microenvironment and therapy can confer selective advantages that allow successful clones to eventually supersede one another [1,2]. Two of the major findings of next generation sequencing studies are that biopsies taken at the time of tumor presentation already contain a significant amount of genetic heterogeneity [5,6] and importantly that often a rare subclone in the primary tumor subsequently becomes the founding clone of a metastatic or relapse tumor [5,6,9,10]. The remarkable accuracy and read coverage depth of whole genome sequencing technology now enables the inference of intra-tumor heterogeneity and cancer evolutionary history that is encrypted in its mutational profile [2,9,11–14]. The study of the evolutionary history of a cancer provides insight into which mutations are cancer-driving from the numerous passenger mutations and sheds light on the mechanism of tumorigenesis.

Rhabdomyosarcoma (RMS) is the most common soft-tissue sarcoma of childhood. Despite a growing understanding of the molecular mechanism underlying RMS, the disease continues to have significant mortality and morbidity especially when the tumor is metastatic or recurrent [15–17]. RMS tumors can be subdivided into two major subtypes: PAX fusion positive (PFP) and fusion negative (PFN), characterized by the presence or absence of a oncogenic fusion between the PAX3 or PAX7 and FOXO1 genes [18]. Fusion positive tumors tend to occur in adolescence and are associated with an adverse outcome. Fusion negative tumors typically occur at a younger age and have been associated with significant aneuploidy, loss of heterozygosity (LOH) at chromosome 11p15.5 [19] and mutations of NRAS, KRAS, HRAS [20], PIK3CA, CTNNB1 [21,22] and FGFR4 [22]. We have recently reported that PFP- and PFN-RMS have distinct landscapes of somatic genomic alterations that activate a common molecular pathway [23]. Another effort, using high coverage targeted re-sequencing of two relapsed fusion negative tumors, showed that indeed RMS relapse tumors are derived from a minor subclone discovered in the primary tumor [24].

In this study, we decipher for the first time, the evolutionary history of RMS using high coverage whole-genome sequencing of 44 primary tumors with their matched normal samples. We applied a framework of algorithms that enabled the prediction of the sequential order of mutational events and subclonality using chronological molecular information encoded in somatic mutations and allelic copy number (including copy number changes and allelic imbalance). The major inputs to our method include the variant allele fraction (VAF) of somatic mutations, the variant allele fraction of germline single nucleotide variants in tumor samples and somatic copy number status. We verified the accuracy of mutation detection and VAF estimation by deep sequencing on independent platforms. Our framework first estimates the rate of normal cell contamination and corrects its effect on VAF and copy number. Tumor subclones with different genomic profiles are then identified using the distribution of VAF of somatic mutations as well as the allelic copy number status. Finally, we estimate the temporal progression of the observed somatic events in each tumor genome based on the fact that somatic mutations that occur before and after aneuploidy events will have different values of VAF.

Through our analysis, we discovered that the initiating common lesions of PFN-RMS are the combination of loss of heterozygosity of chromosome 11p15.5 and point mutations in members of the RAS pathway in majority of cases. In a small number of PFN-RMS tumors where no obvious RAS pathway mutations were present, we identified alternative genes that were mutated early in the tumors progression, indicating their potential roles in oncogenesis. In PFP-RMS, we discover that a whole genome duplication event which results in tetraploidy consistently occurs in the middle or late stage of the development of these tumors and that the PAX3-FOXO1 fusion event occurs prior to the whole genome duplication event, making it probable that the fusion is an early event in the evolutionary history of these fusion-positive RMS. Finally, we find that in general RMS tumors are universally composed of a dominant clone, although each tumor contains subclonal populations with a unique mutational profile, which may provide selective advantage to a relapse or metastatic tumor. Because we sequenced only one tumor sample for each patient, a limitation of our study, despite a moderately high coverage of the whole genome sequencing (average 105X), is that we may miss some subclonal mutations and clones especially for subclonal mutations with a VAF<0.1. Nevertheless our findings allow us to propose a developmental model of how this devastating pediatric cancer initiates and evolves prior to presentation.

Results

Large-scale deep WGS analysis of RMS

We used data from whole-genome sequencing (WGS) of 44 primary RMS tumors (19 PAX- fusion-positive, and 25 PFN tumors) with paired blood samples (providing germline status) with an average of 105x coverage per genome base. The accuracy of somatic mutation detection has been estimated as 93% by experiments in which 604 non-silent mutations were verified by whole-exome sequencing and targeted sequencing on independent platforms [23]. To verify the observed VAF of somatic mutations, we re-sequenced the somatic mutations discovered in two RMS samples using multiplex PCR and deep sequencing (1997x coverage) with an orthologous platform. The verification rate of VAF was high by the targeted sequencing (90% accuracy respectively, S1 Fig.). In order to accurately estimate the timing of somatic alterations and dissect intra-tumor subclonality, we first estimated the portion of normal cells contaminating the tumor sample by surveying allelic copy number status and VAF of somatic mutations across the genome (S1 Text). The normal cell contamination rate was generally low in the 44 tumor samples, with a range of 0–33%, a median of 16% and a standard deviation of 9% (S1 Table). The effect of normal cell contamination was numerically corrected for in all subsequent analyses so that the VAF and copy number status are purely for tumor population. The efficacy of normal cell contamination correction is illustrated in S2–S4 Figs., where the observed non-integer allelic copy number was corrected to integers and the observed VAF distribution of somatic mutations was corrected to the expected distribution, e.g. VAF = 0.5 for heterozygous mutations on chromosomes without aneuploidy.

Subclonality analysis

Multiple studies have shown that subclones exist in individuals with cancer which may become a major clone at relapse or progression and this phenomenon has been reported in two fusion negative RMS tumors [24]. To identify those events which may provide resistance to therapy and allow recurrence, we searched for genomic alterations that specifically were confined to the subclonal population. Subclonal copy number alteration events and mutational events were identified using an in-silico method (for which the workflow is shown in S5 Fig.). Fifty point mutations were identified where the variant’s allele frequency was greater than 10% from the expected full clonal allele frequency (S2 Table). Included in this gene list were 7 COSMIC genes including ABL1, BUB1B, CDK12, ERBB2, IGF2, KDR, and SMARCA4. Gene ontology showed enrichment for genes involved in cell adhesion (GO:0007155 p = 0.002), negative regulation DNA repair (GO:0006281, p = 0.05) and cell migration (GO:0030335 p = 0.026). Across the population, several chromosomal events were recurrently found in the subclonal population, including gain of chromosomes, 5, 8, 11, 13, 14, 18 and 19 (S1 Table).

We illustrate our method in more detail using one PAX-fusion negative RMS tumor (RMS2110). Genome wide, this tumor had 3,889 somatic mutations including oncogenic mutation KRAS G13D and FGFR4 V550L [23]. Using the procedure described in the previous section, we estimated that this sample had 12% normal cell contamination (S1 Table). Copy number analysis showed that RMS2110 had large-segment allelic imbalance (including copy neutral LOH) and copy number alteration on eleven chromosomes (Fig. 1A).We employed the combination of somatic mutations, allelic imbalance and copy number alterations to dissect tumor clones that have unique genomic profiles. First, VAF of somatic mutations was determined to find subclonal mutations, based on the fact that a subclonal mutation (present only in a part of the tumor cells) usually has lower VAF than full-clonal mutations (present in all tumor cells) [2]. A typical scenario is on chromosomes without aneuploidy or allelic imbalance, where heterozygous full-clonal mutations have VAF equal to 0.5 while subclonal mutations have VAF<0.5. Analysis of RMS2110 showed that more than half of the chromosomes are without aneuploidy or allelic imbalance (Fig. 1A and 1B). On these chromosomes, a small portion of somatic mutations have observed VAF significantly lower than others indicating the presence of subclones (Fig. 1C, where there are two distinct VAF clusters, one is centered on VAF = 0.5 and the other is centered on VAF = 0.2; for individual chromosomes see Fig. 1D). Using a clustering algorithm with cluster-number-selection procedure (S1 Text) we can identify subclonal mutations as those with lower VAF. Given the depth of sequencing coverage, we estimated that we could detect subclonal mutations with VAF as low as 0.1 (S6 Fig.). Second, allelic copy number status is used to detect subclonal copy number alterations, based on the fact that subclones of different copy number will result in a non-integer allelic copy number for the whole tumor sample. The joint status of total copy number and lesser allele fraction (LAF—the ratio between the less allelic copy number and total copy number, estimated by germline single nucleotide variants, see Methods) reflects whether the allelic copy number is an integer—whether the observed allelic copy number (red dots) is on the expected position (blue crosses) in Fig. 1D. Therefore our approach predicts the subclonal copy number alterations and the fraction of tumor cells possessing these changes (S1 Text).

We thus predict that the subclonal somatic mutations made up only 5% of all the somatic mutations in RMS2110 tumor. Subclonal copy number alterations were detected on chromosomes 5, 7, 8, 12, 13, 17 and 20.

The method was then applied to all the 44 RMS tumors. Looking across all 44 genomes revealed that a dominant clonal lineage was present in each tumor sample (Fig. 2A). The dominant clone carried a large proportion (from 81–96%) of somatic mutations regardless of PAX-fusion status of the tumor. Meanwhile, a small percentage of the cells in each tumor did display evidence of subclonal changes (Fig. 2B). In addition, we observed more subclonal aneuploidy events in PFN tumors (more than half of the tumors have detectable subclones) than PFP tumors (Fig. 2C). The inferred subclonal aneuploidy and mutational events were listed in S1 and S2 Tables, respectively.

Inference of the RMS cancer evolutionary history

The timing of genomic alteration events can be derived using chronological molecular information encoded in the somatic mutations and copy number alterations, because a specific order of genomic event will result in different VAF [2]. To explain the workflow of our method, we again take sample RMS2110 as an example.

The majority of this cancer genome was the expected diploid status with LAF = 0.5, however regions with aneuploidy or allelic imbalance provided the opportunity to identify the timing of genomic events, by comparing the allelic copy number status with the VAF distribution of somatic mutations (S7 Fig.). For example, chromosome 9p and 11 have 2 copies with LAF of 0 (Fig. 1B), which is likely the loss of one allele followed by a duplication of the remaining allele although we cannot formally exclude the possibility of the duplication of both alleles with subsequent loss of the 2 copies of one allele. This observation was confirmed by the fact that most germline single nucleotide variants on chromosome 9p and 11 had VAF near 1 (Fig. 3A). In this case, the somatic mutations occurring before the “LOH+duplication” event must be present on both copies, with an expected VAF of 1, whereas those occurring after the “LOH+duplication” event would be present on only one copy, with an expected VAF of 0.5. The data confirmed this prediction with the VAF displaying a bi-modal distribution with two peaks at 0.5 and 1, respectively (Fig. 3B). The ratio between the numbers of mutations in the two clusters reflects the fraction of “molecular time” it undergoes to accumulate mutations before and after the “LOH+duplication” event, assuming a constant accumulation rate [2]. We acknowledge that the somatic mutation accumulation rate varies among small genomic segments [5,25–27], but for chromosome-level segments used in this study, the average accumulation rates of somatic mutation were observed to be consistent with one another (r² > 0.98, S8 Fig.). Therefore, the molecular timing inferred for different aneuploidy events were comparable among the segments within the same tumor sample.

This analysis was extended genome-wide and allowed the inference of the "phylogenetic tree" of tumor RMS2110 (Fig. 3C). Somatic mutations occurring before and after chromosome aneuploidy events had distinct VAF values, which allowed inference of the timing of two candidate driver mutations in this tumor, the nonsynonymous mutation of FGFR4 V550L and the nonsynonymous mutation of KRAS codon G12. The mutation of FGFR4 V550L has a VAF = 0.95, indicating that the mutation happened before the uniparental disomy event of chromosome 5, which is estimated to occur at 26% of molecular cancer lifetime (S1 Table). Therefore FGFR4 mutation should happen between 0 to 26% of molecular time. Similarly, the mutation of KRAS codon G12 has a VAF of 1, indicating that it occurred before the copy gain of chromosome 12p, timed at 21% of molecular cancer lifetime. Other early events included trisomy of chromosome 8 which happened at 8% of the molecular cancer lifetime and uniparental disomy of chromosomes 11 and 9p which occurred at 14% of the molecular cancer lifetime. Note that on both copies of chromosome 9p, there is a somatic focal deletion of CDKN2A (chr9:21981721-21952010), which is also verified by RNAseq (FPKM = 0).

Other aneuploidy events in this sample included: trisomy of chromosomes 12 and 19 which occurred at 21% of molecular cancer lifetime; uniparental disomy of chromosomes 5, 7 and 17 which happened at 26% of molecular cancer lifetime; tetrasomy of a segment of chromosome 3q which occurred around 50% of molecular cancer lifetime; trisomy of chromosome 13 which occurred at 77% of molecular cancer lifetime; hexasomy of 12p which happened at 85% of molecular cancer lifetime; and trisomy of chromosome 20 which occurred near the time of tumor presentation (Fig. 3C).

This method was applied to all the 44 tumor-normal sample pairs to build the evolutionary history of RMS (Fig. 4, S9 Fig., and S1 Table). In summary, our results showed three major findings. First, LOH of 11p15.5 was a consistent early founding event (in average occurred at 35% of molecular cancer lifetime) in PFN-RMS. In comparison, other highly recurrent aneuploidy events such as the copy gain of chromosome 8 and 2 were not consistent early occurring events (time ranges from 1%∼95% and 16%∼96%, respectively). Second, mutations in RAS pathway genes, including FGFR4, KRAS, NRAS and HRAS, were recurrent early events in PFN-RMS. In addition, mutations in other genes (PKN1, CCND1, CUL2, and TTK) occurred early suggesting their role in tumorigenesis. Third, PFP-RMS tumors in general had much fewer somatic alterations and few of them occurred early in the tumor’s molecular lifetime. Of note, a whole-genome duplication event consistently occurred at the middle or late point in the molecular lifetime of these tumors. The high recurrence suggests that this event might be crucial for the presentation of this cancer subtype. We will discuss these findings in more details in the following sections.

Loss of Heterozygosity of 11p15.5 is a common founding event in PAX fusion negative RMS

The LOH of 11p15.5 is a critical event in fusion-negative tumors [19,28] and frequently results from uniparental disomy and trisomy. In our study, the LOH of 11p15.5 occurred in a total of 26 RMS tumors. As previously described, there was a distinct enrichment for the event in the PFN-RMS population (24/25 PFN samples). Of note, the remaining PFN sample had a small deletion event (3 bps) in the 3-prime non-coding region of IGF2 of undetermined significance. In contrast, only 2 out of 19 PFP-RMS samples had 11p15.5 LOH. The observed LOH usually was accompanied by a copy gain (23 out of 24 PFN samples with LOH on 11p15.5) of the remaining allele resulting in uniparental disomy (n = 8 in PFN and n = 1 in PFP), trisomy (n = 4 in PFN), tetrasomy (n = 8 in PFN and n = 1 in PFP), or pentasomy (n = 3 PFN) (S10 Fig.).

Applying our evolution-history-inference method across the RMS population revealed that the 11p15.5 LOH event universally occurred at an early time point in these tumors development (S11 Fig.). In these tumors, 11p15.5 had a LAF of 0 (S12 Fig. (A)-(C)), meaning that most germline single nucleotide variants present in the tumor, were homozygous. In contrast, the somatic mutations can be homozygous (before 11p15.5 copy gain) or not (after 11p15.5 copy gain). For instance, on chromosome 11 (uniparental disomy) of sample RMS2110, the majority of somatic mutations distributed around a VAF of 0.5 with a minority distributed around VAF of 1 (Fig. 3B). The cluster of somatic mutations discovered around a VAF of 0.5 indicated that these mutations occurred after chromosomal duplication. Comparison of the relative size of the mutations at VAF of 0.5 to those at VAF of 1 indicated that the duplication of chromosome 11 happened early in the cancer evolutionary history. On average, the “LOH+duplication” event of 11p15.5 occurred at 35% of molecular time of the 23 PFN-RMS tumors that had this event. The early occurrence of 11p15.5 “LOH+duplication” event was especially prominent in the patients diagnosed at a later age—the mean occurrence molecular time is 29% of lifetime and the variance is 10% in ≥3 year old patients (S1 Table).

In contrast to the high recurrence (92%) of 11p15.5 LOH in PFN tumors, 11p15.5 LOH was found in only 11% of our PFP tumors. Moreover, in contrast to the consistent early occurrence in PFN tumors, the LOH event did not consistently happen at early lifetime of PFP tumors—37% and 64% of lifetime for two PFP tumors, respectively.

Beyond 11p15.5, the LOH of chromosome 17 (6 out of 25 PFN-RMS) was identified as a recurrent alteration event, typically occurring at an early time point in the evolutionary history of RMS. The LOH region on chromosome 17 includes two small genomic regions that encompass the TP53 and NF1 locus, respectively (S10 Fig.). The observed lesions occurred at an average molecular timing of 38% of cancer evolutionary history. The LOH of 9p, was present in 3 out of 25 PFN-RMS and was found to occur at an average 41% of cancer evolutionary history. Additional accumulation of aneuploidy, such as gain of chromosome 2, 8 and 13, frequently followed the founding events and had no consistent timing pattern.

Whole genome duplication in PAX fusion positive tumors

Unlike PFN-RMS tumors which have significant chromosomal rearrangements, PFP tumors typically have an LAF at the expected 0.5 and a consistent copy number status across the genome. The VAF distribution of somatic mutations revealed that many samples have genome-wide tetraploidy as evidenced by a cluster of somatic mutations with VAF around 0.25 (Fig. 2B, PFP tumors). In total, we found similar VAF distributions in 17 out of 19 PFP tumors. While it is possible that there were two distinct subclones each occupying 50% of tumor cells, given that the cluster with VAF around 0.25 is consistently found across multiple tumors, tetraploidy is a much more likely assumption. Interestingly, when we analyzed these 17 PFP cancer genomes chromosome by chromosome, we found that the inferred tetraploidy occurred in one apparent event (S13 Fig.). This event typically occurred around the mid-point of the molecular lifetime of the tumor (62%±16%) (S1 Table).

In order to identify the sequential order of the PAX3-FOXO1 fusion event and the whole genome duplication event we interrogated the number of the sequencing reads across the PAX3-FOXO1 junction and found that the junction is duplicated. In contrary to PAX7-FOXO1, PAX3-FOXO1 is not known to be focally amplified [29], thus the observed duplication of PAX3-FOXO1 is likely due to the whole genome duplication. Our analysis demonstrated that the PAX3-FOXO1 fusion was consistently duplicated in the whole genome duplication and thus was likely to have occurred prior to the whole genome duplication event (S4B Fig., S3 Table).

Timing the evolutionary history of 44 RMS tumors demonstrates that RAS pathway mutations are critical events in the development of fusion negative RMS

In an effort to obtain a comprehensive picture of the evolutionary history of the RMS tumor genomes, we applied our method to estimate the timing of all the aneuploidy events and mutational events that occurred prior to each tumor’s presentation (Fig. 4). In addition to the recurrent and early 11p LOH event in fusion negative tumors, mutation of FGFR4, KRAS, NRAS and HRAS frequently occurred at an early time point in the cancers evolutionary history (S1 Table). Interestingly, many of the early-mutated genes belong to the RAS pathway (Fig. 4). In total 13/25 fusion negative samples carried the combination of early loss of LOH of 11p15.5 and the mutation of RAS, NF1 or FGFR4 which form the founding events in the evolutionary history of these tumors.

Timing and expression analysis of all the observed somatic mutations was used to discover potential driver mutations that occurred early during tumor development (S4 Table). These mutations included recurrent alteration in two PFN RMS of PKN1 (<30% and <42% of molecular lifetime) which encodes a kinase belonging to the protein kinase C superfamily, CCDN1 at <33% of molecular lifetime in RMS2032, CUL2 at <33% in RMS2034, PAX2 at <77% in RMS 2117, TTK at 17–78% in RMS202.

The E216K mutation of PKN1 inhibits terminal differentiation of C2C12 cells

Recurrent and early mutation of PKN1 in two of the evaluated tumors led us to hypothesize that these changes were involved in myogenic differentiation. To test our hypothesis and assess the potential of our method in identifying driver mutations, we conducted a cell line study described in this section. PKN1 is a member of the AGC-subfamily of serine/threonine kinases. The protein product of the PKN1 gene is composed of a C-terminal kinase domain with significant homology to that of the protein kinase C isoforms, but a unique auto-inhibitory N-terminus made up of 3 homologous stretches of anti-parallel coiled-coil folds (ACC1–3) known to bind to Rho-family GTPases in a nucleotide-dependent manner followed by a C2-like region, known to bind phospholipids and fatty acids. Rho-GTPase binding to the ACC regions causes a conformational change that allows PKN1 to be phosphorylated and activated by PDK1. Active PKN1 plays a role in diverse cellular processes such as regulation of the actin cytoskeleton, cell adhesion, vesicular transport and glucose metabolism [30]. PKN1 represses WNT/CTNNB1 signaling [31] and stimulates the ATF2 and MEF2A transcription factors via a signaling pathway that involves MAP2K3/MAP2K6 and MAPK12 [32]. The observed mutations in RMS202 (E216K) and RMS2035 (A298T) occurred prior to 30% and 42% of molecular lifetime of the tumors, respectively. Both mutations occur in the region of the third ACC domain and could potentially interfere with regulation of the kinase activity of PKN1. To test the functional consequences of the PKN1 E216K mutation, wild type and mutant PKN1 viral constructs were made and transduced into the mouse skeletal muscle precursor cell line C2C12. Using a differentiation assay described by [33], defects in terminal differentiation, reflected in expression of myosin heavy chain (MHC) were observed with the mutated version of PKN1, (Fig. 5A-C), suggesting that the mutant PKN1 prevented the C2C12 cell differentiation. C2C12 cells expressing wild type PKN1 could be induced to express MHC although cell fusion, as determined by the number of nuclei per MHC positive cell, was significantly inhibited. Since PKN1 is known to regulate the activity of several transcription factors known to play a role in myogenic differentiation, we performed expression analysis of the constructed cell lines to determined differentially expressed genes among the constructs. Interestingly, gene set enrichment analysis showed that YAP1 target genes were induced in myoblasts and skeletal muscle genes were repressed in myotubes (Fig. 5D) when the PKN1 mutation was present.

Discussion

While a growing number of studies have used next generation sequencing to examine the evolutionary development of adult tumors [2,9,13], this study is the first application of whole-genome sequencing to the study of a pediatric solid tumor’s subclonality and evolutionary history. The implication from adult studies is that progression of a tumor cell follows a long course, with subsequent lineages of cells acquiring mutations until a particular alteration allows expansion of one clone; ultimately leading to clinical detection. Studies of multiple cancer types including melanoma [34] and colon [35] have shown that even potent oncogenic mutations may be relatively quiescent and thus remain in an undetected premalignant state until released by an additional genetic hit. The patients in our study had a median age at diagnosis of 6.5 years, indicating that the tumors in these children were developed over a relatively short time course, perhaps hinting at the relative potency of the observed somatic changes. Despite a shorter window of development, our study shows that the basic pattern of sequential mutation accumulation is maintained.

Given that similar analyses of adult tumors find the presence of a dominant clone [2], we anticipated that RMS tumors might have a similar finding. Our results show that a dominant clone accounts for more than 80% of cancer cells in a typical tumor sample and genome wide carries many hundreds or even thousands of accumulated point mutations. Even with our genome-wide sequencing and deep coverage, the subclones are frequently observed only at the limits of detection. Admittedly, our study may miss part of the subclonal mutations due to the limit of coverage and the complex nature of genomic subclonality, especially for subclonal mutations with a VAF<0.1. Nevertheless GO analysis of subclonal mutations showed enrichment for potential resistance and progression pathways including cell adhesion, DNA repair and cell migration. While our DNA extraction included relatively large tumor sections, it is also possible that further sampling of the tumor may yield additional subclonal populations as have been appreciated in the study of distant metastases [10]. Even with these constraints, we can conclude that one clone typically dominates these samples, and the majority of tumor cells share most of the detectable genomic alterations. In addition, our findings indicate that a significant proportion of the clonal heterogeneity in RMS is found in gain or loss of copy number, mirroring similar findings in a study of pediatric acute lymphoblastic leukemia [36]. Future efforts will be needed to determine if these alterations might play a role in defining a "cancer stem cell" which resists the selective pressures of the microenvironment and therapy to survive as a recurrence or metastatic lesion.

From the analyses described here, we can begin to understand the dynamics of RMS development. LOH of 11p15.5 has long been described as recurrent feature of several pediatric conditions including the overgrowth phenotype of Beckwith-Wiedemann Syndrome [37], Wilms tumor [38] and embryonal RMS. The proposed mechanism for oncogenesis of this lesion is loss of imprinting control over the IGF2 locus resulting in over-expression of this developmentally regulated growth factor. In this study, we find that not only is this lesion highly recurrent (>90%), it also appears to be the key early landmark in the evolution of fusion-negative tumors. The discovery of a somatic mutation of IGF2 within a fusion negative sample that does not harbor LOH of 11p15.5 (RMS2037) provides additional support to the role of dysregulation of IGF2 in PFN-RMS. In evolutionary terms, the presence of 11p15.5 LOH defines the “most recent common ancestor” when combined with a mutation in a gene in the RAS pathway (NRAS, KRAS, HRAS, FGFR4). While the progression we describe in this study indicates a possible common sequence of events (S1 Table), in some tumors it is equally as likely that the oncogenic mutation of a RAS pathway gene is the founding lesion. This corroborates the finding that patient’s with Costello Syndrome (HRAS germline mutation) [39], Noonan Syndrome (NRAS, KRAS, PTPN11 germline mutations) [40] and Neurofibromatosis (NF1 germline mutation) [41] all have increased risk of developing fusion-negative RMS. Interestingly, while RMS is certainly described in patients with Beckwith-Wiedemann Syndrome (germline uniparental disomy of 11p15.5), these patients appear to have a higher relative risk of developing Wilms’ tumor and hepatoblastoma than RMS [42]. Regardless of which lesion comes first, the combination of LOH of 11p15.5 with a RAS pathway mutation appears to set a clone on the course towards developing a fusion negative RMS tumor (Fig. 6A).

By examining other mutations that occur early in the molecular history of fusion negative tumors, especially those tumors that had no mutation in a candidate “RAS” pathway member, it is possible to nominate other potential founding lesions. In this analysis, mutation of PKN1 stands out as both recurrent (occurring in tumor RMS202 and RMS2035) and occurring early in the evolutionary history of these tumors (<30% and <42% respectively). PKN1 is a member of the protein kinase C family and has been implicated as a repressor of WNT/CTNNB1 signaling [43], a mediator of insulin signaling to the actin cytoskeleton [44] and an activator of MEF2A dependent transcription [32]. The observed mutations (E216K and A298T) occurred within the third ACC domain and may relieve N-terminal auto-inhibition of PKN1 kinase activity (S14 Fig.). Our functional assessments of the PKN1 mutation demonstrated that the observed mutation can inhibit terminal differentiation of skeletal muscle in a dominant-negative fashion. Gene expression analysis of PKN1 E216K expressing C2C12 cells grown in low serum media reveals that this differentiation block occurs as a result of repression of skeletal muscle specific genes and genes with MEF2A binding sites in the promoter. Expression of MEF2 proteins is induced by myogenin during normal skeletal muscle differentiation, and MEF2 factors, in complex with other transcriptional activators, are known to play a role of in myogenic differentiation [45]. One possible explanation for these results is that PKN1-dependent signaling leads to replacement of MEF2A at skeletal muscle specific genes with a transcriptional repressor. The identity of this transcriptional repressor is not yet known.

GSEA analysis of PKN1 E216K expressing C2C12 myoblasts reveals that the gene expression signature of these cells is enriched for a YAP gene signature. The YAP gene signature is composed of genes activated by YAP overexpression in human mammary cells, mouse liver tissues and immortalized mouse fibroblasts, and thus represents a list of evolutionary YAP target genes [46]. YAP1 has recently been identified as a potent ERMS oncogenic driver, and YAP1-TEAD complexes repress expression of genes important for myogenic differentiation, in part by impinging upon the binding of myogenic transcription factors including MEF2 [47]. Intriguingly, YAP1 is able to functionally compensate for oncogenic KRAS in colorectal and pancreatic cancer mouse models [48,49], which is consistent with the fact that RAS mutation and PKN1 mutation are mutually exclusive in the set of ERMS tumors presented here. The mechanism by which PKN1 leads to increased YAP1 activity is currently unknown.

With whole genome sequencing from a single biopsy, we were unable to reliably estimate the molecular timing of some of the recurrent genetic lesions observed in RMS. These events include mutations in PIK3CA, FBXW7, NF1, and BCOR, and the focal amplification of MDM2, CDK4, and MYCN. Moreover, somatic events such as the copy gain of chromosomes 2, 8 and 13 were estimated to occur in an inconsistent temporal pattern among RMS patients (occurring at multiple time points throughout the evolutionary history). Certainly, the recurrence of these lesions in multiple large sample sets speaks to their importance, however since they can occur either early or late in the cancers evolutionary history they do not represent the initiating lesion in the tumor but rather modifiers of the tumor as it grows towards presentation.

Early flow cytometry studies to evaluate DNA content [50] and cytogenetic analysis [51] in ARMS noted the frequent presence of tetraploidy. We find evidence that fusion positive tumors demonstrate a very high recurrence of a whole genome duplication event that results in a tetraploid genome. In general, the duplication event occurs around the midpoint of the tumor’s evolutionary history and prior to additional accumulation of aneuploidy (Fig. 6B). The mechanism of a cell moving from a diploid state to a tetraploid state can be due to a cytokinesis failure, a cell fusion event or mitotic slippage and it has long been speculated that tetraploidy is an intermediate in cancer progression [52]. This is born out in at least one premalignant state (Barrett’s esophagus), where tetraploidy represents an early intermediate as the cancer precursor develops towards esophageal adenocarcinoma [53]. Further experimental evidence of a tetraploidy intermediate state is demonstrated in TP53 null mouse mammary cells where tetraploid but not diploid cells promote tumorigenesis in nude mice [54]. Of note, in the mouse model of PFP ARMS both the increase in allelic copy number of the PAX3-FOXO1 fusion and the loss of TP53 or CDKN2A were found to be a critical component in increasing the penetrance of the phenotype [55]. Given the role of loss of TP53 in allowing a permissive environment for the tetraploid cell to escape cell cycle arrest, this is consistent with our findings. It is interesting to speculate that in our tumor series the whole genome duplication event is an attempt by the tumor to increase the allelic dose of the PAX fusion. Another intriguing possibility is that the whole genome duplication event alters the telomere length and activates telomerase; a critical cooperating genetic event in the temporal sequence that produces alveolar rhabdomyosarcoma from human myoblasts [56].

In conclusion, based on deep whole genome sequencing, we developed a systematic method to infer the evolutionary history and identify the causative lesions of this pediatric solid tumor. From our data, we provide a model for how these tumors develop. Our results demonstrate two distinct evolutionary paths resulting in a convergent phenotype of this soft tissue cancer: 1) genomic loss of heterozygosity on 11p15.5, nonsynonymous mutations on RAS pathway and cell cycle genes, including FGFR4, KRAS, NRAS, HRAS and CCDN1, as well as several other genes, including CUL2, TTK and PKN1, occur early in the evolution history of PFN-RMS; 2) recurrent whole genome duplication occurs in the middle or late stage of the tumors evolution in the PAX-fusion positive RMS tumors. PAX3-FOXO1 fusion occurs before the whole genome duplication event. Intriguingly, a recent report of the clonal evolution in two cases of relapsed fusion-negative RMS tumors after chemotherapy and radiation demonstrated that recurrent tumors are derived from the founding tumor’s minor clone [24]. Nevertheless given that survival following relapse is 30% [57], our findings that LOH of 11p15.5 and mutation of RAS genes form the “trunk” of the fusion negative RMS tumors has important therapeutic implications; at least in theory it is preferable for treatment to target mutations that are present in all of the tumor cells. In addition, while further biologic validation is required, our analysis identifies candidate lesions that may underpin the ability of the minor clone to survive therapy and proliferate during relapse or metastasis. Thus given clonal evolution and heterogeneity, our data suggests that future efforts to understand the emergence of resistant clones should include deep sequencing analysis on all patients at presentation and relapse. These efforts will add to the growing understanding of the biology of RMS and identify actionable genetic aberrations for targeted therapies.

Methods

Patient samples

All patient sample collection was approved by the institutional review board of the participating facility. All de-linked and de-identified patient sample information was collected under the approved institutional review board of the National Cancer Institute protocol 10-C-0086. Samples were assembled from collections at the Pediatric Oncology Branch of the National Cancer Institute, Children’s Oncology Group, and the Tumor Bank at The Children’s Hospital at Westmead (New South Wales, Australia). All tumors were collected at initial diagnosis and prior to any therapy with the exception of samples NCI0040 and NCI0080 which were collected at relapse. Samples were de-identified and histologic diagnosis and clinical information were compiled. Quality control genotyping for the whole genome samples was performed to ensure the match of tumor normal pairs.

Whole genome sequencing

Approximately 6 micrograms of DNA was sequenced using the paired end sequencing method of Complete Genomics. Data analysis was accomplished using CGA tools package v2.0 as well as a number of in house tools described in [23].

Targeted resequencing

A custom panel of oligonucleotides designed to incorporate somatic single nucleotide variants discovered in samples RMS2110 and RMS2107 was generated using Ion Ampliseq designer software. 150 base pair amplicon libraries were generated using multiplex PCR according to the Ion Torrent Ampliseq Library 2.0 kit. Individual samples were barcoded and the generated libraries were sequenced using a 318 chip on a Personal Genome Machine (Life Technologies). Compiled reads were mapped to Hg19 and the expected variants were analyzed for coverage and VAF. The targeted sequencing has a depth of 1997x. We used it to verify the VAF estimates of somatic mutations across whole genome of two samples. The result shows VAF estimation accuracy is as high as 90% (S1 Fig.).

SNP array

Illumina Omni 2.5M (97 paired plus 30 unmatched tumors) or 5M (10 paired samples) were performed according to the standard procedure from the manufacturer (Illumina, San Diego, CA) at the National Cancer, Cancer Genomics Research Laboratory. The data were previously reported in [23].

Whole transcriptome sequencing

PolyA selected RNA libraries were prepared for RNA sequencing on Illumina HiSeq2000. 100 bases long paired-end reads were assessed for quality and reads were mapped using CASAVA (Illumina, San Diego, CA). The generated fastq files were analyzed by TopHat2 [58] and Cufflinks [59] The data were previously reported in [23].

Somatic mutation calling

We used the somatic score from Complete Genomics (CG), a highly sensitive and specific somatic mutation-calling criteria [23,60,61], combined with principles previously used [25,62] to detect somatic signal nucleotide variants (SNVs). The somatic score is designed by systematically considering sequencing error, mapping error and read count and has been validated as an effective criterion by our previous studies as well as studies from other groups [23,61,63]. To ensure sensitivity and specificity of mutation calling, we set the somatic score cutoff as ≥0, with a set of additional filters for removing system artifacts and mapping errors, to select somatic mutations, based on the verification results from an independent sequencing platform [23]. By SOLiD whole exome sequencing (35x coverage) and Ion Torrent Ampliseq targeted sequencing (300x coverage) on 30 RMS tumor-normal sample pairs sequenced by CG, we verified that our filtering removes 99% of false positives while maintaining 80% of true positives [23] in CG’s comprehensive somatic mutation pool (called in a loose criterion to be as inclusive as possible) [23,60], which is premium given the less somatic mutations in pediatric tumors [3].

Copy number profiling

For the 44 matched tumor and normal sample pairs with whole genome sequencing, we used somatic copy number segmentation profile provided by Complete Genomics, with customized corrections. For a larger cohort, with 120 matched tumor and normal sample pairs measured by SNP array, we used the copy number profile provided by NEXUS copy number analysis.

In copy number profiling, an important quantity, lesser allele fraction (LAF), is estimated to represent the fraction of copies coming from each parental allele. Along with total copy number, LAF tells us the allelic copy number status. LAF is defined as the ratio between the copy number of the lesser allele (with fewer copies than the other allele) versus the total copy number. For example, a trisomic chromosome with allelic type “AAB” has LAF equal to 1/3. LAF is estimated based on germline heterozygous single nucleotide variants that are present in cancer genome. For each genomic segment of tumor sample, sites that have heterozygous single nucleotide variants in matched normal samples were selected and LAF was estimated by the ratio between the read count of lesser allele and the total read count for these sites.

Statistical analyses

In order to study the evolutionary history of RMS, we performed an integrated statistical analysis, using information including germline single-nucleotide variants, somatic mutations (single-nucleotide variants), somatic copy number alterations (CNA) and junctions (mapped breaking points), to estimate the normal cell contamination, intra-tumor heterogeneity (due to subclones) and the timing of somatic variants. The methods are detailed in the S1 Text.

Cell lines and retroviral constructs

Mouse myoblast cell line C2C12 was obtained as a generous gift from Dr. Marc Landanyi. The cells genotype was performed by the NCI Core Genotyping facility and the cell line was confirmed to be myocoplasm negative. All cell culture was performed in DMEM supplemented with 10% FBS. The differentiation assay was performed using 2% horse serum (Life Technologies) as previously described [33]. For retrovirus production, the pBabe vector system (Addgene) was used.

Retrovirus construction and transduction of C2C12 cells

Plasmids encoding human PKN1 cDNA was purchased from Addgene, and constructs were subcloned into pBABE containing a N-terminal HA tag. A cDNA encoding the PKN1 E216K mutant was generated using the GeneART Site-Directed Mutagenesis kit (Life Technologies) and was subcloned into pBabe. The generated mutation was confirmed by Sanger sequencing. Retroviruses were generated by contransfection of pBabe constructs with pCL-10A1 into 293T cells (American Type Culture Collection (ATCC), CRL-3216) subsequently used to infect the C2C12 cell line as previously described; infected cells were selected with 2 μg/ml puromycin (Life Technologies).

Immunofluorescence

Cells grown on Nunc chamber slides were fixed with 4% paraformaldehyde and permeabilized in PBS containing 0.5% Triton X-100 and blocked in Block-Aid (Life Technologies) for 1 h at room temperature. Cells were then incubated with MF20 monoclonal antibody (DSHB) against MHC (1:40 dilution;) overnight at 4 degrees. Secondary antibody Alexa Fluor 488-conjugated secondary antibody (1:200 dilution; Life Technologies) for 1 h at room temperature. Cells were mounted with ProLong Gold antifade reagent with DAPI (4′,6-diamidino-2-phenylindole; Life Technologies).

Immunoblotting

Cells were lysed in M-PER lysis buffer (Pierce Biotechnology). Lysates were denatured in 4× sample buffer at 70°C for 10 min, resolved on 4–12% NuPAGE gels (Life Technologies) and transferred onto PVDF (polyvinylidene fluoride) membranes. Membranes were blocked in 5% nonfat milk in TBST buffer (TBS with Tween-20) for 1 h at room temperature and probed with primary anti-HA antibody obtained from Covance (1:2,000 dilution). Bound antibodies were detected with peroxidase-labeled horse antibody to mouse IgG and visualized using enhanced chemiluminescence reagents (ThermoScientific).

Expression analysis

Total cellular RNA was isolated using the RNAeasy mini kit (Qiagen). Cellular RNA (250ng) was in vitro transcribed, fragmented, hybridized and applied to Affymetrix Mouse 430A arrays according to the standard operating procedure of the Laboratory of Molecular Technology core facility (http://atp.ncifcrf.gov/genetics-and-genomics/laboratory-of-molecular-technology) and the manufacturer’s instructions (Affymetrix, Santa Clara, CA). For gene set enrichment analysis (GSEA), the normalized gene expression data were z-scored and ranked according to absolute fold-change expression over the control. GSEA analysis (http://www.broadinstitute.org/gsea/index.jsp) was performed using default parameter settings.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15

Attachment 16

Attachment 17

Attachment 18

Attachment 19

Attachment 20