Whole-genome sequencing of recurrent neuroblastoma reveals somatic mutations that affect key players in cancer progression and telomere maintenance

Neuroblastoma is the most common and deadly childhood tumor. Relapsed or refractory neuroblastoma has a very poor prognosis despite recent treatment advances. To investigate genomic alterations associated with relapse and therapy resistance, whole-genome sequencing was performed on diagnostic and relapsed lesions together with constitutional DNA from seven children. Sequencing of relapsed tumors indicates somatic alterations in diverse genes, including those involved in RAS-MAPK signaling, promoting cell cycle progression or function in telomere maintenance and immortalization. Among recurrent alterations, CCND1-gain, TERT-rearrangements, and point mutations in POLR2A, CDK5RAP, and MUC16 were shown in ≥ 2 individuals. Our cohort contained examples of converging genomic alterations in primary-relapse tumor pairs, indicating dependencies related to specific genetic lesions. We also detected rare genetic germline variants in DNA repair genes (e.g., BARD1, BRCA2, CHEK2, and WRN) that might cooperate with somatically acquired variants in these patients with highly aggressive recurrent neuroblastoma. Our data indicate the importance of monitoring recurrent neuroblastoma through sequential genomic characterization and that new therapeutic approaches combining the targeting of MAPK signaling, cell cycle progression, and telomere activity are required for this challenging patient group.

Constitutional alterations in cancer-predisposing genes. To detect possible genetic predisposition to NB, constitutional DNA was analyzed, focusing only on variants in 565 genes with established associations with sporadic and/or hereditary cancer (Supplemental Table 2). This analysis did not detect any activating point mutations in hot-spot regions of ALK or loss-of-function variants in PHOX2B, the two genes mainly connected to NB predisposition. However, a novel ALK frameshift variant (i.e., ALK p.V1471fs*6) of unknown significance was detected in one patient. This variant is predicted to result in a truncated protein that retains the kinase domain but lacks a portion of the C-terminal domain containing multiple phosphorylation sites crucial for ALK functions including tyrosine residue 1604. Thus, it is unlikely that this variant has the transforming potential seen for of the ALK hot-spot mutations common in NB. Three of the seven patients had germline variants likely to affect the function of DNA repair genes. These include a nonsense mutation in BARD1 (p.Q545), a missense variant in WRN (p.S1141L) affecting an ATM phosphorylation site important for WRN protein function, and two rare missense variants in CHEK2 (p.I157T) and BRCA2 (p.R2991C), respectively, detected in one patient. These variants and other variants associated with loss-of-function or predicted to be damaging by SIFT and PolyPhen2 are summarized in Table 2.
Overall comparison of somatic variants between time of diagnosis and time of recurrence. In total, 529 nonsynonymous somatic SNVs were detected in the analyzed tumors (Table 3 and Supplemental  Table 3), with an average of 28 SNVs (range 8-63) in primary tumors, increasing to an average of 61 SNVs in relapsed tumors (range 51-76), corresponding to 0.7 and 2.1 variants per megabase of the coding genome, respectively. The mutational burden at time of relapse/progression was significantly higher than at time of diagnosis (p = 0.002) (Fig. 1A,B). Overall, the mutational load was in the range of that previously reported in NB. The median number of somatic SNVs shared by each primary and relapse/refractory pair was 14 (range 0-40), corresponding to 50% of all mutations detected in primary tumors and to 23% of all mutations detected in relapsed/ refractory tumors. However, for patient NB60R6, for whom the investigated biopsy specimens were from a metastasis at time of diagnosis and a metastasis at time of relapse, no somatic nonsynonymous SNV shared by the two samples could be detected, not even after including mutations below the 10% variant allele fraction.
The number of somatic mutations at time of diagnosis was correlated with patient age at time of sampling, at time of both diagnosis and relapse/progression. Pearson correlations were calculated after removing case NB67R9, which appears to have had a markedly higher mutation rate than the other samples (Fig. 1C). For the remaining tumor specimens, Pearson correlation coefficients indicate statistically significant correlation between the total number of non-synonymous SNVs and patient age at time of diagnosis or relapse occurrence (R = 0.915, p = 3.0 × 10 -5 ) (Fig. 1C).

Copy number alterations and structural variants.
Although some persistent copy number alterations (CNAs) were detected in matched primary and relapsed tumors, the relapsed tumors had increased numbers of structural variants (SVs) at time of relapse (Fig. 1B, Supplemental Fig. 1, and Supplemental Table 4). At time of diagnosis, a median of 38 somatic SVs (range 5-122) were detected, the number later increasing to a median of 49 SVs (range 14-89) at time of relapse when excluding the excessive number of SVs within highly amplified regions (i.e., the MYCN locus in NB63R4 and the 12q-amplicons in NB59R9). Persistent CNAs present in both diagnostic and relapsed/refractory tumors include: 17q-gain, detected in all tumors except one that instead displayed gain of whole chromosome 17; 11q-deletion, detected in four of six patients; MYCN amplification, detected in one patient; and CDK4/MDM2/FRS2 amplification, detected in one patient ( Fig. 2 and Supplemental Fig. 1). Loss of chromosome 1p, also associated with poor prognosis in NB 21 , was seen in six patients: three cases (i.e., NB67R9, NB63R4, and NB69R8) in which 1p-deletion was detected in both primary and relapse tumors, one case (i.e., NB60R6) in which 1p-deletion was detected only in the diagnostic specimen, one case (i.e., NB53R2) in which 1p-deletion was detected only in the relapsed specimen, and one case (NB67R5) in which 1p-deletion was detected in both primary and relapsed tumors, albeit with different breakpoints (Fig. 2 and Supplemental Table 4). On the genomic level, additional recurrent CNAs and structural variants were seen in association with TERT (in three patients, one relapse specific), CCND1 (in two patients, one relapse specific), and LSAMP (in two patients, both relapse specific) (Fig. 1D,E) www.nature.com/scientificreports/ the difference in translocation partners and breakpoints relating to these genomic alterations (Supplemental Table 4). In patient NB60R6 we also detected CNAs in both primary and relapse tumors, resulting in gain of www.nature.com/scientificreports/ chromosomal regions 7q and 6p; however, this also emerged through parallel events, as the two tumors have different breakpoints for these CNAs and also gains of different parental chromosomes as judged by the allele fraction (Figs. 1D and 2). The paired primary and relapse tumors of patient NB60R6 shared only one common SV, a translocation between chromosomes 11 and 17 that ultimately resulted in 11q-deletion and 17q-gain. However, as the position of the translocation was identical in both samples, this indicates that both analyzed specimens share a common descent. Patient NB59R9 had the highest number of SVs due to high genomic instability, with a cataclysmic event resulting in extensive focal genomic shattering and rejoining. This affected regions in chromosome 5 and amplified regions in chromosome 12, connected to MDM2 and CDK4, were present in both samples for this patient, while similar focal genomic shattering was seen in chromosomes 7, 18, and 21 in primary tumors and in chromosome 4 in relapse tumors (Supplemental Table 4).

Recurrent alterations and mutations in cancer-associated pathways.
Telomere maintenance is crucial for sustained proliferative capacity and is associated with poor prognosis in NB. In our cohort, tumors from three patients had TERT-associated SVs and tumors from two patients had tandem duplications causing disruption of ATRX (Fig. 1D,E). All three patients with TERT alterations had genomic profiles that included 11q deletions, though in two patients we did not detect any alterations in genes directly connected to telomere maintenance. Both these patients displayed high levels of genomic amplification: patient NB59R9 for several 12q loci, including CKD4, MDM2, and FRS2, and patient NB63R4 for the MYCN locus. Interestingly, patient NB67R5, with two different TERT-associated structural variants in primary and relapsed tumors, also had two mutations in the chromatin-associated genes ARID1A and ARID1B. The primary tumor of NB67R5 had a structural variant causing disruption of ARID1B, a variant that was lost in the relapse sample, which instead carried a missense mutation in ARID1A (Fig. 1D, and Supplemental Tables 2 and 3). Previous studies by us and others have indicated that mutations in ALK and other genes of the RAS-MAPK pathway are connected to relapsed NB 18,19 . In our cohort, an activating somatic ALK mutation was seen in the tumors from one patient (i.e., NB69R8). This is an ALK Q1275R alteration in which the tumor displayed increased variant allele frequency from time of diagnosis to time of relapse. This increase was likely due to a new copy-neutral loss of heterozygosity (CN-LOH) of chromosome 2 that preserves the mutated allele, as judged by B-allele frequency. In addition to ALK, somatic variants were also detected in other genes associated with RAS-MAPK signaling. This includes mutations in KRAS, FGFR1, and MAP3K14 (NIK), the first two genes being detected only in the relapsed samples of two different patients, while the last was shared by the primary and relapse tumors in one patient (Fig. 1D). We also detected a somatic frameshift mutation in PHOX2B, a gene important for proper neural crest development and present in both the primary and recurrent tumors of patient NB67R5. Somatic alterations were seen in several genes connected to cell cycle progression and/or DNA repair, such as CCND1, CDK4, MDM2, BRCA1, BRINP1, and WRN (Fig. 1D).   www.nature.com/scientificreports/ Besides CCND1 and TERT, 30 additional genes were altered in at least two tumor samples in our cohort, either through SNV or SV. Recurrent mutated genes present only in tumors at time of diagnosis were limited to ABCA5, CNTNAP2, CTNND2, and PTPRM, while we detected recurrent mutations in relapsed samples of 16 genes either already present in primary tumors or acquired in relapse. Of these 16 genes, eight were found to be relapse specific (i.e., CDK5RAP1, CFAP44, LSAMP, MUC16, NUBPL, POLR2A, PRR12, and SLC6A18) (Fig. 1D and Table 3). Two tumors carried missense mutations predicted to be damaging by SIFT and PolyPhen2 in the tumor-suppressor gene FAT1 (Fig. 1D). The FAT1 mutation in NB63R4 was acquired during relapse, while the FAT1 mutation in NB67R5 was present at the subclonal level with a 2% variant allele fraction in the primary tumor but increasing in relapse (i.e., a 15% variant allele fraction), indicating clonal expansion (Supplemental Table 3 and Supplemental Fig. 2).

Pathway-centered gene set enrichment of relapsed and primary neuroblastoma. Pathway
analyses were performed on somatically mutated genes detected among all relapsed and primary tumors, respectively, to detect differences in pathway signature. Gene set analyses in KEGG for all mutations present in relapse tumors identified pathways connected to the glutamatergic synapse, PI3K-Akt signaling, longevity regulation, and Rap1 signaling among the highest ranked enrichment although not at a statistically significant level after correction for multiple testing (Supplemental Fig. 3 and Supplemental Table 5). Gene set analyses also indicated the prior enrichment of genes involved in PI3K-Akt signaling and in cancer pathways in the diagnostic samples (Supplemental Fig. 3 and Supplemental Table 5).

Discussion
In attempting to identify molecular features related to disease progression in NB, we have characterized the genomes of seven relapsed/progressive NB tumors and corresponding tumor specimens from the time of diagnosis. In this study, we show that relapsed tumors generally have a higher mutational load and have acquired additional segmental variants (Fig. 1A,B, and Supplemental Fig. 1), in line with previous studies of relapsed NB 19,20 . We also show that the number of mutations present in primary tumors and the number of acquired mutations in relapsed tumors correlates well with patient age at diagnosis and recurrence, respectively (Fig. 1C). The coverage depth used for WGS in this study limits the ability to analyze whether relapse-specific mutations were acquired de novo or emerged from subclonal events in the primary tumor. However, we do detect examples of the clonal expansion of specific mutations, for example, as seen for FAT1 (Supplemental Fig. 2), as well as the presence of clonal eradication, as only a mean of 36% of all mutations detected in primary tumors are also present in relapse tumors. For some samples, such as NB60R6 and NB67R9, there is striking in-patient heterogeneity, with very few somatic SNVs shared by different tumor samples (Fig. 1E). This indicates that the longitudinal At diagnosis Common At relapse www.nature.com/scientificreports/ www.nature.com/scientificreports/ www.nature.com/scientificreports/ monitoring of disease burden and treatment response using specific genetic alterations could be challenging for a subgroup of NB patients. Recurrent alterations were seen in 32 genes in tumors from at least two patients; eight of these genes were relapse specific and include MUC16 and LSAMP (Fig. 1D). Tumors from three patients also display alterations connected to cell cycle genes. These include tumors from two patients with CCND1 gain and tumors from one patient with CDK4/MDM2 amplification (Fig. 2).
Mutations in ALK have previously been shown to be enriched in relapsed NB 18 , as have other mutations that activate RAS-MAPK signaling, such as RAS, NF1, BRAF, and FGFR1 19 . Tumors from four of our seven patients carry somatic ALK-MAPK-RAS signaling pathway alterations. These include somatic mutations in ALK and MAP3K14 that were already present at time of diagnosis and mutations in KRAS and FGFR1 that are relapse specific (Fig. 1D). The KRAS K117N alteration, detected in patient NB67R5, has been shown to decrease GTPase activity 22,23 although studies from colorectal cancer indicate this alteration to be less potent than mutations in codons 12 and 13 22 . The tumor with KRAS K117N also carried a FAT1 mutation present at the subclonal level in the primary tumor and later found to be enriched at time of relapse. FAT1 is a tumor-suppressor gene that restrains the oncogenic activity of the Yes-associated protein (YAP) 24 , and studies of other malignancies show that YAP is a critical effector of oncogenic RAS 25,26 . Therefore, it is possible that the malignant capacity of KRAS K117N is further potentiated through increased YAP signaling in relapsed tumors.
Both PHOX2B and ALK have previously displayed recurrent mutations in both familial and sporadic NB 4,5,27 , as they are crucial for proper development of the sympathetic nervous system. In our cohort, the tumors from one patient (i.e., NB67R5) shared a somatic frameshift mutation of PHOX2B (p.G230fs). Truncating PHOX2B mutations have been shown to possess a dominant-negative effect resulting in halted sympathetic nervous system differentiation and increased proliferation, making cells susceptible to additional transforming events 28 . Regarding PHOX2B and ALK alterations in constitutional DNA, we detected a novel ALK frameshift variant already present in the germline (p.V1471fs*6) in one investigated patient. Although ALK is of great interest in NB etiology, it is difficult to envision the role of this particular variant in the oncogenic process, as it is expected to lack phosphorylation sites essential for ALK activation.
However, it has been shown that children and adolescents with cancer have a higher incidence of pathogenic or likely pathogenic germline mutations than do healthy controls 29 . When expanding our investigation of other cancer predisposition genes besides ALK and PHOX2B, we detect rare variants such as a nonsense mutation in BARD1 as well as missense variants in CHEK2, BRCA2, and WRN (Table 2). Germline variants in BARD1 and CHEK2 have indeed been shown to be enriched in NB patients 16,30 , while BRCA2 has been shown to be enriched in pediatric cancer patients 29 . The gene WRN encodes a RecQ helicase important for maintaining genome stability through resolving stalled replication forks in association with double-strand breaks. The missense variant detected in WRN affects an ATM phosphorylation site required for correctly handling stalled forks 31 . The tumors of patients who are heterozygous for the WRN variant display high genomic instability, with extensive focal genomic shattering and rejoining that might be attributed to deficient double-strand break repair.
For one patient (i.e., NB60R6), we analyzed a metastasis at time of diagnosis and a metastasis at time of relapse that provided indications of branched evolution. These two tumor specimens displayed no common nonsynonymous SNVs and only one shared structural variation, a translocation between chromosomes 11 and 17 resulting in 11q-deletion and 17q-gain (Fig. 3A). That the fusion point was exactly the same in both samples indicates that they emerged from a common ancestral tumor cell progenitor and not through multifocal events. However, both lesions acquired 6p-gain, 7q-gain, and CCND1-gain, albeit with different breakpoints and gains of different parental alleles (Figs. 2 and 3B). Additional parallel events between primary and relapse tumors were seen in case NB67R5, which displayed two different SVs close to the TERT locus, and in case NB67R9, which displayed two different mutations in ARID1B and ARID1A, two members of the SWI/SNF chromatin-remodeling complex, in tumors at time of diagnosis and time of relapse, respectively (Fig. 1D). The primary and relapsed tumors from patient NB67R9 also had a common TERT rearrangement that might be further potentiated through the ARID1B and ARID1A mutations, as loss of ARID1A has been shown to enhance TERT transcription 32 .
The findings of alterations of TERT, ARID1A, and ARID1B are related to telomere maintenance, which is a strong predictor of poor outcome in NB 13,33 . Besides the MYCN-amplified NB case, in which sustained telomere length might be achieved through elevated TERT expression driven by MYCN, all six non-MYCN-amplified cases in our cohort displayed genomic aberrations associated with telomere maintenance or chromatin remodeling; three of these cases had TERT rearrangements, two had ATRX alterations, and one had a mutation in ARID2, which is also part of the SWI/SNF chromatin-remodeling complex. This further supports the finding that telomere maintenance is associated with an unfavorable clinical course in NB 13,33 , as most tumors in this study were characterized by different means of sustaining telomeres.
In this study we have identified a diverse set of gene mutations in tumors from time of diagnosis and time of recurrence in this challenging patient group. This set includes mutations in genes promoting cell cycle progression, activating ALK-RAS-MAPK signaling, or important for telomere maintenance and immortalization. Our findings in combination with previous studies of relapsed NB provide a rationale for the comprehensive genomic characterization of recurrent NB in order to optimize treatment and aim for prolonged survival. This information could be used to spur novel combination therapies using compounds abrogating MAPK signaling, cell cycle progression, and telomere activity.

Material and methods
Patient material and ethics. NB  Only high-quality called variants with a minimum 10% variant allele frequency and a total read coverage of ten were considered for further analysis. All synonymous variants or variants in non-coding regions were excluded, except those affecting canonical splice sites and remaining variants were assessed manually using the Integrative Genomics Viewer (IGV) 34 to remove calls due to mapping artifacts and paralogs. The germline variants were analyzed using a more stringent filtering approach excluding all variants with a population allele frequency above 1.0%, in either gnomAD v2.  Table 2. The filtering of somatic and germline variants was done using the ingenuity variant analysis software v5.6 (Qiagen, Hilden, Germany). The Canvas tool (version 1.38.0.1554) 35 was used to call copy number alterations (CNAs), with gains and losses of genomic regions as well as LOH regions being identified using read depth coverage in combination with SNV b-allele frequencies. For the calling of germline and somatic CNVs, the SmallPedigree-WGS and Somatic-WGS algorithms of Canvas were used, respectively. Somatic structural variants (SVs) were called using the Manta tool (version 1.1.1) 36 , which applies information about soft-clipped read ends, disjointed read pairs, and read-pair orientation to identify larger structural variations (e.g., deletions, duplications, inversions, and translocations). Calls from the constitutional DNA were used to filter out germline variation and artifacts caused by problematic regions. In addition, further filtering and removal were done based on presence in the SweGen Variant Frequency dataset (https ://swefr eq.nbis.se/) or in our in-house set of normal controls. Genes affected either directly by SVs or through smaller focal gain/loss were also noted for further genomic analysis.
To verify same individual origin for matched tumor-normal pair, we ran a Python script developed in-house (available on request) that calculated the fraction of shared single nucleotide polymorphisms (SNPs) relative the total number of SNPs using 400,000 SNPs (all present in dbSNP; version 138) from chromosomes 1-20. Using this method, the fraction reliably ends up in the 0.97-1 range for verified tumor-normal matched pairs, whereas unrelated randomly selected samples have a fraction between 0.6 and 0.7, while a parent and child-level of relatedness ends up close to 0.8. All paired samples in this study showed a tumor-normal matched fraction above 0.97 and thus, authenticated sample identity.
Statistics. IBM SPSS software was used to calculate Pearson correlations between mutational loads at diagnosis and relapse, respectively, relative to age at time of sampling and to perform independent t-testing for differences in the number of somatic SNVs between primary and relapsed tumors. Gene set analysis for the enrichment of somatically mutated genes was conducted using the curated Kyoto Encyclopedia of Genes and Genomes (KEGG) cell signaling pathway database 37 (www.kegg.jp/kegg/kegg1 .html) through the Enrichr web server 38 .