Genomic pro�le of MYCN non-amplied neuroblastoma and potential for immunotherapeutic strategies in neuroblastoma

Background: MYCN ampli�cation is the most important genomic feature in neuroblastoma (NB). However, limited studies have been conducted on the MYCN non-amplied NB including low-and intermediate-risk NB. Here, the genomic characteristics of MYCN non-amplied NB were studied to allow for the identi�cation of biomarkers for molecular strati�cation. Methods: Fifty-eight whole exome sequencing (WES) and 48 whole transcriptome sequencing (WTS) samples of MYCN non-amplied NB were analysed. Forty-one patients harboured WES and WTS pairs . Results: In the MYCN non-amplied NB WES data, maximum recurrent mutations were found in MUC4 (26%), followed by RBMXL3 (19%), ALB (17%), and MUC16 and SEPD8 (14% each). Two gene fusions, CCDC32 - CBX3 (10%) and SAMD5 - SASH1 (6%), were recurrent in WTS analysis, and these fusions were detected mostly in non-high-risk patients with ganglioneuroblastoma histology. Analysis of risk-group-specic biomarkers showed that several genes and gene sets were differentially expressed between the risk groups, and some immune-related pathways tended to be activated in the high-risk group. Mutational signatures 6 and 18, which represent DNA mismatch repair (MMR) associated mutations, were commonly detected in 60% of the patients. In the tumour mutation burden (TMB) analysis, four patients showed high TMB (> 3 mutations/Mb), and had mutations in genes related to either MMR or homologous recombination. Excluding four outlier samples with TMB > 3 Mb, high-risk patients had signi�cantly higher levels of TMB compared with the non-high-risk patients. Conclusions: This study provides novel insights into the genomic background of MYCN non-amplied NB. Activation of immune-related pathways in the high-risk group and the results of TMB and mutational signature analyses collectively suggest the need for further investigation to discover potential immunotherapeutic strategies for


Background
Neuroblastoma (NB), the most common extracranial solid tumour in children, accounts for 6 to 10% of all childhood cancers.NB arises from precursor cells of the sympathetic nervous system and adrenal medulla [1].The clinical course is highly heterogeneous, ranging from spontaneous regression without therapeutic intervention to rapid progression to death, despite modern intensive multimodal treatment regimens.Thus, clinical and biological factor-based risk strati cation and tailored treatment approaches have been the mainstay of NB treatment.International Neuroblastoma Risk Group (INRG) de nes the high-risk group to include patients with MYCN ampli ed tumours and patients > 18 months old with metastatic tumours [2].Ampli cation of the MYCN oncogene is the rst genetic marker reported to indicate highly aggressive and advanced-stage NB.It is observed in approximately 20% of cases and remains a powerful prognostic factor, indicating adverse clinical outcomes [3].
The clinical features of MYCN-ampli ed NB have been attributed to the biological consequence of MYCN ampli cation.MYCNampli ed tumours make up about 40% of high-risk NBs [4], indicating that 60% of high-risk NBs are MYCN non-ampli ed tumours.Despite the extensive study of the genomic characteristics of high-risk NB including MYCN-ampli ed tumours [4][5][6], genomic pro ling of MYCN non-ampli ed NB, including low-and intermediate-risk NB, has been limited.Immunotherapy, which includes the use of immune checkpoint inhibitors, has become a potential therapeutic option, especially in adult oncology, and tumour mutational burden (TMB) is known to be a predictive marker for immunotherapy in many studies [7,8].However, except for a monoclonal antibody that acts against the tumour-associated disialoganglioside, GD2 [9], little is known about immunotherapy in NB.Here, we examined the genomic pro les of MYCN non-ampli ed NB and studied risk-groupspeci c biomarkers, TMB, and mutational signature to identify biomarkers for the molecular strati cation of NB.

Characteristics of patients
WES and WTS were conducted for 70 and 63 NB samples, respectively.QC ltering and removal of MYCN-ampli ed NB patient data yielded 58 WES and 48 WTS samples (Figure 1a).All patients in this study were East Asian.Thirty-ve patients (53.8%) were diagnosed as metastatic, and 26 patients (40%) were classi ed into the high-risk group.The median age was 3.1 years, with a range of 0 to 14.9 years.Seven patients (10.8%) experienced recurrence at median 1.7 (0.2-3.8) years (Table 1).S1A).Frequently mutated genes were summarized in Additional le 1: Figure S1B.The most frequently mutated gene, MUC4, was found to be mutated in 26% of samples, followed by RBMXL3 (19%), ALB (17%), MUC16 (14%), and SEPD8 (14%) (Figure 1b).In comparison with 9 MYCN-ampli ed tumours, there was no statistically signi cant difference between the mutation frequencies of single genes (Figure 1c).However, mutations in mucin family genes such as MUC4, MUC16, and MUC17 were more frequent in MYCN non-ampli ed subjects.
An association between risk groups and genomic variants was not observed (Additional le 2: Figure S2A).To identify the associations between altered pathways and risk groups, 17,810 annotated gene sets were analysed, and mutation status was determined in each pathway using MSigDB v6.2.Among 7,044 pathways acquired from the BIOCARTA, KEGG, REACTOME, and Gene Ontology (GO) databases, 48 pathways had a P-value < 0.05 in Fisher's exact test that became insigni cant after multiple FDR corrections (Additional le 2: Figure S2B).Alterations in metabolic pathways were enriched in the high-risk group.

Gene fusions
Gene fusions predicted by three or more algorithms were considered to be true positives.Among 48 WTS samples, 21 gene fusions were detected in 15 samples (Table 2).CCDC32-CBX3 fusion recurred in ve samples, while SAMD5-SASH1 fusion recurred in three samples (Additional le 3: Figure S3A).The existence of recurrent fusions was veri ed by Sanger sequencing.
These two recurrent fusions were detected only in ganglioneuroblastoma (GNB) histology.Most of the recurrent fusions were detected in non-high-risk patients, aside from one patient who had four other fusions, including CCDC32-CBX3 and SAMD5-SASH1 fusions.In patients with SAMD5-SASH1 fusion, SAMD5 and SASH1 were upregulated; however, this correlation was not observed in CCDC32-CBX3 fusion (Additional le 3: Figure S3B).Transcriptome analysis to identify risk-group speci c biomarkers WTS data of 48 patients showed a correlation between gene expression patterns and the risk-group identity (Additional le 4: Figure S4A).Forty-six genes were signi cantly over-expressed in the high-risk group and 40 genes were over-expressed in nonhigh-risk group (P-value < 0.05 and absolute fold change > 2) (Figure 2a and Additional le 4: Figure S4B).FAM153A (SAMD15) and FAM15B (TMED8) were the most signi cantly over-expressed genes in the high-risk group.Figure 2 also suggested that there were subgroups in the non-high-risk group, so unsupervised clustering of non-high-risk group was performed.It revealed two distinct subgroups showing clear differences in age, pathologic differentiation and stage (Additional le 4: Figure S4C).
The GSVA score was computed for 17,810 gene sets, and an association test was performed for differentially expressed gene sets between high-and non-high-risk groups (Additional le 5: Figure S5A).In total, 44 gene sets were signi cantly different between the risk groups (P-value < 0.05 and absolute mean difference > 0.3) (Figure 2b and Additional le 5: Figure S5B).The gene sets did not differ in most of the acquired canonical pathways; only 15 of these pathways showed statistically signi cant differences (Additional le 5: Figure S5C).In the high-risk group, the pathways of ketone body metabolism and mitochondrial fatty acid beta oxidation were inactivated, while the pathways of TALL-1, regulation of MHC class II biosynthesis, and regulation of interferon gamma secretion were activated (Additional le 5: Figure S5D).The ganglioside biosynthesis pathway showed correlations with risk group identity and GNB histology (P-value= 0.0002) (Additional le 5: Figure S5E).
Analyses of mutational signature and tumour mutation burden (TMB) A mutational signature analysis was performed with WES data (Additional le 6: Figure S6A).Each sample was assigned to the most predominant signature among the 30 signatures (MS-1 to MS-30).Sixty percent of samples were assigned to either MS-15 or MS-6, which were denoted as MMR signatures (Figure 3a).Association was not found between MMR signatures and mutation incidence in MMR-related genes such as MLH1, MSH2/6, or PMS2.Only ve samples (5.58%) were assigned to MS-18, a known NB signature.The association between MS-15 and GNB histology was identi ed (P-value = 0.0012) (Figure 3b).
The median TMB was 0.66 Mb (Figure 3c) in all patients; the speci c median values were 1.04 Mb in the high-risk group and 0.53 Mb in the non-high-risk group (P-value = 0.195).Excluding four outlier samples with TMB > 3 Mb, high-risk patients had signi cantly higher TMB compared with the non-high-risk patients (0.95 Mb vs 0.60 Mb, P-value = 0.001) (Figure 3d).Among the four patients whose tumours had TMB > 3 Mb, two patients belonged to the high-risk group and the other two belonged to the non-high-risk group.One sample from the high-risk group, having two missense mutations (Q1285K and T290K) in POLE and one missense mutation (G357W) in MLH1, showed a high TMB value of 10.66 Mb.Another high-risk patient with a high TMB value of 3.83 Mb had mutations in ATR, ATRX, POLQ, RAD54L, and SPIDR, all of which play roles in DNA repair or homologous recombination.Further, two samples from the non-high-risk group showed relatively high TMB, at 5.1 Mb and 4.33 Mb.Both patients were diagnosed at a very young age (under two months).Although the tumours of both patients showed high microsatellite instability (MSI), their underlying mutation pro les were different-POLE splicing variant was detected in one patient, while the other patient had deleterious mutations in RAD51AP1, RAD51B, and RMI2, which are involved in homologous recombination de ciency.
In terms of the association between TMB and mutational signatures, the sample with extremely high TMB (>10 Mb) showed MS-18, while the three samples with relatively high TMB (> 3 Mb) showed MS-6, one of the four MMR signatures.(Figure 3e and Additional le 6: Figure S6B).

Discussion
The genomic characteristics of MYCN non-ampli ed NB were identi ed in this study using WES and WTS.MYCN-ampli ed NB was excluded, because MYCN ampli cation is a well-known prognostic factor in NB, and gene expression in MYCN-ampli ed NB is quite different from that in MYCN-non-ampli ed NB.A total of 26 high-risk patients were without MYCN ampli cation, and these patients accounted for 40% of all MYCN-non-ampli ed patients in our dataset.
Common variants of high-risk NB, such as ALK, ATRX, PTPN11, NRAS, and MYCN mutations [5], were not noticeable.Instead, after eliminating the effect of MYCN ampli cation and considering low-VAF mutations, a number of novel recurrently mutated genes were found.Although the recurrent mutations did not show strong patterns of association with the different risk groups, their roles in MYCN-non-ampli ed NB warrant further exploration.The mutation pro les of MYCN-ampli ed and non-ampli ed patients were compared, and mutations in mucin family genes were found to be more frequent in the MYCN-non-ampli ed subjects.Although mutations in the mucin gene family have been reported in NB [10], their biological relevance to NB remains unclear.
Two recurrent fusions, CCDC32-CBX3 and SAMD5-SASH1, were newly detected in this study.No recurrent fusion has been reported in NB, with the exception of fusions including the NBAS gene in MYCN-ampli ed tumours [5].All but one of the patients in this study who presented recurrent fusions fell under the non-high-risk category; thus, it is likely that these fusions have not been reported because most of the previous studies included only high-risk patients.These fusions have been detected in other cancers [11,12] and further research is needed to investigate the potential roles of these fusions in the tumorigenesis of NB.
In the analysis of risk-speci c biomarkers, several genes and gene sets were differentially expressed between the risk groups.
Speci cally, some immune-related pathways, such as regulation of MHC class II biosynthesis and regulation of interferon gamma secretion, tended to be activated in the high-risk group.In this study, high-risk patients had higher TMB values compared with the non-high-risk patients (when four outlier samples with TMB > 3 Mb were excluded), which could be a factor that causes the activation of immune-related pathways in the high-risk group.The underlying mechanism of this nding remains to be elucidated.
Notably, this study presents several ndings in support of the possible application of immunotherapy in NB.In the TMB analysis, a subset of patients was found to have much higher TMB values than the other patients.One sample had TMB > 10 Mb, and three more had a moderate threshold of > 3 Mb.All four of these tumours had mutations in DNA mismatch repair de ciencyrelated genes or genes involved in homologous recombination de ciency.MSI was high in the two non-high-risk patients.
Additionally, in the mutation signature analysis, 65% of tumours showed MMR signatures when each sample was designated to the most predominant signature out of all 30 signatures.PD-1/L1 expression, TMB, and MSI have been considered as predictive biomarkers for immunotherapy in many studies [7,8,[13][14][15], and the ndings of this study suggest the possibility of immunotherapy introduction in a subset of patients with NB.
Despite comprehensive analysis, this study has several limitations.In mutational signature analysis, MS-18, a known NB signature, was present in only a few samples.Mutational signatures were calculated based on a pattern of 96 base substitution combinations, so an insu cient number of mutations may have affected the analysis.The median number of mutations, at 34.5, was relatively small.Therefore, the mutational signatures of patients with lower numbers of variants may fail to represent all of the characteristics.Since the number of variants and TMB in childhood cancers are smaller compared to those in adult cancers [16][17][18][19][20], the results of mutational signature analysis need to be interpreted with caution.Furthermore, MSigDB contains pathways with large numbers of genes, and the pathways investigated here had gene sets with up to hundreds of genes.
Therefore, it is necessary to verify the effects of individual mutations.

Conclusions
In conclusion, this study provides novel insights into the genomic background of the MYCN-non-ampli ed NB population.Activation of immune-related pathways in the high-risk group and the results of TMB and mutational signature analyses collectively suggested the need for further investigation to discover potential immunotherapeutic strategies for NB.

Study population and data collection
From November 2015, tissue and blood samples were collected prospectively from NB patients undergoing biopsy.Samples from patients who were diagnosed before November 2015 that had been deposited at the Samsung Medical Center Biobank were also included.Medical records regarding age, sex, stage, risk group, pathology, and outcome were collected.Tumour staging was determined by following the International Neuroblastoma Staging System standards [2].MYCN ampli cation was determined by performing interphase uorescence in situ hybridization on tumour tissues.Patients older than 18 months and in stage four malignancy and patients with MYCN-ampli ed tumours were strati ed as high-risk patients.

DNA and RNA extraction
All tumour specimens were reviewed by a pathologist to determine the percentage of viable tumours and their adequacy for sequencing.Genomic DNA from the tissue and blood was extracted using a QIAamp DNA Mini Kit (Qiagen, Valencia, CA, USA).
The total RNA from the same fresh frozen tumour tissues was extracted with an RNeasy Mini Kit (Qiagen, Valencia, CA, USA), according to the manufacturer's instructions.The quality and quantity of extracted nucleic acids were evaluated using Nanodrop 8000 UV-Vis spectrometer (NanoDrop Technologies Inc., Wilmington, DE, USA), Qubit ® 3.0 Fluorometer (Life technologies Inc., Carlsbad, CA, USA) and 4200 TapeStation (Agilent Technologies Inc., Santa Clara, CA, USA).Specimens with a yield over 100 ng were selected for whole exome sequencing (WES) and whole transcriptome sequencing (WTS).Those with a median DNA fragment size of 350 bp and an RNA integrity number (RIN) of 5 were selected.

WES and variant calling
Tumour and matched normal DNA were enriched for exon regions, using the SureSelect XT regent kit (Agilent Technologies Inc., in Illumina HiSeq 2500 (Illumina Technologies Inc., San Diego, CA, USA).Unresolved bases in FASTQ les were trimmed, reads were aligned to the human reference genome, hg19, using TopHat version 2.0.6 [29], and reference-guided assembly of transcripts was performed using Cu inks version 2.1.1 [30].Alignment quality was veri ed with SAMtools version 0.1.19[22].
Gene expression was estimated from the RNASeq data of 56 patients using a count-based method with RSEM [31].In total, 20,345 protein-coding genes were selected.Further, genes that were expressed in at least three samples were retained.A total of 16,120 genes were analysed.Gene counts were used as the input for Trimmed Mean of M value (TMM) normalization in the R package, edgeR [32], and normalized counts were transformed to log2-counts per million (logCPM) using the voom application in the R package, limma [33].

Validation of fusions by RT-PCR and Sanger sequencing
The putative gene fusions, detected by RNA-Seq, were veri ed by reverse transcription PCR (RT-PCR), followed by Sanger sequencing.cDNA was synthesized from 2 µg total RNA using a QuantiTect Reverse Transcription kit (Qiagen Inc., Hilden, Germany) with primers that ank the breakpoint of the fusion, in DNA Engine Tetrad 2 Peltier Thermal Cycler (BIO-RAD, Hercules, CA, USA) with the following cycling conditions: one cycle of 5 minutes at 95 °C, followed by three-step cycles of 30 seconds at 95 °C, 30 seconds at 62 °C, 10 minutes at 72 °C, and a nal extension for 20 minutes at 72 °C.PCR products were puri ed using a Multiscreen lter plate (Millipore Corp., Bedford, MA, USA) and sequenced in an ABI prism 3730XL Analyzer (Thermo Fisher

Supplementary Files
This is a list of supplementary les associated with this preprint.Click to download.

Table 1
Demographics of study populations for WES and WTS analysis Mutation pro les of MYCN non-ampli ed NB WES data of 58 patients were analysed.The median number of variants per sample was 34.5, with a range of 11-537 (Additional le 1: Figure