Genomic, Epigenomic, and Transcriptomic Inter- and Intratumor Heterogeneity in Desmoid Tumors

Desmoid tumors account for ∼3% of soft-tissue tumors and are most commonly diagnosed at the median age of 37 to 39 years (1). Most desmoid tumors are sporadic, but a small subset is associated with familial adenomatous polyposis (1). Desmoid tumors usually occur in the extremities, abdominal wall, and abdominal cavity (2). Histologically, desmoid tumors consist of a bland proliferation of spindle cells resembling fibroblasts, lacking nuclear pleomorphism and regional morphologic variability seen in high-grade sarcomas. Despite their bland cytologic appearance, desmoids aggressively invade surrounding tissues. On MRI, desmoid tumors typically show homogeneous isointensity on T1-weighted images and heterogeneously high signals on T2-weighted or short tau inversion recovery (STIR) images (3). Desmoid tumors have a high rate of local recurrence even after complete resection (2, 4). Treatment options include radiotherapy, chemotherapy, and surgical resection (5). Recently, the γ-secretase inhibitor nirogacestat was approved as a new pharmaceutical therapy for desmoid tumors (6).

Desmoid tumors are considered neoplastic monoclonal proliferations based on chromosome X inactivation patterns (7–10). However, it remains unknown whether they undergo clonal evolution, as seen in malignant high-grade tumors. Prior studies indicated mosaicism in desmoid tumors, with the variable ratio of CTNNB1-mutant neoplastic cells and normal fibroblasts (11, 12). Mutations in the CTNNB1 gene, mainly in exon 3 at T41A and S45F/S45P, play a key role in desmoid tumorigenesis and are generally mutually exclusive (12, 13), although rare “double-hit” mutations have been reported in two cases of desmoid tumors (14, 15). Also, nuclear immunoreactivity for β-catenin protein that results from CTNNB1 mutations can vary in distribution and intensity within a tumor (13, 16). Together, these observations suggested clonal heterogeneity within individual desmoid tumors. In addition, recent transcriptomic studies identified distinct molecular subtypes of desmoid tumors and proposed prognostic transcriptomic signatures to predict recurrence (17–20). However, neither of these signatures were introduced into clinical practice, and identifying markers predictive of clinical outcomes and recurrence after surgery remains a subject of active investigation.

Here, we sought to develop a new prognostic signature for desmoid tumors and validate the previously published signatures (19, 20). This analysis revealed substantial variability in the expression of distinct transcriptomic signatures within individual tumors both in patients who did and did not progress. This observation led us to a comprehensive multiomic investigation of molecular intra- and intertumor heterogeneity of desmoid tumors.

We obtained archival formalin-fixed, paraffin-embedded (FFPE) specimens of primary and recurrent desmoid tumors from 20 patients (Table 1). These patients were treated at Stanford Health Care and Vanderbilt University Medical Center. The study was approved by the Stanford University Institutional Review Board (IRB-39030). Written informed consent was obtained from the patients, and the study was conducted in accordance with recognized ethical guidelines (US Common Rule).

Forty-four FFPE tumor specimens were used for extraction of DNA and total RNA (Tables 1 and 2). Diagnosis was confirmed by a surgical pathologist (M. van de Rijn). Multiple 2-mm cores from FFPE tumors were re-embedded horizontally, and sections were reexamined by light microscopy to select cases that contain at least 90% lesional spindle cells. DNA and RNA were extracted from re-embedded tumor cores using AllPrep DNA/RNA FFPE Kit including RNase A and DNase treatment steps, respectively, according to the manufacturer’s protocol (Qiagen).

Seventy-five nanograms of DNA from 23 desmoid tumor specimens from three patients was used for genome-wide copy-number and allelic ratio profiling using the OncoScan FFPE Assay performed by Affymetrix Research Services Laboratory (Affymetrix). Results were visualized and analyzed using Nexus Express for OncoScan3 software (BioDiscovery, RRID: SCR_004557) with the SNP-FASST2 algorithm. All specimens had the median absolute pairwise difference value below 0.3 and the ndSNPQC value above 26, indicating good quality of data (Supplementary Table S1). Copy-number gain and loss were defined as the log₂ ratio >0.2 and <−0.2, respectively (Supplementary Table S2). A minimum of 20 probes in the segment were used to call copy-number alterations. DNA copy-number alterations reported in the Database of Genomic Variants in the general population were excluded (RRID: SCR_007000, http://dgv.tcag.ca/dgv/app/home, accessed March 22, 2023; ref. 21). Heatmaps representing copy-number alterations in each patient were visualized using the ClustVis package in R (version 0.0.0.9000, RRID: SCR_017133; ref. 22).

Whole-exome sequencing libraries were prepared from DNA extracted from 22 desmoid tumor specimens from three patients and from patient-matched normal tissues as germline controls. Normal tissues not involved by tumor were selected by a surgical pathologist (M. van de Rijn), which included the intestinal mucosa, skin, and a scar for patients 1, 2, and 3, respectively. One hundred nanograms of DNA was sheared (Covaris S2) to ∼250-bp fragments. Libraries were prepared using KAPA LTP Library Preparation Kit (KAPA Biosystems). Hybridization-based enrichment was performed using the SeqCap EZ Human Exome v3.0 Library (Roche) at 47°C for 64 to 72 hours. Pooled indexed libraries were sequenced using 2 × 101 bp mode on HiSeq 4000 (Illumina). Reads were aligned to the reference genome (GRCh37/hg19) using the BWA-MEM algorithm (BWA version 0.7.13, RRID: SCR_010910) with default settings (23). SAMtools (version 1.3.1, RRID: SCR_002105) was used for sorting and indexing (24). Picard (version 1.141, http://broadinstitute.github.io/picard, RRID: SCR_006525) was used to remove duplicate reads. GATK (version 3.5.21, RRID: SCR_001876) was used for local realignment and base-call recalibration (25). Detailed metrics are included in Supplementary Table S1. Somatic single-nucleotide variants and indels in tumor specimens that were absent in patient-matched normal DNA were identified using Mutect (version 1.1.7, RRID: SCR_000559; ref. 26). Variants were annotated using a Variant Effect Predictor (version 109, RRID: SCR_007931; ref. 27). Variants reported in 1000 Genomes (Phase 3, RRID: SCR_006828) and gnomAD (version r2.1, RRID: SCR_014964) databases were excluded from the analysis. Nonsynonymous missense, frameshift, stop-gain, stop-loss, and start-loss variants affecting the protein-coding and splicing regions were retained in the analysis.

To reconstruct phylogenetic trees, we used all somatic CTNNB1 mutations regardless of the coverage, and all other nonsynonymous variants identified with ≥50 read coverage in the tumor. We required 0% variant allele fraction (VAF) in the patient-matched normal tissue and ≥10% VAF in the tumor. With these criteria, we identified a median of 149 somatic mutations per specimen (range: 1–594; Supplementary Table S3). These variants were used for phylogenetic tree reconstruction using Lineage Inference for Cancer Heterogeneity and Evolution (LICHeE; version 1.0; ref. 28).

We used the OncoPrinter tool on cBioPortal version 5.4.2 (cbioportal.org/oncoprinter, RRID: SCR_014555) to visualize somatic nonsynonymous mutations that occurred in at least two distinct specimens in patients 1 and 2 (29).

Five hundred nanograms of DNA from 20 desmoid tumor specimens from three patients was profiled using Infinium MethylationEPIC microarrays (Illumina). We required that >90% of probes were measured with P < 0.05 on at least three beads [detectionP function, minfi package (30) version 1.32.0, RRID: SCR_012830; and beadcount function, wateRmelon package (31) version 1.30.0, RRID: SCR_001296; Supplementary Table S1]. We also required log₂ median intensity of the methylated and unmethylated channels >9 (shinySummarize function, minfi package; ref. 30). The data were normalized using Illumina preprocessing (preprocessIllumina function, minfi package; ref. 30). We removed poor quality probes (measured with P ≥ 0.01), probes associated with SNPs (dropLociWithSnps function, minfi package; ref. 30), cross-reactive probes (xreactive_probes function, maxprobes package version 0.0.2, https://github.com/markgene/maxprobes), and probes associated with sex chromosomes. Exploratory analysis was done using 1% of probes with the most variable β values. Heatmaps that visualize unsupervised clustering (Euclidean distance and Ward linkage) were plotted using the ClustVis package in R (version 0.0.0.9000; ref. 22).

Whole-transcriptome sequencing [RNA sequencing (RNA-seq)] was performed on 44 desmoid tumor specimens from 20 patients. Libraries were prepared from 1 μg of total RNA using TruSeq Stranded Total RNA Library Prep Kit (Illumina). Pooled indexed libraries were sequenced on HiSeq 2000 (Illumina, RRID: SCR_020130) in 2 × 101 bp mode. The RNA-seq statistics are summarized in Supplementary Tables S1 and S4. Reads were mapped to the reference genome (GRCh37/hg19) using Rsubread in R (version 1.26.1, RRID: SCR_016945; ref. 32). The transcripts assigned to 19,727 protein-coding genes were quantified using featureCounts function (Rsubread package version 1.26.1). The data were normalized with variant stabilizing transformation using the DESeq2 package in R (version 1.26.0, RRID: SCR_015687; ref. 33). Exploratory analysis was done using normalized expression values of the 10% most variable genes (Supplementary Tables S5–S7). Heatmaps that visualize unsupervised clustering (centered and scaled normalized counts, Euclidean distance, and Ward linkage) were plotted using the ClustVis package (22). The FactoMineR package in R (version 2.10, RRID: SCR_014602) was used to confirm that the main source of variance in RNA-seq data was not associated with any technical features. Enriched pathways were identified using the Molecular Signatures Database (MSigDB, RRID: SCR_016863) and the Gene Set Enrichment Analysis (www.gsea-msigdb.org, RRID: SCR_003199). RNA-seq reads were used to validate clonal and subclonal mutations identified by exome sequencing. Sorted and indexed BAM files from RNA-seq data were inspected in the Integrative Genomics Viewer (IGV; version 2.14.1, RRID: SCR_011793; ref. 34). For validation, we required >10 read coverage and the variant to be present in both forward and reverse reads.

Five-micrometer whole sections of FFPE tissues or tissue microarrays containing duplicate 1-mm cores from each tissue were used for hematoxylin and eosin staining and IHC staining by the avidin–biotin–peroxidase complex method (antibodies listed in Supplementary Table S8). Appropriate positive and negative controls were analyzed in parallel. The intensity of staining was evaluated at 400× magnification.

Sanger sequencing was performed to detect CTNNB1 mutations in all specimens. Primers were designed using Primer3 (forward: 5′-GAC​ATG​GCC​ATG​GAA​CCA​GA-3′ and reverse: 5′-CCC​ACT​CAT​ACA​GGA​CTT​GGG-3′; ref. 35). PCR was performed on 50 ng of DNA. PCR products were obtained using AmpliTaq Gold DNA Polymerase (Life Technologies), purified using QIAquick PCR Purification Kit (Qiagen), and sequenced using BigDye Terminator v3.1 (Life Technologies) on the 3130 × l Genetic Analyzer (Life Technologies) in the forward and reverse directions. 4Peaks software (Nucleobytes, RRID: SCR_000015) was used for the analysis. CTNNB1 mutations were further verified using RNA-seq data in the IGV. Two patients had familial adenomatous polyposis syndrome, and the APC mutations were determined in their tumors using RNA-seq data in the IGV.

Sequencing data were deposited in the Sequence Read Archive (exome sequencing data: PRJNA1052459; RNA-seq data of 24 specimens from 3 patients: PRJNA1052707; and RNA-seq data of primary tumors from 20 patients: PRJNA1099876). Microarray DNA copy-number data and methylation data were deposited in the Gene Expression Omnibus (GSE250149 and GSE250403, respectively). Analysis was done in R version 3.6.3 (RRID: SCR_001905). The code for data analysis is available at: https://github.com/przybyllab/dtf. The data used to plot each figure are available from the corresponding author upon request.

We performed RNA-seq profiling of 31 specimens from 20 primary desmoid tumors to identify genes that may help predict relapse after surgery. Thirteen of these patients developed local recurrence after surgical resection of the primary tumor. We profiled between 2 and 6 different regions of the same primary tumor from six patients to examine the possible intratumor heterogeneity of transcriptomic signatures. We identified 14 differentially expressed genes between patients who did and did not develop local recurrence (Supplementary Table S9). However, we noticed that the expression of these genes was highly variable within the compared groups of patients and in different regions within the same tumors (Fig. 1A).

Next, we examined whether our transcriptomic data can be used to validate the gene expression signatures associated with progression-free survival that were identified in two previously published studies (19, 20). Salas and colleagues identified a 36-gene prognostic signature based on microarray studies, and Kohsaka and colleagues identified a different three-gene prognostic signature based on RNA-seq profiling. The comparison of patient cohorts included in our study and in the studies published by Salas and colleagues and Kohsaka and colleagues is included in Supplementary Table S10. In the signature proposed by Salas and colleagues, 18 transcripts were upregulated in primary desmoid tumors that recurred (genes associated with poor outcomes) and 18 transcripts were upregulated in primary tumors that did not develop recurrence (genes associated with favorable outcomes). Of these 36 genes, 34 genes could be measured in our dataset (Fig. 1B). Among the three prognostic genes identified by Kohsaka and colleagues, two genes could be measured in our dataset (IFI6 and LGMN; Fig. 1C). We observed remarkable interpatient and intratumor heterogeneity of expression of genes included in these two prognostic signatures both in patients who did and did not progress (Fig. 1B and C).

From the 36-gene signature, in our dataset, only the FECH gene showed significant differential expression between primary tumors that did and did not recur (unpaired two-sided t test, P = 0.01). In the study of Salas and colleagues (20), this gene showed higher expression in primary tumors that recurred. In contrast, in our dataset, the FECH gene showed higher expression in tumors that did not recur (Fig. 1B). We also observed highly variable expression of selected genes in different regions of primary tumors of patients 1, 2, 3, 6, 7, and 20. For example, in patient 1 who developed recurrence, we observed highly variable expression of genes associated with both favorable prognosis (e.g., TFG and ST7L) and unfavorable prognosis (e.g., ENY2, RAN, and ARF6) across six regions of a single primary tumor (Fig. 1B).

In the study of Kohsaka and colleagues (19), elevated IFI6 and LGMN expression was associated with high risk for recurrence, but in our dataset, these genes were not differentially expressed between tumors that did and did not recur (unpaired two-sided t test, P = 0.8 and P = 0.4, respectively; Fig. 1C). We also noted substantial intratumor heterogeneity of expression of these genes (e.g., IFI6 in patients 2 and 20 and LGMN in patient 1; Fig. 1C).

We validated variable expression of these markers in desmoid tumors also at the protein level by IHC using a tissue microarray that contains multiple cores sampled from primary tissue specimens from 13/20 patients analyzed by RNA-seq. Cores for the tissue microarray were taken from tissue directly adjacent to the cores taken for RNA-seq. We selected one marker from our signature (EFL1) and three markers from the study of Salas and colleagues (HDAC6, TFG, and hnRNPC). TFG and hnRNPC were strongly expressed in all specimens (Supplementary Table S11). EFL1 and HDAC6 showed variable expression (Supplementary Fig. S1; Supplementary Table S11), which was not associated with development of local recurrence (two-sided χ² test, P = 0.78 and P = 0.89, respectively). We also observed intratumor heterogeneity of protein expression in different regions of the same tumor (Supplementary Fig. S1).

This analysis demonstrates that prognostic classification based on gene expression signatures and proteins encoded by these genes may vary significantly with the area of the tumor that is analyzed. These findings prompted us to further explore multiomic intra- and intertumor heterogeneity in primary and recurrent tumors of three patients, who had desmoid tumors at different anatomic locations, developed either a single or multiple local recurrences, and underwent different treatment regimens.

Patient 1 was diagnosed with a 19-cm mesenteric desmoid tumor that locally recurred 1 year after surgical resection with positive margins. The recurrent tumor measured 4.5 cm. The patient had not received any adjuvant treatment. We sampled cores from six different regions of the primary tumor and three regions of the recurrent tumor (Fig. 2A).

By exome sequencing, we detected two different somatic mutations in the CTNNB1 gene (Fig. 2B and C; Table 2). All regions of the primary and recurrent tumors carried the T41A mutation, but two regions of the recurrent tumor also carried the S45P mutation, at a lower VAF than the T41A mutation (Table 2). We validated the presence of these distinct CTNNB1 mutations in DNA by Sanger sequencing (Fig. 2D) and in the RNA-seq reads (Supplementary Tables S12 and S13). A phylogenetic tree reconstructed using all nonsynonymous somatic mutations revealed the presence of three distinct subclones in the primary tumor and two distinct subclones in the recurrent tumor that were not present in the primary tumor (Fig. 2B and C). All sampled regions of the recurrent tumor shared a mutation of the PTP4A1 gene, which corresponded with significantly lower RNA expression compared with the sampled wild-type areas of the primary tumor (Mann–Whitney test, P = 0.02; Fig. 2E). Among the clonal and subclonal mutations in patient 1, we validated 23/32 of these mutations in the RNA-seq data, 7/32 variants were not detected in the RNA-seq data (wild-type allele), and the genomic position of 2/32 variants had no coverage in the RNA-seq data (Supplementary Table S12). In addition, inspection of the RNA-seq reads demonstrated that selected subclonal mutations that were detected only in DNA in the recurrent tumor were detectable in RNA in one region of the primary tumor (PTP4A1 and RNPEPL1; Supplementary Table S13).

All sampled regions of primary and recurrent tumors showed normal diploid karyotypes with only few DNA copy-number alterations (Fig. 2F; Supplementary Table S2). DNA methylation and transcriptomic profiles were clearly distinct between the primary and recurrent tumors (Fig. 2G and H). Interestingly, one region of the primary tumor had a divergent transcriptomic profile from the other regions within the same tumor, showing enrichment of 12 genes in the β-oxidation pathway in mitochondria (MSigDB KEGG_OXIDATIVE_PHOSPHORYLATION pathway, adjusted P = 1.5 × 10⁻⁸), which may indicate areas of distinct metabolic activity within the primary tumor. We also noticed variable expression of selected genes of the Notch pathway (e.g., high expression of NOTCH2 in one region of the recurrent tumor; Supplementary Table S5). Among the most variably expressed genes, HTRA3 and TCTEX1D1 showed an inverse correlation between gene expression and methylation in the promoter region (Pearson correlation r = −0.79 and r = −0.73, respectively, P < 0.05), with significantly lower expression and higher methylation in the promoter region of these genes in the primary tumor compared with the recurrent tumor (normalized gene expression: unpaired two-sided t test, P = 0.0128 and P = 0.0102, respectively; DNA methylation: unpaired two-sided t test, P < 0.0001 and P < 0.0001, respectively; Fig. 2I and J).

These molecular profiles showed intertumoral heterogeneity between the primary and recurrent tumors, as well as intratumoral heterogeneity within both the primary and recurrent tumors. Our analysis also demonstrated the coexistence of two distinct mutations in the CTNNB1 gene in selected regions of the recurrent tumor.

Patient 2 was diagnosed with a 6.5-cm desmoid tumor of the buttock that was resected with positive margins and received radiotherapy. Two years later, the patient was treated with imatinib and later had several recurrent tumors excised at 4, 10, and 11 years after initial resection. The patient also received vinblastine and methotrexate between the second and third recurrences. The first recurrence was not available for molecular profiling. We sampled two different regions of the primary tumor, three regions of the second recurrence, and two regions of the third recurrence (Fig. 3A).

All seven samples from the primary and recurrent tumors carried the same CTTNNB1 T41A mutation in the exome sequencing and RNA-seq data (Table 2; Supplementary Tables S12 and S14). The third recurrent tumor showed three new mutations in DNA affecting the MKI67, MYOF, and HAUS1 genes that were not present in the primary tumor and second recurrence (Fig. 3B and C). RNA-seq data confirmed the MYOF mutation in one region of the third recurrence and showed that the subclonal HAUS1 mutation was present also in one region of the second recurrence (Supplementary Table S14). In this patient, the primary and recurrent tumors carried several clonal and subclonal DNA copy-number alterations (Fig. 3D; Supplementary Table S2). Based on DNA methylation profiles, the primary and recurrent tumors formed distinct clusters (Fig. 3E). Transcriptomic analysis demonstrated highly distinct profiles of the primary tumor and the second and third recurrent tumors (Fig. 3F). Interestingly, one region of the second recurrent tumor had a highly similar transcriptomic profile to the primary tumor (Fig. 3F). In conclusion, the multiomic analysis of the primary and recurrent tumors from patient 2 revealed diverging genomic, epigenomic, and transcriptomic programs in recurrent tumors that developed only 1 year apart.

Patient 3 initially underwent an excisional biopsy of a 3.5-cm primary desmoid tumor of the arm with positive margins, followed by radiotherapy. Afterward, the patient developed multiple local recurrences 1, 4, and 5 years after primary resection. The patient was treated with toremifene and meloxicam between the second and third recurrences. From this patient, we sampled three regions of the primary tumor, three regions of the first recurrence, and a single region from each of the second and third recurrences (Fig. 4A).

All but one specimen analyzed by whole-exome sequencing carried the S45F mutation in the CTNNB1 gene (Fig. 4B). We did not detect a CTNNB1 mutation in the single sample taken from the second recurrence, but the sequencing coverage of the genomic position that corresponds to the S45F mutation was lower compared with the other specimens (17 reads compared with a median of 57 reads in the remaining specimens). However, we confirmed the same CTNNB1 mutation in the RNA-seq reads from the same specimen (VAF = 0.45, 204 read coverage; Supplementary Fig. S2). There were no other genes with mutations shared by different regions of the primary or recurrent tumors in this patient (Fig. 4B).

One of the regions of the first recurrence and the second recurrence showed a high number of focal DNA copy-number alterations that were not present in other specimens (Fig. 4C; Supplementary Table S2). DNA methylation analysis demonstrated that the second recurrence had a highly distinct epigenomic profile from all other primary and recurrent specimens obtained from this patient (Fig. 4D). RNA-seq revealed that the primary tumor and the first recurrence showed a high degree of intratumoral heterogeneity at the transcriptomic level (Fig. 4E). The third recurrence showed significant enrichment of 12 genes upregulated in response to hypoxia (MSigDB HALLMARK_HYPOXIA pathway, adjusted P = 3 × 10⁻³; Fig. 4F). Overall, the multiomic analysis of specimens obtained from patient 3 revealed a high degree of intertumoral heterogeneity between subsequent recurrent tumors and substantial intratumoral heterogeneity within the primary and recurrent tumors, exemplified by variable genomic complexity and focal enrichment of the genes upregulated in hypoxia.

This study represents the first attempt to characterize the molecular heterogeneity of desmoid tumors, which are bland fibroblastic tumors that can show local variability in cellularity but no significant nuclear pleomorphism. First, we demonstrated substantial interpatient and intratumor heterogeneity of expression of transcriptomic prognostic signatures in a series of 20 patients. Next, our multiomic and multiregion analyses of primary and recurrent tumors from three patients demonstrated remarkable and unexpected intra- and intertumor heterogeneity in this disease. Heterogeneous patterns of point mutations, DNA copy-number alterations, DNA methylation, or gene expression were identified in tumors from all three patients.

Recognizing molecular heterogeneity of tumors is vital for tailoring treatment, understanding disease progression, and developing effective therapies. Clonal mutations, present in most neoplastic cells within a tumor, are considered driver events most suitable as therapeutic targets. Subclonal alterations that are heterogeneous within tumors are also important for understanding tumor evolution and may contribute to relapse and therapeutic resistance. Intra- and intertumor heterogeneity has been characterized in multiple types of cancer [reviewed by Greaves and Maley (36)]. We also previously demonstrated genomic intra- and intertumor heterogeneity of leiomyosarcoma (37).

Intratumor heterogeneity of desmoid tumors included in our series was evident in the variable expression of three distinct gene signatures associated with disease recurrence. We also confirmed the heterogeneous expression of selected markers from these prognostic signatures at the protein level (Supplementary Fig. S1). This demonstrated that the prognostic gene expression signatures may have limited clinical utility in desmoid tumors because they are highly vulnerable to tumor sampling bias, as previously reported in several types of cancer, including clear-cell carcinoma, non–small cell lung cancer, breast cancer, and hepatocellular carcinoma (38–41). For example, Gyanchandani and colleagues (39) compared the recurrence scores assigned with five different gene expression panels in breast cancer and demonstrated that profiling a single tissue core can under- or overestimate prognostic risk. These authors proposed that assessing multiple representative and atypical tissue cores may help fully account for the heterogeneity-driven variation in risk prediction (39).

There are several key differences between our study and the studies of Salas and colleagues (20) and Kohsaka and colleagues (19) that might have contributed to the discrepancies between the identified gene expression signatures (Supplementary Table S10). We analyzed FFPE tissues, whereas the other studies used fresh frozen tumor specimens. Our study and Kohsaka and colleagues applied RNA-seq, whereas Salas and colleagues used cDNA microarrays. Our cohort and that of Salas and colleagues included only primary untreated tumors, whereas Kohsaka and colleagues had 14/64 (22%) patients treated with NSAIDs as first-line treatment. The anatomic distribution of tumors was similar between our study and that of Salas and colleagues, but Kohsaka and colleagues did not include intra-abdominal tumors. Other patient characteristics were comparable across the cohorts (Supplementary Table S10). The discrepancies may also originate from different statistical methods used to develop prognostic gene signatures. We applied differential expression analysis between patients who did and did not progress after surgery. Salas and colleagues also applied differential expression analysis but only between patients with no recurrence versus those with multiple local recurrences, excluding those with only one recurrence. Kohsaka and colleagues applied univariate and multivariable Cox regression analyses. Although all these methods are valid for identifying prognostic gene signatures, these differences might have contributed to the lack of concordance between them.

Mutations of the CTNNB1 gene are well-established driver events in desmoid tumors (12, 42). Our study confirmed the presence of CTNNB1 mutations in DNA of all but one of the regions of the primary and recurrent tumors. By applying the LICHeE algorithm for phylogenetic inference using somatic single-nucleotide variants (28), we identified a few clonal molecular alterations, other than CTNNB1 mutations, that were shared by all regions of the primary or recurrent tumors. We also discovered new mutations in recurrent tumors, which may be the driver events of relapse (e.g., MKI67, MYOF, and HAUS1 in patient 2). In patient 3, we identified DNA copy-number alterations shared by different regions of the primary and recurrent tumors (present in six of seven specimens), including loss of the CHL1 gene that acts as a tumor suppressor in neuroblastoma and nasopharyngeal carcinoma (43, 44). Although these alterations may be potential driver events in desmoid tumors, none of these alterations are currently actionable by targeted pharmacotherapy.

We also identified subclonal CTNNB1 mutations in different regions of a recurrent tumor from patient 1 (Table 2). The rare instances of “two-hit” mutations in the CTNNB1 gene described in the literature have not been associated with any specific clinicopathologic features or an unusually aggressive phenotype. We also identified subclonal mutations in other genes, including TRIP11, SOAT1, GGCX, TOP3B, and MYH6 genes in the primary tumor of this patient and SETD1A and RNPEPL1 genes in the recurrent tumor. These genes have not been directly implicated in the oncogenic processes in cancer, and their significance in desmoid tumors remains unknown. Among these genes, a missense mutation of the MYH6 gene was previously reported in a single desmoid tumor (12). It remains to be confirmed in future studies whether these mutations occur in a significant subset of desmoid tumors. A possible confounder in our clonality analysis is sampling bias as we did not analyze multiple regions from all recurrent tumors because of lack of available material. Also, our clonality analysis might have been affected by the depth of sequencing and quality of DNA from FFPE specimens.

Overall, desmoid tumors carry limited DNA copy-number alterations. However, we observed an unusually high number of subclonal genomic alterations in specimens from two distinct recurrent tumors from patient 3. One region of the first recurrence showed multiple copy-number losses, and the second recurrence showed multiple copy-number gains. We speculate that selected regions of the recurrent tumors may demonstrate variable, perhaps temporary, activation of mechanisms contributing to genomic instability and disease relapse.

We sought to verify whether the variable number of mutations and DNA copy-number alterations in different regions of primary and recurrent tumors is potentially associated with markers of genomic instability. We analyzed the expression of mismatch repair proteins (MLH1, MSH2, MSH6, and PMS2) in a tissue microarray with specimens of primary and recurrent tumors from patients 2 and 3, which showed intact expression in all specimens (Supplementary Table S15). DNA repair pathways cataloged in MSigDB (Hallmark, Reactome, and WikiPathways collections) were not differentially enriched between tumor regions with varying numbers of genomic alterations (not shown). Additionally, there was no correlation between global DNA hypomethylation and the number of somatic mutations or DNA copy-number alterations (not shown).

We also investigated a possible association between immune cell infiltrates and heterogeneity of desmoid tumors. However, we did not detect mutations in any immune escape genes in the exome sequencing data and did not identify transcriptomic signatures of any immune cells using CIBERSORT (45) in the RNA-seq data (not shown). Thus, we do not suspect a substantial presence of immune cell infiltrates in desmoid tumors.

Interestingly, we identified a distinct extent of heterogeneity in different types of data. Overall, transcriptomic profiles of desmoid tumors showed the highest degree of variability within and between individual tumors. DNA methylation profiles were generally similar between different areas within individual primary or recurrent tumors but showed substantial differences between the primary and recurrent tumors. In addition, in patient 1, we identified two candidate genes (HTRA3 and TCTEX1D1) in which transcription may be differentially regulated between primary and recurrent tumors by methylation of the promoter region. The changes in DNA methylation profiles between primary and recurrence could indicate a possible contribution of epigenomic mechanisms to desmoid tumor progression.

Different regions of the recurrent tumor from patient 1, who did not undergo any adjuvant therapy, showed highly similar profiles in the multiomic analysis. In contrast, we demonstrated divergent molecular profiles of subsequent recurrent desmoid tumors that developed after different modalities of adjuvant therapy in patients 2 and 3. This indicates that desmoid tumors alter their molecular profile over time and after exposure to treatment, which may contribute to therapeutic resistance. These findings suggest that each subsequent tumor recurrence may require different therapy. The Notch pathway is a therapeutic vulnerability in desmoid tumors that can be targeted by γ-secretase inhibitors (6, 46, 47). Interestingly, we observed high expression of the NOTCH2 gene in one area of the recurrent tumor in patient 1 (Supplementary Table S5). This suggests a possible variable expression of the components of the Notch pathway in desmoid tumors, which may potentially affect response to γ-secretase inhibitors. Our findings indicate that clonal evolution may contribute to the variability in the clinical behavior of these tumors.

In summary, our study demonstrates for the first time a remarkable intra- and intertumor heterogeneity of desmoid tumors using a multiomic approach. Our analysis shows that in desmoid tumors, similar to other types of solid tumors, analysis of a single-tumor biopsy may substantially underestimate the extent of molecular alterations and that prognostic transcriptomic signatures for desmoid tumors are vulnerable to sampling bias. Our study also shows that recurrent desmoid tumors acquire multiple molecular alterations that are not present in primary tumors. This molecular heterogeneity is an important consideration in drug development and in validation of prognostic and predictive biomarkers for desmoid tumors.

J. Cates reports service on the scientific advisory board of Eluciderm, Inc. There are no relevant conflicts of interest or patents relevant to this submission. B. Haibe-Kains reports personal fees from Code Ocean, Consortium de recherche biopharmaceutique, Break Through Cancer, Commonwealth Cancer Consortium, and Cancer Grand Challenges outside the submitted work. No disclosures were reported by the other authors.

C. De Bellis: Formal analysis, visualization, writing–original draft. S. Vennam: Data curation, formal analysis, methodology, writing–review and editing. C. Eeles: Formal analysis, methodology, writing–review and editing. P. Rahimizadeh: Validation, methodology, writing–review and editing. J. Cates: Resources, data curation, funding acquisition, writing–review and editing. T. Stricker: Resources, data curation, writing–review and editing. J. Hoffman: Resources, data curation, writing–review and editing. K. Ganjoo: Resources, data curation, writing–review and editing. G.W. Charville: Formal analysis, validation, methodology, writing–review and editing. B. Haibe-Kains: Formal analysis, methodology, writing–review and editing. M. van de Rijn: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, methodology, writing–original draft, project administration, writing–review and editing. J. Przybyl: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

This research was supported by research grants from The Desmoid Tumor Research Foundation, Inc. (J. Cates, M. van de Rijn, and J. Przybyl) and The Desmoid Tumour Foundation of Canada (J. Przybyl).