Abstract Neuroblastoma (NB) is a childhood solid malignant tumor originating from precursor cells of the peripheral nervous system. We have previously established a risk classification system based on DNA copy number profiles. To further explore the pathogenesis of NBs in distinct risk groups, we performed whole-exome sequencing analysis of 57 primary and 7 recurrent/metastatic tumors with unique chromosomal aberration profiles as categorized by our genomic sub-grouping system. Overall, a low frequency of somatic mutations was found. Besides ALK (4/64, 6.3%), SEMA6C, SLIT1 and NRAS, genes involved in the axon guidance pathway, were identified as recurrently mutated in 6 of 64 tumors (9.4%). Pathway enrichment analysis revealed enrichment of 25 mutated genes in the mitogen-activated protein kinase (MAPK) pathway, 13 genes in the Wnt pathway, and 12 genes in the axon guidance pathway. Genomic analyses demonstrated that primary and matched recurrent or metastatic tumors obtained from sporadic and monozygotic twin NBs were clonally related with variable extents of genetic heterogeneity. Monozygotic twin NBs displayed different evolutionary trajectories. These results indicate the involvement of the axon guidance, MAPK and Wnt pathways in NB and demonstrate genomic diversity with NB progression. Keywords: whole-exome sequencing, array CGH, neuroblastoma, familial neuroblastoma, genomic diversity INTRODUCTION Neuroblastoma (NB) is one of the most common solid tumors in children, which originates from precursor cells of the peripheral (sympathetic) nervous system. One striking characteristic of NB is the heterogeneous biological and clinical behaviors exhibited by individual tumors. Patients over 1 year of age usually suffer from disseminated disease at diagnosis and have a poor outcome despite intensive multimodal treatment [[62]1–[63]3]. Currently, NB patients are stratified into low-, intermediate-, or high-risk categories based on a risk assessment of well-defined prognostic factors including the age at diagnosis, International Neuroblastoma Staging System (INSS) stage, tumor histopathology, MYCN amplification status and tumor DNA ploidy. The 5-year survival rates (SRs) for NB are approximately 95% for the low-risk group and 80–90% for the intermediate-risk group, but only 30–50% for high-risk NB. The extreme clinical heterogeneity of NB reflects the complexity of genetic and genomic events associated with disease development and progression. To date, several genomic changes that correlate with an unfavorable prognosis for NB have been identified. MYCN oncogene amplification is present in approximately 20% of NB patients and correlates with tumor progression, advanced disease stages and poor outcome [[64]1, [65]2]. Additionally, chromosomal aberrations such as 1p loss, 11q loss and 17q gain have been shown to predict poor patient outcome [[66]1, [67]2]. Several genetic mutations associated with NB have also been identified. Germ-line mutations in paired-like homeobox 2B (PHOX2B) have been discovered in a minority of familial NB patients and account for their predisposition to NB [[68]4–[69]8]. Activating mutations in the tyrosine kinase domain of the anaplastic lymphoma kinase (ALK) oncogene are acquired hereditarily or somatically in familial and sporadic NB cases, respectively, and account for disease susceptibility as well as disease progression [[70]9–[71]12]. Recent studies using next-generation sequencing (NGS) approaches have identified mutations in the ATRX gene that are strongly associated with late onset age of NB [[72]13, [73]14], a high frequency of chromothripsis in stage 3 and 4 tumors (18%) and frequent mutations in the Rac/Rho pathway genes that guide neuritogenesis [[74]15, [75]16], and recurrent TERT rearrangements in high-risk NBs (24%) [[76]17, [77]18]. Chromosomal deletions and sequence alterations in the ARID1A/1B genes, which regulate chromatin remodeling, have also been detected in 11% of NB patients [[78]19]. We have previously established a genomic subgrouping system based on array CGH analysis for the risk stratification of neuroblastoma [[79]20, [80]21]. To further identify novel somatic mutations linked to tumorigenesis and different outcomes in NB, in this study we focused on the two major subgroups, Ss and P1a, with favorable and unfavorable clinical outcomes, respectively, to perform whole-exome sequencing (WES) analysis using the Illumina platform. Our results indicate that the axon guidance, MAPK and Wnt pathways may be involved in molecular behavior of NB. Furthermore, we describe the genetic heterogeneity within NB tumors and evolutionary trajectories of monozygotic twin NBs. RESULTS Genomic subgrouping and risk stratification of NBs We have previously established a genomic sub-grouping system based on array CGH analysis for the risk stratification of neuroblastoma in 236 primary NBs [[81]18, [82]19]. Using this system, we have further investigated an additional independent cohort of 107 primary NBs for their risk classification. Array CGH analysis revealed three major chromosomal aberration profiles in the 343 total NBs examined: a silent (S) pattern almost without any chromosomal aberrations; a partial gains and/or losses (P) pattern; and a whole gains and/or losses (W) pattern. Each group was further classified into several subgroups on the basis of clinical outcome and known genomic signatures, including MYCN amplification, 1p loss, 11q loss and 17q gain. Consequently, the S group was divided into Ss (s, MYCN non-amplification) and Sa (a, MYCN amplification) subgroups according to the MYCN copy number, and the P group into P1 (1p loss and 17q gain), P2 (1p loss, 11q loss and 17q gain), P3 (11q loss and 17q gain) and P4 (17q gain) subgroups, each of which was further split into s and a subsets in light of MYCN amplification status. The W group was characterized by whole chromosomal gains and/or losses, especially the gain of chromosome 17, and rarely exhibited MYCN amplification. Among them, the major subgroups with high frequency are Ss (n = 67, 20%), P1a (n = 58, 17%) and W4s (n = 87, 25%). Notably, the Ss subgroup had an overall 8-year SR of 82%, whereas the P1a subgroup, which possessed genomic signatures predicting unfavorable outcome in NB including 1p loss, 17q gain and MYCN amplification, had an overall 8-year SR of 33% (Figure [83]1). Moreover, P1a tumors harbored a high frequency of aberrations in the ALK gene (9/58, 15.5%; including mutation (n = 6) and amplification (n = 3)). Figure 1. Genomic profiles of Ss and P1a subgroups and their association with overall survival in NB. Figure 1 [84]Open in a new tab (A) Genomic aberration profiles of Ss and P1a subgroups identified by array CGH analysis. The 8-year survival rate was calculated using a cohort of 343 NBs. (B) Kaplan–Meier survival curves drawn for Ss versus P1a subgroups. Survival distributions were compared using the log-rank test. Whole exome sequencing analysis To identify novel somatic mutations linked with tumorigenesis and different outcome in NB, in this study we employed a WES approach to 56 paired NBs (57 primary tumors and 7 recurrent/metastatic tumors; median age, 19 months; [85]Supplementary Table 1) with a focus on Ss (n = 21) and P1a (n = 17) subgroups. All samples were selected on the basis of the availability of high-quality and sufficient DNA for WES. Our sample set contained 51 sporadic and 5 familial background NBs. Additionally, 7 recurrent/metastatic tumors from individual cases were examined so as to investigate genome heterogeneity between primary and recurrent/metastatic tumors. In total, 64 tumors and their matched blood-derived DNA samples were examined by WES using an Illumina Hiseq2000 platform with a minimum sequencing depth of 75X per case. We identified 1685 candidate somatic mutations in 1263 genes from 64 tumors, including 417 silent mutations, 1132 missense, 60 nonsense, 19 splice-site changes and 57 insertions or deletions (indels). On average, each tumor harbored 26.3 coding mutations (19.8 non-silent mutations comprising 18.9 substitutions and 0.9 indels) with a wide range of 3–204 (Table [86]1 and [87]Supplementary Table 2). Mutation spectrum analysis showed that single-nucleotide mutations were predominantly C–T substitutions (Figure [88]2A), similar to previous observations for NB [[89]22] and adult tumors [[90]23]. Moreover, the mutation rate is strongly correlated with MYCN amplification (P = 1.10×10^−4; Figure [91]2B). Table 1. Summary of somatic mutations identified in 64 NB tumors. Mutation type Number Description Silent 417 Mutation causing no change in the amino acid Missense 1132 Mutation encoding a different amino acid Nonsense 60 Mutation encoding a premature stop codon Splice-site change 19 Mutation within 2bp from splicing site Indel 57 Small insertion or deletion [92]Open in a new tab Figure 2. NB somatic mutation characteristics. Figure 2 [93]Open in a new tab (A) The mutation spectra of somatic single nucleotide mutations. (B) Average number of somatic alterations in NB patients with MYCN amplification (MYCN Amp) and MYCN non-amplification (MYCN Sing). (C) The number of somatic mutations in Ss, P1a and recurrent/metastatic tumors. (D) Average number of somatic alterations in Ss and P1a subgroups. The number of patients per group is shown in parentheses. Box plot in R software was used to create the graph. Permutation test in R software was used to assess statistical difference. The recurrent/metastatic tumors harbored the most somatic mutations, as compared with Ss and P1a primary tumors (Figure [94]2C). P1a tumors possessed significantly more somatic mutations than Ss tumors (permutation test, P = 6.35×10^−6; Figure [95]2C and 2D). Of the 1268 non-silent somatic mutations, only 211 were discovered in 21 primary and 2 recurrent Ss tumors (mean, 9 mutations per tumor; range, 2–18). In contrast, 520 mutations were detected in 17 primary and 2 recurrent/metastatic P1a tumors (mean, 27 mutations per tumor; range, 2–160), whereas recurrent and metastatic tumors harbored the highest number of somatic mutations (mean, 57 mutations per tumor; range, 0–174). Notably, four recurrent/metastatic tumors show fewer mutations that are comparable to that of the Ss primary tumors. Of these, two tumors (NB56-T2 and NB57-T2) were obtained from the metastatic lesions along with the primary tumors from the twin patients with stage 4s NB. Their genomic profiles are both W4s characterized by whole chr17 gain without 1p loss, 11q loss and MYCN amplification. Considering that the overall 5-year SR of the W4s subgroup in our dataset is 97% [[96]20, [97]21] and that the two patients are alive, such a lower number of somatic mutations even in the metastatic lesions account for their favorable outcomes. The other two tumors (NB37-T2 and NB51-T2) are recurrent ones that were extracted after neoadjuvant therapy. In view that there are significant differences in genomic profiles between the primary and recurrent tumors (P1a in NB37-T1 and Ss in NB37-T2, Sa in NB51-T1 and Ss in NB51-T2), fewer mutations in the two recurrent tumors may be attributed to the loss of viable tumor cell content due to neoadjuvant therapy. Potential genes and pathways We detected 27 genes with non-silent somatic mutations occurring in at least two NB cases (Table [98]2). The mRNA expression levels of these genes and their association with overall survival in NB were verified using our expression data and two independent NB microarray analysis datasets (Kocak - 649 - custom - ag44kcwolf and SEQC - 498 - RPM - seqcnb1) in the R2: microarray analysis and visualization platform ([99]http://r2.amc.nl), except OR52R1 whose data were absent in one of the datasets. Of these genes, ALK mutations, which were all restricted to the kinase domain, were found in 7.0% (4/57) of the primary tumors, consistent with the frequencies identified in several large NB cohorts using conventional sequencing method [[100]9–[101]12]. Apart from ALK, mutations in TTN and EMR1 were detected in 5 and 3 out of 64 NB tumors (7.8% and 4.7%), respectively. PCDHGA4, PLCH2, NRAS and SEMA6C harbored mutations at the same sites in two independent cases. Gain-of-function mutations to NRAS have been described in many adult tumors as well as stage 4 NBs recently [[102]14]. In our dataset, the p.Gln61Lys mutation was identified in two P1a tumors: one primary tumor from patient NB43 and the other a recurrent tumor from NB53 (-T2). Both individuals died of disease. Notably, this mutation was also identified in 3 of 25 NB cell lines. Somatic mutations in PLCH2 (p.Thr177Met), PCDHGA4 (p.Asp536Tyr) and SEMA6C (p.Leu348Arg) commonly appeared in NB56-T1/T2 and NB57 (patients NB56 and NB57 are twins) at stage 4s with the W4s genomic profile. None of the 27 genes were commonly mutated in two Ss tumors. However, SLIT1, SYNE1 and VPS13C were mutated in two primary P1a tumors, respectively. Table 2. Summary of non-silent mutations occurring in ≥ 2 cases across 64 NB tumors. Gene symbol Mutation site Number of cases mRNA expression* TTN 7 5 H>L MUC16 5 5 HL RAD54L2 2 2 HL VPS13C 2 2 H>L UTRN 2 2 H>L TRIM42 2 2 H=L SYNE1 2 2 H>L SLIT1 2 2 HL ENTHD1 2 2 H=L CYP4F11 2 2 H=L CMYA5 2 2 H=L LILRA3 2 2 H=L NPAP1 2 2 H=L POLN 2 2 H=L PCDHGA4 1 2 H>L PLCH2 1 2 H>L SEMA6C 1 2 H>L ZNF304 1 2 H=L NRAS 1 2 H=L [103]Open in a new tab *indicates the correlation of mRNA expression level with overall survival as determined by Kaplan-Meier analysis using two independent microarray analysis datasets from the open-access R2 interface (Kocak - 649 - custom - ag44kcwolf and SEQC - 498 - RPM - seqcnb1). For PCDHGA4, the correlation was evaluated using Versteeg - 88 - MAS5.0 - u133p2 dataset instead of Kocak - 649 - custom - ag44kcwolf dataset, since this gene does not exist in the Kocak - 649 dataset. H, high expression; L, low expression; H=L, there is no significant difference in overall survival between low and high expression subsets; H>L, high expression is significantly associated with good overall survival; H