Genomic characterization of metastatic ultra-hypermutated interdigitating dendritic cell sarcoma through rapid research autopsy

Interdigitating dendritic cell sarcoma (IDCS) is an extremely rare cancer of dendritic cell origin that lacks a standardized treatment approach. Here, we performed genomic characterization of metastatic IDCS through whole exome sequencing (WES) of tumor tissues procured from a patient who underwent research autopsy. WES was also performed on a treatment-naïve tumor biopsy sample obtained from prior surgical resection. Our analyses revealed ultra-hypermutation, defined as >100 mutations per megabase, in this patient's cancer, which was further characterized by the presence of three distinct mutational signatures including UV radiation and APOBEC signatures. To characterize clonal heterogeneity, we used the bioinformatics tool Canopy to leverage single nucleotide and copy number variants to catalog six subclones across various metastatic tumors. Truncal alterations, defined as being present in all clonal tumor cell populations, in this patient's cancer include point mutations in TP53 and CDKN2A and amplifications of c-KIT and APOBEC3A-H, which are likely driver mutations. In summary, we have performed genomic characterization evaluating tumor mutational burden (TMB) and heterogeneity in a patient with metastatic IDCS. Despite ultra-hypermutation, this patient's cancer was not responsive to treatment with PD-1 inhibition. Our results underscore the importance of characterizing clonal heterogeneity in TMB-high cancers.


INTRODUCTION
Interdigitating dendritic cell sarcoma (IDCS) is an extremely rare malignancy of dendritic cell origin with approximately 100 cases reported to date [1,2]. Due to its rarity and challenging diagnosis, genomic characterization of this neoplasm has not been previously reported. Furthermore, no standard therapy exists for IDCS, which Research Paper www.oncotarget.com tends to affect middle-age adults with median age of diagnosis of 56.5 years [3]. Localized IDCS constitutes 47% of cases and manifests as painless lymphadenopathy, most commonly involving the cervical and axillary nodes. Isolated extra-nodal disease occurs in 25% of cases, involving the liver, lung, spleen, bone marrow and gastrointestinal tract [3]. Distant metastases occur in 39% of cases and most frequently involved lymph nodes, lung, liver and bone marrow [3].
Here we performed whole exome sequencing (WES) of multiple tumors obtained through rapid research autopsy of a patient with metastatic IDCS. To our knowledge, this is the first time IDCS has been extensively sequenced and analyzed. Our findings demonstrate a rare ultra-hypermutated genotype, with >100 somatic mutations per megabase of genome (mut/Mb) [4], characterized by distinct mutational signature profiles and clonal diversity that occurred early in the evolution of this patient's cancer. Phylogenetic analysis revealed inactivation of TP53 and CDKN2A (p16 INK4a /p19 ARF ) as well as amplification of c-KIT and PDGFRα as truncal alterations. We further detected amplification of APOBEC3A-H encoding cytidine deaminases on chromosome 22q that may have attributed to the ultrahypermutation. In summary, our work describes the genomic landscape and clonal heterogeneity of an ultra-rare cancer through the innovative approach of research autopsy.

Diagnosis and treatment
Our patient was a 57-year-old Caucasian male who developed an isolated FDG-avid 2.4 × 2.4 cm right-sided neck mass on PET scan in 2016. Biopsy result of this mass showed a "pleomorphic/spindle cell neoplasm". Immunohistochemical (IHC) analyses demonstrated focal CK7 staining. Clinical evaluation revealed no primary lesion involving the skin or oropharynx. In August 2016, he underwent modified radical right neck dissection and partial submandibular gland excision with one level I lymph node demonstrating complete tumor involvement and suspicious for extranodal extension; remaining lymph nodes (0/29) from levels II, III and IV were negative for tumor infiltration. IHC on the surgical specimen showed positivity for S100, SOX10 and CK7. A diagnosis of IDCS was issued. Given his localized presentation, he received adjuvant radiation to the right level IB-III nodes with a boost to the primary site of disease. In November 2016, he initiated adjuvant nivolumab (3 mg/kg every 2 weeks) given FoundationOne ® report showing >100 mut/Mb in his tumor. Following adjuvant immunotherapy, he developed metastatic recurrence that failed to respond to subsequent treatments including combination immunotherapy (CTLA-4/PD-1 inhibitors), chemotherapy and molecularly targeted therapies ( Figure 1). His clinical course leading up to the research autopsy is summarized in Figure 1A.

Rapid research autopsy
Prior to his death, the patient was consented to an Institutional Review Board (IRB)-approved research autopsy study for patients with advanced cancers (Figure 2). CT scans obtained prior to hospice enrollment showed tumor involvement of multiple organs ( Figure 1B) and were used to guide sample procurement at time of research autopsy, which was performed within three hours postmortem. A total of twenty-four metastatic tumor samples were procured from involved organ sites ( Figure 1B

Genomic characterization of IDCS
We performed WES on nine metastatic tumors and the resected pre-treatment tumor ( Table 1). Calculation of tumor mutational burden (TMB), including single nucleotide variants (SNVs) and insertions/ deletions (indels), demonstrated ultra-hypermutation (130.1-167.0 mut/Mb) in all tumors analyzed (Table 1;  Supplementary Table 1). Interrogation of microsatellite status utilizing MANTIS [5] showed that all tumors were microsatellite stable (MSS) despite their ultrahypermutation (Table 1). To further characterize this patient's cancer, we investigated whether specific substitution mutational signatures, which arise due to different mutagenic processes [6], may be present. This analysis revealed the presence of Signature 2 (APOBEC), Signature 7 (UV light) and Signature 11 (alkylating agent), which are all characterized by C>T substitutions, in all tumor samples (Figure 3A-3B; Supplementary Figure 1 and  Supplementary Table 2). Of the three mutational signatures, Signature 7 had the highest prevalence ( Figure 3A). We next classified the 7,417 somatic variants (SNVs and indels) in all ten tumor samples of this patient into three categories: Ubiquitous, Shared and Private ( Figure 4A-4B). Ubiquitous variants are variants detected in all tumor samples and included driver mutations in TP53 and CDKN2A, while shared variants are those detected in subsets of tumor samples. Finally, private variants are variants unique to each tumor sample. A tumor-centric evolutionary tree was constructed based on the pairwise genetic distance, as measured by number of discordant SNVs, between the wildtype cells (normal) and tumor samples in this patient ( Figure 4C). Finally, we interrogated the prevalence of copy number variations (CNVs) using the allele-specific CNV caller FALCON [7]. After manual review and curation, 29 unique CNVs including amplification (including c-KIT, PDGFRα, and APOBEC3A-H), deletion and loss-ofheterozygosity events were detected (Supplementary

Reconstructing clonal evolution
Integrating the above genomic data (SNVs, indels and CNVs), we used the program Canopy [8] to synthesize a hypothetical model of the clonal evolution of this patient's cancer ( Figure 5A). Canopy identified six genetically distinct clonal populations of tumor cells, present in various proportions within each tumor ( Figure 5B; Supplementary Table 5). These six clones are differentiated by ten unique groups of genomic alterations ( Figure 5A Table 2). Of note, group a contains truncal alterations ancestral to all tumor cells including two well-documented driver mutations, TP53 P278L and CDKN2A R80X, which were classified as ubiquitous ( Figure 4C). Interestingly, we also detected a truncal variant RPA2 V108L. RPA2 encodes a subunit of the heterotrimeric Replication Protein A complex that is critical for DNA replication, repair, recombination and DNA damage response [9].

DISCUSSION
Extensive multi-regional sequencing of tumors revealed the complex genetic landscape of treatmentrefractory cancers by demonstrating intratumor heterogeneity [10,11], which underscores the limitations of tumor profiling with tissue derived from a single biopsy specimen. From a clinical perspective, tumor heterogeneity contributes to the incomplete therapeutic responses seen in patients receiving different types of anticancer therapies including molecularly targeted therapies. Mechanistically, tumor heterogeneity drives acquired therapeutic resistance by facilitating the selection and expansion of therapy-resistant clones. Research autopsy has emerged as a powerful approach to characterize tumor heterogeneity at different metastatic sites in the individual patient. Together with information from treatment-naïve samples, sequencing of metastatic tumor samples from autopsy can provide a comprehensive molecular portrait of advanced cancers that reflects changes in tumor cell populations and genomes over time [12,13]. Here we performed research autopsy on a patient with widely metastatic IDCS. Genomic profiling of his tumors revealed somatic ultra-hypermutation and tumor heterogeneity in the form clonal diversity. Recently, Campbell et al. demonstrated a prevalence of ultra-hypermutation (TMB >100 mut/Mb) of 0.6% in a cohort of 78,452 adult cancers [4]. Therefore, the rare histologic classification of IDCS in this patient is accompanied by a rare 'genotype' of ultra-hypermutation.
Hypermutation in cancer may arise from intrinsic defects in DNA damage repair pathways or extrinsic mechanisms such as mutagenic exposure [4]. The majority of ultra-hypermutation detected by Campbell et al. was attributed to mutations in mismatch repair (MMR) genes and replication-associated DNA polymerases Polε or Polδ. While we did not detect genomic alterations involving MMR or Polε/Polδ to explain this patient's ultrahypermutated cancer, we did detect a truncal amplification of a region on chromosome 22q containing genes encoding the APOBEC3A-H (or A3 subfamily) cytidine deaminases. The deregulated expression and activity of A3A and A3B have been linked to cancer development through enhanced DNA hypermutation and unfaithful RNA editing [14][15][16]. Furthermore, we identified a truncal RPA2 variant with a predicted (by DUET [17]) destabilizing missense mutation in the C-terminal winged helix domain important for protein-protein interaction that could have further contributed to defective DNA damage repair and hypermutation. Together with driver mutations in tumor suppressors TP53 and CDKN2A, the above events could have produced ultra-hypermutation seen in this patient's cancer.
Different models of mutation accumulation have been proposed for hypermutated cancers and could be classified as 'steady' versus 'dynamic' hypermutation [4]. In the steady model, mutations accumulate gradually due to continuous mutagenic exposure or germline MMR deficiencies, leading to hypermutation over time. The dynamic model, seen in cancers with Polε or Polδ mutations, is characterized by 'bursts' of mutations followed by rise in genome-wide TMB [4]. Through analyses of clonality and mutational signatures, we hypothesize a gradual mode of hypermutation consistent with the former model in our patient's cancer. Genomic profiling of treatment-naïve tumor revealed high TMB at baseline that may have been the result of UV exposure, producing the driver mutations in TP53 and CDKN2A. TMB analysis of his autopsy tumor samples demonstrated an increase of ~20-37 additional mut/Mb that could be attributed to the truncal amplification of the A3 subfamily of cytidine deaminases, as indicated by the presence of APOBEC-specific Signature 2, and/or potentially therapy-induced mutations. Interestingly, the brain metastasis (T10) was the only autopsy tumor sample that had similar TMB as the treatment-naïve tumor (T1) while still retaining the APOBEC signature. A technical explanation underlying this finding may be due to lower TC contents in both samples affecting variant detection. Alternatively, this may reflect similar biology between the treatment-naïve tumor sample obtained from neck surgery and the brain metastasis; this similarity is demonstrated by the shorter genetic distance between 'normal' and T1 or T10 samples relative to other tumor samples in the NJ tree. Finally, the genetic dissimilarity between the brain metastasis and tumor samples from other organs also suggests the presence of unique genetic features of cancer cells with increased propensity for invasion of the central nervous system. Although high TMB has been established  as a clinically useful biomarker of response to checkpoint blockade in multiple human cancers [18,19], our patient developed progressive cancer while receiving adjuvant therapy with the PD-1 inhibitor nivolumab and had disease progression while being treated with dual CTLA-4/PD-1 blockade. Therefore, his case highlights the need for identification of additional biomarkers to predict clinical benefit for immunotherapy in patients with hypermutated cancers.
In summary, we present the first genomic characterization of metastatic IDCS, an extremely rare neoplasm of dendritic cell origin that lacks any standard therapy. This patient had metastatic IDCS characterized by ultra-hypermutation and clonal heterogeneity, likely through a combination of chronic mutagen exposure (UV), acquired defects in pathways important for DNA repair (TP53, CDKN2A, RPA mutations), and gain of genes that promote DNA hypermutation (APOBEC3A-H). His cancer was aggressive and refractory to multiple anticancer therapies including molecularly targeted agents and immunotherapies. In the future, it will be important to study additional patients with this and other rare cancers with hypermutation. Finally, the broader clinical implication of our results is that although patients with hypermutated cancer, originating from either somatic or germline genomic aberrations, are more likely to benefit from checkpoint inhibition, research is still needed to stratify these patients to maximize therapeutic efficacy and identify the different genetic determinants of primary or acquired resistance to immunotherapy.  Table 1. www.oncotarget.com

Research autopsy
The patient was consented to an IRB-approved clinical study for tumor profiling and body donation. At time of death, the patient's next-of-kin (and/or hospice agency) notified members of the research team.
The deceased was transported to the OSU Medical Center and a limited research autopsy was conducted within 3 hours after patient's passing. Metastatic tumors and adjacent normal tissues were sampled and immediately archived as fresh frozen specimens in OCT medium. After conclusion of the autopsy, the deceased was transported back to the designated funeral home within 24 hours.  Table 1.

Whole exome sequencing
We extracted genomic DNA from frozen tumors and normal (blood) samples and prepared sequencing libraries using an established protocol [20] that included enrichment with the xGEN ® Exome Research Panel v1.0 from Integrated DNA Technologies. Sequencing was performed on an Illumina HiSeq4000 at The Genomics Services Laboratory at Nationwide Children's Hospital (Columbus, Ohio) and achieved a median depth of 100-286x.

Alignment, variant calling, and annotation
Sequencing reads were aligned to hg19 using bwa [21] version 0.7.14. Deduplication was performed using Picard [22] version 2.3.0. Quality recalibration and local realignment around indels was performed with Picard and GATK [23] version 3.5. Somatic SNVs (single-nucleotide variations), somatic indels, and germline SNVs were called using VarScan2 [24] version 2.3.9. Somatic SNVs were filtered using bam-readcount as follows: minimum average base quality of variant-supporting reads 22, average variant distance to 3' read end of 0.24, Fisher's exact test P ≤ 0.05, maximum average sum of base quality mismatches 100, maximum average mismatches per base of variant-supporting reads 0.04, and minimum VF (variant fraction) 6%. Germline SNVs were filtered as above, except for Fisher's exact test P ≤ 0.1 (as recommended by VarScan2 documentation). Somatic indels were filtered with Fisher's exact test P ≤ 0.05 and VF ≥ 6%. In addition, a list of repetitive regions in hg19 was generated using the RepeatFinder tool included in MANTIS, modified to output all repeats of 1-mers to 5-mers spanning at least 4 bp. Indels falling within this list were excluded from further analysis.
Somatic SNVs and indels were annotated with ANNOVAR (revision #11f4bb, 2016-02-01) [25]. Mutational signatures from the COSMIC Mutational Signatures set [26] were called with deconstructSigs [27] version 1.8.0 with default settings and exome2genome trinucleotide frequency correction, run on R version 3.3.2. Microsatellite instability testing was performed using MANTIS version 1.0.3, run with three threads and otherwise default settings.
To compute a neighbor joining tree of per-sample phylogeny, we first generate a distance matrix D ∈ + ( ) × + ( ) For sample i ϵ1..N, define S i as the set of somatic SNVs in that sample called as above. Define S N+1 as the empty set, representing germline (with no somatic SNVs by definition). We compute D as follows: where △ is the set symmetric difference. Neighbor joining was performed over D using RapidNJ [28]

Copy number variation
Copy number variations (CNVs) were called using FALCON version 0.2 with threshold 0.3, and the QC procedure provided with Canopy was run with default length and ΔCN settings. The read depth ratio for each sample was computed as the ratio of aligned reads in tumor to aligned reads in normal. CNVs called by FALCON were manually inspected to identify events with major copy number >2 or minor copy number < 0.5 in at least one tumor sample, and spanned at least 25% of a chromosome. For each of these CNVs, common breakpoints were estimated across all tumors. FALCON was then re-run in each sample with τ chr set to the nearest SNP and threshold 0.2, with manual rescue of CNVs (threshold ≥ 0.1 in all cases).

Subclonal phylogenetic analysis
The reference read count, alternate allele count, and variant fractions of all nonsynonymous somatic variants called within each tumor sample were compiled. Somatic SNVs were further filtered with at least 100x coverage in all samples, minimum of 20 alt-supporting reads in at least one sample, not on a sex chromosome, and a DANN [30] predicted mutational impact score ≥ 0.96 (n = 893). This list of high-confidence mutations was downsampled to a set of 60 mutations with maximal diversity using a high correlation filter, implemented as follows: while |M| > n: i ← arg max x (ρ MxMy ); y > x; x, y ϵ 1.
where M is the set of mutations, with each mutation being a vector of its per-sample VFs, n is the desired final number of mutations (in this case, n = 60), and ρ is the correlation between VFs in each sample, defined as follows: where Cov(x, y) is the covariance between elements of same-length vectors x and y, and σ x is the population standard deviation of elements of vector x. Canopy was then run with an in-house parallelized version of sample.cluster mode with cluster number from 2 to 9, 10 MCMC clustering runs, τ K+1 0.05, number potential subclones from 3 to 9, 50 chains per subclone number, burnin 100, thinning parameter 5, simulation runs from 20000 to 50000, and writeskip 200. As FALCON does not generate standard deviations for regions with identical major and minor copy number, and Canopy does not support sparse ԑ M or ԑ m matrices, the Frobenius norm of non-NA values of the ԑ M and ԑ m matrices was used for Canopy. The optimal subclone number was selected by highest BIC (Bayesian Information Criterion) score. www.oncotarget.com Afterwards, the remaining somatic SNVs (not used for Canopy tree building) and indels were retroactively assigned to the resulting tree using a maximal likelihood model. For an arbitrary somatic mutation v, define N as the number of tumor samples,  r N ∈  as the number of alternate-supporting reads for v in each sample,  x N ∈  as the coverage of the variant position in each sample, and  p N ∈  , a vector of expected variant fractions in each sample (derived from the tree as shown below). Given known  x and  p , and modeling alternative read depth as jointly binomial across tumor samples, the probability and log-likelihood can be computed as follows: To compute  p, we utilize the phylogenic tree and clonal fractions reported by Canopy, along with its persubclone major and minor copy number estimates. Define Z as the number of edges of the tree, and K as the number of tumor subclones + 1 (to account for normal cells). Given the tree, a matrix T ∈{ } × 0 1 , Z K representing the phylogenetic composition of each subclone can be constructed as follows (note that the tree root corresponds to germline): Next, define P ∈ ×  K N as the clonal fraction matrix provided by Canopy (where entries in the first column represent proportion of normal cells in each sample), T as the number of CNV regions used for Canopy,  C M ∈ ×  ( )

T K
as the per-subclone major copy number matrix, and as the per-subclone minor copy number matrix. Any zero values within each column of P were set to 0.0005, with the remaining entries reduced equally to ensure that each column sums to 1. This is first used to compute A ∈ ×  ( ) , T N the total copy number of each CNV region in each sample (independent of any specific mutation), as follows: Next, assume without loss of generality that v is within edge z. From Canopy, if v is within a CNV t, t is within edge y. We can calculate the copy number of v in each subclone,  b K ∈  , as follows: To perform the assignment, ( ) v z ∈ is computed as described above for all edges. If v is within a CNV region and the candidate edge is before the CNV, its possibilities of being on the major or minor allele are both considered. If the candidate edge is after the CNV, only the possibility of it occurring on one copy of the allele is considered, utilizing the simplifying assumptions (also made by Canopy) that each somatic mutation occurs only once and no back-mutations occur. If the candidate edge contains the CNV, all three above possibilities (major/minor/after) must be considered. To avoid numerical errors, if p i = 0 and r i ≠ 0 for a candidate assignment in any sample, that candidate assignment is discarded. The edge with highest log-likelihood is selected, with ties (such as can occur with deletions or on a different branch than a CNV) broken in favor of the edge furthest from the root node. This approach was implemented and run using Python 2.7.1, using scipy [31] 0.10.1.