Lineage-defined leiomyosarcoma subtypes emerge years before diagnosis and determine patient survival

Leiomyosarcomas (LMS) are genetically heterogeneous tumors differentiating along smooth muscle lines. Currently, LMS treatment is not informed by molecular subtyping and is associated with highly variable survival. While disease site continues to dictate clinical management, the contribution of genetic factors to LMS subtype, origins, and timing are unknown. Here we analyze 70 genomes and 130 transcriptomes of LMS, including multiple tumor regions and paired metastases. Molecular profiling highlight the very early origins of LMS. We uncover three specific subtypes of LMS that likely develop from distinct lineages of smooth muscle cells. Of these, dedifferentiated LMS with high immune infiltration and tumors primarily of gynecological origin harbor genomic dystrophin deletions and/or loss of dystrophin expression, acquire the highest burden of genomic mutation, and are associated with worse survival. Homologous recombination defects lead to genome-wide mutational signatures, and a corresponding sensitivity to PARP trappers and other DNA damage response inhibitors, suggesting a promising therapeutic strategy for LMS. Finally, by phylogenetic reconstruction, we present evidence that clones seeding lethal metastases arise decades prior to LMS diagnosis.

Leiomyosarcoma (LMS) is a malignant neoplasm of smooth muscle differentiation, which accounts for 10–20% of sarcoma diagnoses1. The heterogeneity of its site of origin and clinical course, including the development of metastasis and response to therapy, makes the treatment of LMS particularly challenging. Common sites of presentation include the extremities, the uterus, and the abdomen, the latter of which predominately arises from the retroperitoneum, but also from the viscera. While all LMS, regardless of anatomical origin, differentiate along the smooth muscle lineage, the disease site continues to be used to stratify care pathways. Rates of metastasis at 10 years vary by disease site (31% in extremity, 58% in the abdomen, 53–71% in the uterus)2,3. However, it is unknown whether the molecular factors leading to these site-specific survival differences are inherent to the lineage, differentiation status or mutational processes operative in the cells of origin from which the disease develops, or are acquired, involving unique driver gene mutations.

As the primary pattern of failure is the development of metastasis, outcomes for LMS patients have not improved as the overall effectiveness of current systemic treatment remains poor4. Recent studies using targeted sequencing and/or RNA-seq have confirmed frequent alterations of TP53, RB1, PTEN, and ATRX, as well as recurrent amplifications of 17p11.2-p12 (MYOCD)5,6. One of these studies also suggested that LMS may be sensitive to poly (ADP-ribose) polymerase inhibition (PARPi) due to acquired mutations in BRCA1/BRCA2. The “BRCAness” mutational signature was reported in 48/49 (98%) patients in this study5. If validated, this would support the use of PARPi in LMS. Whole-genome sequencing (WGS) typically yields ~100X more mutations than exome sequencing and therefore, is the preferred approach for validating and refining mutational signature analyses, especially for separating related signatures from each other7.

The majority of LMS patients present with non-specific symptoms at 50–60 years of age. These tumors may develop gradually over the course of several years. A better understanding of the timing of mutations in LMS would therefore provide knowledge about tumor development, and perhaps aid in early detection. Similarly, understanding the ongoing mutational dynamics of LMS, including its mutation rate and overall model of tumor evolution, could define rational surveillance strategies for patients prior to the development of metastatic spread.

In this work, we study the mutational processes underlying primary and metastatic LMS by analyzing 70 whole genomes and 130 transcriptomes from 113 patients. This includes samples from multiple areas from a single tumor and/or distant metastatic relapses separated in time. We observe widespread genetic diversity within primary tumors and between metastatic relapses, especially with respect to rearrangements and clustered mutations that are only detectable by WGS. The results provide compelling support for molecular subtyping in LMS, and evidence for early systemic spread years prior to diagnosis. This study also highlights the potential use of a DNA damage inhibitors (DDRi), including PARP trappers, as a promising therapeutic avenue in these patients.

Whole-genome and transcriptome sequencing of a unique LMS cohort

To reconstruct the molecular events underpinning LMS, we carried out an analysis of the patterns, location, and evolution of somatic mutations in LMS, both between and within tumors. In total, 70 genomes and 130 transcriptomes of LMS were analyzed (Supplementary Fig. 1A). Of these, 53 samples from 34 patients were newly genome sequenced. An additional 18 whole genomes, previously reported by The Cancer Genome Atlas (TCGA)6, were included. Only validated LMS samples approved by the TCGA-SARC program were used (Supplementary Fig. 1B, see Methods). Matched-blood or tissue was used as a reference for all patients. We conducted a detailed pathological review of all specimens including those from TCGA. During this review, we detected a pathogenic KIT variant known to be associated with gastrointestinal stromal tumors8 (GIST) in one TCGA LMS sample, prompting its removal from further analysis. All genomes were processed through the same established informatics pipeline9 to detect somatic mutations, including substitutions, small insertions or deletions (indels), copy number changes, structural rearrangements as well as clustered and complex events (Supplementary Data 1–6). For seven patients, we sequenced multiple areas from a single tumor and/or distant metastatic relapses separated in space or time. Sixteen informative samples from five of these patients were selected for additional, targeted, deep sequencing (~700X) to validate substitutions, indels, and structural variants present in these multiregion or matched primary-metastatic samples (Supplementary Fig. 2). Full RNA sequencing analysis for expression and fusion gene analysis was also carried out for 130 (51 Toronto, 79 TCGA) transcriptomes. To our knowledge, this work represents the broadest analysis of somatic mutations in LMS to date.

Three molecular subtypes of LMS, with distinct genomes and transcriptional programs, correlate with patient survival

Although labeled as a single disease characterized by smooth muscle differentiation, the clinical presentation and behavior of LMS are highly variable10. The degree to which this variability is explained by overall genomic and transcriptomic features is largely unknown. As a starting point, we examined the tumors’ whole transcriptomes and found three predominant gene expression subtypes by principal component analysis that were not fully accounted for by disease site (Fig. 1a). We observed three LMS subtypes consistent with previous studies5,6,11,12,13,14 that harbored established gene expression markers of survival (Supplementary Fig. 3A). Subtype 1 was mixed: 43.4% gynecological, 30.4% extremity, 13.0% abdominal, 8.6% metastases, and 4.3% other. Subtype 2 largely consisted of abdominal and to a lesser degree extremity tumors (63.5% abdominal, 17.6% extremity, 1.1% gynecological, 15.2% metastases, and 2.3% other), while subtype 3 was composed almost exclusively of gynecological LMS (90.9%, 9.1% other). The transcriptome-defined subtypes were also associated with differences in genome-wide mutation burdens (Fig. 1b, Supplementary Data 7). LMS subtypes 1 and 3 (mixed and gynecological, respectively) harbored a higher overall burden of somatic mutations and displayed inferior overall and disease-specific survival compared to subtype 2 (abdominal/extremity; Supplementary Fig. 3B, C). This was true for all classes of mutations, including substitutions (2.7 vs 1.9 Mut/Mb, n.s.), indels (605 vs 359 v indels, p < 0.05, Welch’s t-test), and large-scale rearrangements (224 vs 82 rearrangements, p < 0.01, Welch’s t-test).

a Principal Component Analysis (PCA) of LMS transcriptomes (n = 79 TCGA RNA-seq and n = 51 Toronto RNA-seq) leads to three defined subtypes of LMS, in which segregation is broadly influenced by anatomical location. While subtype 1 is a mix of all sites (extremity, abdominal and gynecological), subtype 2 is largely abdominal, with some extremity tumors. Subtype 2 can be further sub-stratified into subtype 2a and 2b. The metastatic lesions in subtype 2 cluster with their matched primary tumors. There are two metastatic tumors in subtype 1 that are from the same patient. Lastly, subtype 3 represents a largely gynecological (uterine, vaginal, fallopian tube) subtype. b Genomic point mutation, indel (insertion/deletion), and structural variant (SV) burdens are lower in subtype 2 than subtypes 1 or 3. Horizontal black lines represent the median values for each subtype. Source data are provided in Supplementary Data 7. c LMS molecular subtypes are of distinct smooth muscle lineages: vertices of the triangular plot represent smooth muscle of vasculature, digestive tissue, and gynecological tissue. Individual dots represent LMS cancers and where they lie in the cluster. Adjacent contour plots illustrate density distribution of LMS molecular subtypes. d Uniform Manifold Approximation and Projection (UMAP) illustrates clustering of 271 muscle-related GTEx normal tissue types (from the Genotype-Tissue Expression Program) and 130 LMS reveals distinct smooth muscle lineages of LMS subtypes. e Boxplots represent the expression (in transcripts per million, TPM), for smooth muscle (SM) genes: LMOD1 (leiomodin 1), MYOCD (myocardin), DES (desmin), and CALD1 (caldesmon) are key smooth muscle genes that are commonly expressed in LMS. The boxes represent the 25th and 75th percentile (bottom and top of box), and median value (horizontal band). The whiskers indicate the variability outside the upper and lower quartiles. These genes are highly expressed in vascular (n = 110), digestive (n = 119), and gynecological (Gyn., n = 42) normal smooth muscle. Genes are also expressed in subtype 2 (n = 85) and subtype 3 (n = 22) LMS, but not as highly in LMS subtype 1 (n = 23).

Full size image

We then validated these expression subtypes using an orthogonal approach based on Uniform Manifold Approximation and Projection (UMAP15) and Density-Based Spatial Clustering of Applications with Noise (DBSCAN16; see Methods). To optimize subtype generation and classification, we took two experimental approaches. First, we included 12,419 additional tumors of various histologies in our analysis. Second, we clustered the 130 LMS tumors with 1735 normal tissues from cancer-free donors17,18. We found the same subtypes using our pan-cancer cluster approach, except we were now able to divide subtype 2 into two additional classes: subtype 2a, which is comprised almost exclusively of abdominal LMS, and subtype 2b, which consists of abdominal and extremity tumors. (Supplementary Fig. 4). Out of the 130 RNA samples, three cases (2.3%) clustered apart from the others––falling outside the three LMS subtypes. We reviewed the histopathology from each case but did not find any further evidence to support a change in diagnosis. In one instance, a gastric LMS resembled stomach cancers from TCGA. Differences in the tissue of origin may explain these few exceptional LMS tumors. In our second approach, matching tumors to normal tissues, we found that three main tissue types strongly resembled LMS (out of 51 types): (1) Digestive smooth muscle (esophagus muscularis/gastroesophageal junction, colon, bladder); (2) vascular smooth muscle (tibial/coronary arteries, aorta); and (3) gynecological/fallopian tube smooth muscle (Fig. 1c, d, Supplementary Fig. 5). As expected, gynecological LMS most resembled normal gynaecological tissue. Subtype 2a (almost exclusively abdominal) clustered with digestive smooth muscle, while subtype 1 (mixed) and 2b (abdominal/extremity) mostly resembled vascular smooth muscle, consistent with extremity and retroperitoneal LMS being associated with vasculature. These analyses define the cellular lineage underlying the LMS molecular subtypes.

Genomic alterations in muscle-related genes segregate in LMS subtypes

Having seen how closely the molecular subtypes of LMS adhere to their respective smooth-muscle lineages and knowing certain muscle marker genes are overexpressed in some LMS11,12, we next searched for alterations in genes related to smooth-muscle differentiation and function. We looked for mutations (point mutations, indels, structural variants, or focal copy number changes) in smooth muscle marker genes: LMOD1, CALD1, MYOCD, DMD, ACTG2, DES, CFL2, and SLMAP. The most frequent alterations detected were recurrent intragenic deletions of the dystrophin (DMD) gene in 8/51 (16%) and focal myocardin (MYOCD) amplifications in 20/51 (39%) patients’ genomes (Supplementary Fig. 6A–C). In fact, DMD deletions were as frequent as ATRX deletions and among the most recurrent structural variant changes in LMS, after RB1 (Supplementary Data 3). Intriguingly, DMDdeletions occur predominantly in LMS subtype 1, and to a lesser extent subtype 3, while MYOCD amplifications occur more commonly in subtypes 2 and 3 (Supplementary Fig. 6d). Of note, MYOCD (17p12) amplifications frequently accompanied TP53 (17p13) losses in a joint segmental gain/loss pattern (Supplementary Fig. 6B), and in five cases, MYOCD was amplified via complex copy number patterns (Supplementary Fig. 6C). With respect to DMDalterations, it is well-established that dystrophin deficiency causes Duchenne and Becker muscular dystrophies; however, its role in cancer is not well defined. Recent data have emerged suggesting that dystrophin can act as a tumor suppressor and limits metastasis in myogenic cancers19. Furthermore, it was reported that DMD deletions, that were detected by SNP array, occurred in 3/7 primary LMS and 8/13 metastatic LMS that abrogated expression of full-length dystrophin transcript, which encodes Dp427m19. To investigate this further, we searched for the loss of expression of the Dp427m transcript in samples where a genomic DMD deletion was detected and found either no expression or lower expression of its transcript in 7/8 (88%) cases (Supplementary Fig. 7A). Surprisingly, and irrespective of molecular subtype, the majority of gynecological LMS (21/28, 75%) had a reduced or complete lack of expression of the Dp427m transcript even when a DMD deletion was not detected, suggesting that there may be alternative mechanisms of dystrophin inactivation. This appears to be a unique feature of gynecological LMS that is not seen in soft-tissue LMS without a DMD deletion or normal uterine tissues (Supplementary Fig. 7B). To confirm that the dystrophin isoforms seen in LMS resemble those seen in Duchenne muscular dystrophy (DMD) patients, we RNA sequenced muscle biopsies from six patients with DMD and observed that DMD in dystrophin-deleted LMS transcriptomes do have similar expression to muscular dystrophy patients (Supplementary Fig. 7A). Taken together, we find frequent alterations in the DMD gene and loss of full-length DMD expression associated with LMS subtypes 1 and 3 that are consistent with inferior outcomes.

Subtype 1 LMS are associated with myogenic dedifferentiation and high immune infiltration

Given the prevalence of deletions of dystrophin––a critical muscle protein––in subtype 1, we wondered if subtype 1 constituted an overall less differentiated form of LMS, which is known to be a prognostic marker of poor outcome20. Consistent with this notion, most markers of muscle differentiation, including leimodin1 (LMOD1), caldesmon (CALD1), and MYOCD were diminished in subtype 1 (Fig. 1e, Supplementary Fig. 8). We then compared LMS to undifferentiated pleomorphic sarcoma (UPS), obtained from the TCGA, as these two entities can be difficult to distinguish histologically21. Most of the UPS tumors clustered with LMS subtype 1––the only subtype to have non-LMS cancers grouped with it (Supplementary Fig. 9). This suggests that LMS subtype 1 may represent an aggressive dedifferentiated form of LMS.

For many solid tumors, dedifferentiation towards ‘stemness’ is immunosuppressive22,23, however, it is unclear how this relates to LMS and its molecular subtypes. To determine the relationship between the immune microenvironment and the state of dedifferentiation in LMS, we obtained leukocyte proportion and cell type for the 79 TCGA LMS samples24 and performed an in silico immune analysis. Consistent with previous work25,26,27, M2 macrophages were the most prevalent immune cell in LMS. Interestingly, we found that the dedifferentiated subtype 1 LMS harbored a higher leukocyte content, specifically M2 macrophages, than subtypes 2 or 3 (Supplementary Fig. 10A, B, see Methods). Of note, the established marker for this subtype, ARL4C expression, was overexpressed in subtype 1 LMS, as expected11 (Supplementary Fig. 10C). Tumors with higher leukocyte fractions are known to be the most responsive to immune checkpoint inhibition, as recently demonstrated in UPS and dedifferentiated liposarcoma28. In LMS, the infiltrate composition reflects an M2 macrophage-dominated, low lymphocytic profile that is consistent with an immunosuppressed tumor microenvironment. This may render LMS, and in particular poorly differentiated subtype 1 LMS, more responsive to macrophage-focused immunomodulatory agents29,30. However, these findings warrant further investigation.

Genomic substitution and indel signatures identify prior therapy exposure and a treatment strategy for LMS

Next, we performed a mutational signature analysis to determine if LMS subtypes are driven by different mutagenic processes. We began with an updated de novo extraction approach using an established signature extraction tool, SigProfiler7,31 which differs from previous ‘refitting’ approaches used on LMS exomes5. Using SigProfiler for LMS mutational signatures, we optimized our approach for the increased number and variety of mutations detectable by WGS, including signatures based on single-base substitutions (SBS), double-base substitutions (DBS), as well as insertions and deletions (ID). Owing to the larger number of data points provided by WGS, previously correlated signatures could now be separated. One such signature is the HR-deficiency or ‘BRCAness’ signature (SBS3). We validated our power to discriminate SBS3 from other ‘flat’ signatures by using two approaches. First, we demonstrated a strong correlation (>99% Pearson coefficient) between our output and the ground-truth catalogue from ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Network (Supplementary Fig. 11, see Methods). Second, we removed SBS3 from the reconstruction of the mutational profile and evaluated the changes in cosine similarity (Supplementary Fig. 12, see Methods). Thus, using a de novo extraction approach with WGS data points we can confidently identify specific intrinsic sources of mutation.

27 established signatures found in the COSMIC database were identified in this group of 70 LMS tumor genomes, including fourteen SBS, five DBS, and eight ID signatures (Fig. 2a, Supplementary Fig. 13). Of these, 13 have a known or proposed etiology. Four signatures of sequencing artifacts (SBS45, SBS46, SBS49, SBS52) were removed (see Methods). Two ID signatures and one DBS not found in COSMIC were also identified and are of unknown cause (Supplementary Fig. 14). Of note, ID8 was detected in every sample known to have been treated with radiation therapy (9/9), consistent with the suggested etiology of this signature. This included two cases where patients received prior radiation for previous cancer, suggesting LMS is secondary cancer in these individuals. Overall, mutational signatures highlighted the importance of DNA repair defects in LMS. We found SBS8, linked to a deficiency in nucleotide excision repair32, as well as SBS3 and ID6, which are signatures of defective homologous recombination-based DNA damage repair31,33. Either SBS3 or ID6 was detected in 39/61 (63.9%) of LMS samples (13/65 SBS3 alone, 18/65 ID6 alone, and 8/61 in both). Thus, we are able to confirm the existence of HR-deficiency substitution signatures, albeit at a lower frequency than previously reported5, as well as the existence of a corresponding HR-indel signature in LMS.

a Non-negative matrix factorization (NMF)-extracted and decomposed single-substitution (SBS), indel (ID) and double-nucleotide signatures (DBS) are illustrated in the heatmaps. Common substitution signatures include SBS1, SBS5, SBS8, and SBS40. SBS3 and ID6 (HR-deficiency) are found in 64% of samples. SBS2, SBS13, and DBS11 reflect localized hypermutation events, also called ‘kataegis’. ID8 represents a radiation signature, commonly seen in patients treated with radiation therapy. ‘Other’ substitution signatures, present in less than 5% of samples, can be found in Supplementary Fig. 13. Color refers to signature activity. b Evaluation of sensitivity to DNA damage response pathway, including PARPi, in soft tissue (ST) LMS cell lines (STS39, STS54, STS137, STS210, and STS551) and gynecological LMS cell lines (SKLMS-1, SK-UT-1 and SK-UT-1B). c Representative boxplots of EC50 from LMS (n = 8) and UPS (n = 5) cell lines treated with the PARP inhibitors, talazoparib, and olaparib. The boxes represent the 25th and 75th percentile (bottom and top of box), and median value (horizontal band). The whiskers indicate the variability outside the upper and lower quartiles. For olaparib treatment, boxplots were generated for seven LMS and three UPS cell lines only, as growth suppression failed to occur in the remaining one LMS and two UPS cell lines along with the RPEΔp53 control. In contrast, all LMS cell lines are responsive to talazoparib (median EC50 0.37 µM) compared to UPS cell lines (median EC50 6.26 µM, p = 0.072, one-sided Welch’s t-test). Detailed information for all patient derived cell lines (LMS and UPS) can be found in Supplementary Data 8 and 9. d The Traffic Light Reporter (TLR) assay uses a fluorescent-based system (GFP and mCherry) to determine Homologous Recombination (HR) and Non-homologous End Joining (NHEJ) efficiencies, upon induction of a double-strand break (DSB). Stable LMS-TLR (STS39-TLR, STS137-TLR, and STS210-TLR) and control cell lines (RPEΔp53-TLR and RPEΔp53ΔBRCA1-TLR) containing a single copy of the TLR I-SceI target site were generated. An I-SceI tagged with BFP was introduced to evaluate repair efficiencies. Repair of the DSB by HR generates distinct fluorescent signals (GFP+), compared to NHEJ (mCherry+). LMS cell lines demonstrate HR-deficiency comparable with the RPEΔp53ΔBRCA1 control cell line. In contrast, GFP+ cells were detected in the HR proficient RPEΔp53 cell line. e The bar plot illustrates quantification of GFP to mCherry signal in each LMS cell line and controls. Intact HR (GFP+) is 6X higher in RPEΔp53, compared to LMS cell lines or the RPEΔp53ΔBRCA1 control. Data are derived from eight cell lines examined over three independent experiments and the error bars represent the standard deviation. Source data are provided as a Source Data file.

Full size image

Functional validation of DNA damage repair and homologous recombination deficiency (HRD) in LMS

To validate if SBS3, ID6, and increased levels of DNA repair and damage pathway dysregulation identified in LMS primary tumors and metastasis represent a targeted therapeutic opportunity, we functionally tested a panel of LMS cell lines (five LMS primary lines, three ATCC) to evaluate their sensitivity to DNA damage and PARP inhibitors (Fig. 2b, c, Supplementary Fig. 15). These cell lines are reflective of the diversity of LMS sites of origin, as lines were derived from abdominal, extremity, and gynecological sites (Supplementary Data 8). First, we tested the efficacy of other DDRi to inhibit LMS cell growth. These inhibitors target key kinases that play a pivotal role in the DNA damage response pathway. All LMS cell lines were responsive to the CHK1 inhibitor, Prexasertib HCl (LY2606368). Similarly, 7/8 LMS cell lines were responsive to the WEE1 inhibitor, Adavosertib (AZD1775), while only 5/8 LMS cell lines were partially or fully responsive to the ATR inhibitor, Ceralasertib (AZD6738) (Fig. 2b). To determine the effectiveness of PARPi in LMS, we compared LMS cell line responsiveness to primary UPS cell lines (n = 5, see Supplementary Data 9) and controls: (1) a CRISPR TP53 deletion (RPEΔp53, HR-intact) cell, (2) a TP53 and BRCA1 deletion (RPEΔp53ΔBRCA1, HR-deficient34), and (3) Hs 789.Sk. LMS cell lines were then treated with 24-point concentrations of talazoparib and olaparib (0.013–10 μM) in three independent experiments. All LMS cell lines were highly responsive to the PARP trapper, talazoparib, with a median EC50 of 0.37 μM (range 0.04–0.8 μM), with the most potent inhibition occurring in gynecological LMS cell compared to soft tissue LMS cell lines (0.055 μM vs 0.51 μM) (Fig. 2b, c). In contrast, most LMS cell lines were not as responsive to olaparib, given our strict threshold criteria (EC50 < 1 μM is considered responsive). Taken together, it is evident that there are many promising therapeutic avenues for LMS patients that largely rely on LMS cells’ inability to repair DNA efficiently.

To directly investigate whether HR was defective in LMS cell lines, we used a nuclease-induced genome engineered reporter system (Traffic Light Reporter (TLR) assay) to monitor HR and non-homologous end joining (NHEJ) activity in response to DNA damage35,36. In brief, if double-strand break DNA repair occurs via an HR-dependent mechanism in this assay, a GFP open reading frame is restored. Conversely, if NHEJ is the operative DSB repair mechanism, a frameshift places a mCherry coding sequence in-frame. In LMS cells, there was a marked increase in mCherry signal in response to DNA damage, similar to the RPEΔp53ΔBRCA1 HR-deficient control (Fig. 2d, e). Given the near-universality of LMS cell line responsiveness to single agent talazoparib, a PARP trapper, and to WEE1 and CHEK1 inhibitors, these DDRis are promising therapeutic agents for LMS patients.

Genomic mutations highlight early evolutionary divergence of primary and metastatic LMS

We next asked whether an LMS tumor’s subtype is fixed or whether it can change at relapse. In the five primary-distant metastatic relapse tumor pairs with RNA data, all relapses maintained similar gene expression programs to the matched primary tumors. In general, all tumors from the same patient (whether occurring at different locations or time points) maintained their transcriptional subtype (see Supplementary Fig. 4A). This indicates that the LMS subtype is established early during tumor evolution and retained thereafter––even for tumors that arise from different organs separated by many years.

In addition to maintaining their defining molecular subtype, multiple samples from the same patient had shared early missense mutations in recurrent drivers (such as TP53)(Supplementary Fig. 16). To time the remaining mutations across the full cohort (70 genomes), we took advantage of the fact that many LMS undergo whole-genome duplication (WGD), thus providing a common temporal landmark to which we could compare. We detected WGD in 31/70 (44%) LMS whole-genome sequences. Using these WGD cases, the principles from Jolly and Van Loo37 and MutationTimeR38, a tool used to time somatic mutations relative to clonal and subclonal copy number states (see Methods), we found that only 37% of substitutions arise before WGD, whereas the remaining are late events (Supplementary Fig. 17). Similarly, we found that putative driver deletions preceded WGD, including loss-of-heterozygosity (LOH), copy-neutral LOH, copy-gain LOH events and homozygous deletions. These early arm-level or whole-chromosome deletions involved well-known tumor suppressor genes, such as RB1 (93%), BRCA2 (87%), CDH1 (67%), FANCA(60%), TP53 (87%), BRCA1 (42%), and PTEN (71%) (Supplementary Fig. 18A, B, Supplementary Data 4). Therefore, in addition to somatic TP53 substitutions, which is likely the first alteration in the genome (Supplementary Fig. 17A, C), tumor suppressor chromosomal losses are early and initiating events in leiomyosarcomagenesis, which precede DNA amplifications (Supplementary Fig. 18D).

After having established their subtypes and acquired canonical mutations, subclones within LMS tumors begin to diverge. This was particularly apparent when considering somatic structural rearrangements, a mutation type that is missed by exome or targeted sequencing. While >60% of clonal substitutions and indel mutations were shared between primary and metastatic tumor pairs, only 20% of structural rearrangements were common, suggesting the majority of these events occur after divergence (Fig. 3a, Supplementary Figs. 19A). The high frequency of clonal mutations unique to primary or metastatic tumors indicated that they evolved in parallel. This was further supported by the presence of regional hotspots of mutations, also known as kataegis33, present distinctly in either primary or metastatic tumors (Fig. 3b). First described in breast cancer33, this pattern of clustered substitutions has not yet been reported in LMS (Supplementary Data 5). Furthermore, one tumor (Ab17) also exhibited separate chromothriptic events between primary and both metastatic relapses (Supplementary Fig. 19B, C). Together, these data indicate that, after having acquired key drivers, LMS tumors undergo a late burst of additional changes in the form of WGD, kataegis, and/or chromothripsis. These events occur separately between primary and metastatic pairs, suggesting continued genetic diversification of both the primary and relapse after seeding of the metastases.

a The clinical course of patient Ab17 with a primary (Dx) and two metastatic relapses (MR1 and MR2) is shown (far left, n = 3 samples). Structural variant (SV) overlaps (middle, top) and the cancer cell fraction (CCF) of single-nucleotide variants (SNVs) (middle, bottom) illustrate that there are many SVs and clonal variants that arise independently in the primary tumors and metastatic relapses. Phylogenetic reconstruction of Ab17’s tumors can be seen on the far right. The founder clone harbors a pathogenic TP53 substitution, whole-genome duplication (WGD), as well as loss-of-heterozygosity (LOH) events encompassing TP53 and RB1. The color of each circle represents a distinct clone population. The clonal trajectory and final composition are shown per sample. Branch lengths are proportional to Treeomics mutation assignments. b Rainfall plots of patient Ab17 in diagnosis and their first metastatic relapse illustrate differential kataegis events at different chromosomes between the two time points. Targeted sequencing data were used to confirm kataegis events were unique to each specimen. c The clinical course for patient Ab12 is depicted, which involved no prior treatment and only surgery. The primary specimen at diagnosis (Dx) was located in the inferior vena cava. The tumor was bisected and punch-hole biopsied in three physically distant multiregion (MuRe) locations (n = 3 samples). The phylogenetic reconstruction of this tumor is shown on the right of the schematic and a photo of the resection. The founder clone harbors a pathogenic TP53 substitution, as well as LOH events encompassing TP53, RB1, and PTEN. Larger circles represent major clones, whereas smaller circles represent subclones. The color of each circle represents a distinct clone population. The clonal trajectory and final composition are shown per sample. Branch lengths are proportional to Treeomics mutation assignments, except for clones 8,10 and 6,7 where DPClust mutation assignments were used to stratify the sample.

Full size image

Metastatic LMS branches from the primary tumor decades prior to diagnosis

We then examined the absolute timing of mutations in LMS, both within tumors for which multiple regions have been sampled, and between paired primary-relapse tumor pairs. Multiregion sequencing on three tumors from three patients was performed (two primaries and one metastatic relapse), in addition to bulk sequencing of other tumors from the same patients (Figs. 3c and 4). We reconstructed the phylogeny of each patient’s cancer using an established method (Treeomics39), based on the principles that shared mutations are early and private mutations are late. We validated shared or private mutations using our targeted deep sequencing panel, designed to sequence >75% of point mutations, all non-synonymous indels, and structural variants with breakpoints within known cancer genes (Supplementary Fig. 20). In addition, we corroborated our phylogenies and identified subclonal populations within LMS tumors by using a Dirichlet process to model cellular fractions of variants, as previously described40,41 (See Methods, Supplementary Figs. 21–23). In all patients, TP53mutations and chromosomal losses of chr 10 (PTEN), 13 (RB1) and 17 (TP53) were part of the common trunk of the phylogenetic tree, thus supporting our bulk, single tumor sequencing data that these are the earliest events in the genome. We also noted alterations in the Wnt/β-catenin signaling pathway that were part of the trunk of our phylogenies (Fig. 3a and Fig. 4). Strikingly, we found that metastatic relapses branched off early from the primary tumor in both multiregion cases (Fig. 4). This finding was also supported in the (bulk) primary-metastasis pairs. In patient Ab6, the first metastatic relapse (MR1) appears to have branched off before the primary (Dx) and second metastatic relapse (MR2) separate. These data support a model of branching evolution where clones diverge from a common ancestor and evolve in parallel resulting in multiple clonal lineages42.

The clinical courses of two patients with LMS are shown (samples per patient >3). For the phylogenies, larger circles represent major clones, whereas smaller circles represent subclones. The color of each clone represents a distinct clone population. The clonal trajectory and final composition are shown per sample. Branch lengths are proportional to Treeomics mutation assignments. (a) Patient Ab6 was treated with radiofrequency ablation (RFA) and chemotherapy (chemo). They had three tumors at three separate time points (T1, T2, and T3). The primary tumor at diagnosis (Dx) was located in the small intestine, while the first metastatic relapse (MR1) was located in the liver. The second metastatic relapse (MR2) was multifocal and detected in the vastus lateralis (thigh muscle). MR2 was bisected and biopsies were taken from five distinct sections from both foci (Regions/Re 1–5). Following bulk (Dx, MR1) and multiregion (MR2) sequencing, phylogenetic reconstruction can be seen on the right. Early substitutions in TP53, RB1 and CREBBP, as well as LOH events of chromosomes 10 (PTEN), 13 (RB1), and 17 (TP53) are observed in the founder clone of this patient’s tumors. Chromosome 11 kataegis events and an SPEN deletion were common to Dx and MR2, but not MR1. Genome doubling and chromosome X kataegis occurred only in MR2. Metastatic multifocal nodes greatly resemble each other. Using clock-like mutagenesis, the Dx and MR1 diverge approximately 30 years pre-diagnosis in these patients, while Dx and MR2 diverge approximately 25 years pre-diagnosis. (b) Patient Ab11 was treated with radiation therapy (RT) and had two tumors at two separate time points (T1 and T2). The Dx was located at the posterior aspect of the right kidney and inferior vena cava. This tumor was bisected en face and biopsies were taken from four regions. The metastatic relapses (MR1-3) were taken from the liver. Following bulk (T2: MR1-3) and multiregion (T1: Dx, Regions/Re1-4) sequencing, phylogenetic reconstruction can be seen on the right. Much like Ab6, early losses of chromosomes 10, 13, and 17 are observed. Also seen early are AXIN1 and TET2 point mutations, an ATRX deletion, and a RB1 translocation. Chromosome 2 kataegis events are unique to Dx, while chr 6 chromothripsis events are unique to MR1-3. A chromosome12 kataegis event occurs only in 2/3 liver metastases. Using clock-like mutagenesis, Dx and MR1-3 diverge approximately 15 years pre-diagnosis in this patient.

Full size image

To time the divergence of primary and metastatic tumors we established the rate of mutations in LMS and its subtypes (expressed as substitutions per gigabase per year; see Methods and Anderson et al.9). Signature 1 mutations (SBS1) are proportional to chronological age and so can approximate the number of mitoses a cell has undergone during the lineage of cell divisions since the fertilized egg43 (Supplementary Fig. 24). Using SBS1 as a molecular stopwatch, we determined that the primary and relapsed LMS cancers can be estimated to have diverged from one another 10–30 years prior to diagnosis––providing many years to accumulate the private mutations we observed (Fig. 3a). The divergence of multiple metastases from each other also occurred in the patients’ distant past, over a similar time period (10–30 years pre-diagnosis). This leads to the question: how long ago did adjacent tumor foci or adjacent tumor regions––which are separated only by a few centimeters––diverge from one another? As expected, regions at an increased spatial distance were also more genetically disparate. This is particularly evident in patient Ab6, from whom we sequenced five regions from two tumor foci from their second metastasis (three regions from focus 1 and two regions from focus 2), in addition to both their primary tumor and first metastasis (Fig. 4a). Taken together, from earliest to latest divergence: ~30 years prior to the initial diagnosis of a large abdominal LMS, Ab6’s primary tumor and first metastatic relapse (MR1) diverged; five years later (~25 years prior to initial diagnosis), the ancestor of metastatic relapse 2 (MR2) then diverged. After having established itself in the patient’s vastus lateralis, MR2 diverged again ~6–7 years pre-diagnosis leading to two independent foci located only ~1.5 cm from each other; and then, within each parallel evolving focus, further genetic diversification led to heterogeneous tumor regions separated by <1 cm apart from one another, which diverged ~2–4 years after initial diagnosis. A similar analysis in patient Ab11 confirmed the early origins of LMS and the early divergence of its lethal relapse (Fig. 4b). Thus, both multiregion and paired primary-metastatic relapse sequencing establish the early and parallel genomic evolution of LMS cancers.

By taking a broad view of LMS using both whole-genome and RNA sequencing of primary, metastatic and multi-site sampled tumors, we identified the genomic basis of molecular subtypes, found clinically relevant mutational signatures, and outlined the disease’s evolutionary trajectory.

The overarching theme from this study is that the seeds of a LMS cancer’s aggressiveness are planted very early. The tumor’s originating tissue type, either digestive, vascular, and/or gynecological smooth muscle, combined with its degree of dedifferentiation, is one of the primary determinants of LMS molecular subtypes and the patient’s ultimate survival. Overall, the three subtypes we describe bear similarity to the transcriptomic subgroups previously described11,13, although a meta-analysis is warranted to assess the effect of different NGS platforms. Thus, there is a framework emerging from LMS multiomic efforts and ongoing collaborative initiatives that will generate consensus molecular definitions of LMS subtypes, including the resolution of nomenclature, clinical information, and biomarkers that will aid patient stratification and/or prognostication.

In this study, we describe two LMS subtypes with poorer survival, subtypes 1 (dedifferentiated) and 3 (primarily gynecological), that are also defined by a prevalence in somatic dystrophin deletions and high immune infiltration. Consistent with this, dystrophin deletions are also found in non-myogenic cancers at higher or similar frequencies to other well-known tumor suppressor genes and associated with significantly poorer overall survival44, thus highlighting dystrophin’s role as an emerging tumor suppressor. Importantly, these subtypes have the potential to help differentiate the prognosis of tumors from the same anatomic site. Implementing LMS molecular subtyping, including matching tumors to normal smooth muscle tissue, can be used to reveal its originating tissue at diagnosis. Specifically, an LMS tumor in our cohort resected from an extremity site, clustered with both the subtype 1 group and normal gynecologic tissue, indicating it was in fact a uterine metastasis, which correlated directly with the patient’s clinical history. The ability to define the origins of this “extremity” tumor may provide different, more effective therapeutic options for this patient. This approach could also aid in resolving discrepancies when the previous pathology was inconclusive and unable to document a primary site of LMS. This is, for example, an ongoing conundrum with the diagnosis of uterine leiomyomas45.

The molecular subtypes of LMS, as defined by the originating tissue types and dedifferentiation, then determine the overall burden of secondary mutations that are acquired (across all variant classes)––with the most dedifferentiated tumors and gynecological LMS acquiring the highest burden. Once the cancer is clinically diagnosed, our data suggest that many patients will already have another tumor growing in parallel. Again, one sees that the seeds of LMS are established much earlier––in LMS, lethal metastatic disease can arise up to three decades prior to clinical diagnosis, consistent with other adult cancers38. The exact timing of primary-relapse divergence requires validation in larger cohorts of LMS that include all molecular subtypes and anatomic sites.

With its molecularly defined background established, additional genetic events are layered on top of the LMS tumor genome. Early mutations in TP53 appear to be near-universal in LMS. The LMS founder clone is also characterized by arm-level or whole-chromosome LOH events, most frequently involving the RB1 and TP53 genes. Other patient-specific alterations also arise early, including ATRX deletions and Wnt/β-catenin alterations. From then, mid-late events account for the majority of LMS mutations, contributing to the tumors characteristic karyotypic instability––including dramatic genome-wide duplications, kataegis, and chromothripsis. These events generate widespread genetic diversity, as indicated by the presence of multiple subclonal and clonal populations within a single primary or relapse LMS tumor.

Regardless of molecular subtype, more effective therapies are urgently required in LMS, especially in patients with advanced disease. We show and validate that simple substitution combined with indel signatures reflect HR-deficiency in LMS, which can be exploited therapeutically in the form of PARP trappers and DNA damage response inhibition therapy. Together, HR-deficiency signatures were present in 64% of LMS, lower than previously reported by exomes5. Furthermore, we demonstrate that LMS cell lines are highly responsive to talazoparib treatment, in contrast to other sarcoma and HR proficient cell lines. Here, we use stricter sensitivity thresholds to previous reports, thus further refining the sensitivity of LMS cell lines to PARPi5. Additionally, LMS cell lines are responsive to DNA damaging agents CHK1 and WEE1. The present study suggests that LMS patients should be considered for comprehensive HR profiling for use of mutational signatures as biomarkers of therapy.

To conclude, LMS, the catch all-label for any primary malignant smooth muscle neoplasm, encompasses at least three distinct entities. We report that the utility of DNA damage response inhibitors is promising and warrants further exploration. Finally, the molecular evolution and timing of LMS provide a broad window for early intervention and risk stratification for these patients.