Association of high-risk neuroblastoma classification based on expression profiles with differentiation and metabolism

Neuroblastoma, the most common extracranial solid malignancy among children, originates from undifferentiated neural crest cells (NCC). Despite recent intensified treatment, high-risk patients still have a high mortality rate. To explore a new therapeutic strategy, we performed an integrated genomic and transcriptomic analysis of 30 high-risk neuroblastoma cases. Based on the expression profiling of RNA sequencing, neuroblastoma was classified into Mesenchymal (MES; n = 5) and Noradrenergic (ADRN; n = 25) clusters, as previously reported in the super-enhancer landscape. The expression patterns in MES-cluster cases were similar to normal adrenal glands, with enrichment in secretion-related pathways, suggesting chromaffin cell-like features built from NCC-derived Schwann cell precursors (SCPs). In contrast, neuron-related pathways were enriched in the ADRN-cluster, indicating sympathoblast features reported to originate from NCC but not via SCPs. Thus, MES- and ADRN-clusters were assumed to be corresponding to differentiation pathways through SCP and sympathoblast, respectively. ADRN-cluster cases were further classified into MYCN- and ATRX-clusters, characterized by genetic alterations, MYCN amplifications and ATRX alterations, respectively. MYCN-cluster cases showed high expression of ALDH18A1, encoding P5CS related to proline production. As reported in other cancers, this might cause reprogramming of proline metabolism leading to tumor specific proline vulnerability candidate for a target therapy of metabolic pathway. In ATRX-cluster, SLC18A2 (VMAT2), an enzyme known to prevent cell toxicity due to the oxidation of dopamine, was highly expressed and VMAT2 inhibitor (GZ-793A) represented significant attenuation of cell growth in NB-69 cell line (high SLC18A2 expression, no MYCN amplification) but not in IMR-32 cell line (MYCN amplification). In addition, the correlation of VMAT2 expression with metaiodobenzylguanidine (MIBG) avidity suggested a combination of VMAT2 inhibitor and MIBG radiation for a novel potential therapeutic strategy in ATRX-cluster cases. Thus, targeting the characteristics of unique neuroblastomas may prospectively improve prognosis.

Introduction Neuroblastoma, the most common extracranial solid malignancy among children under 15 years of age, accounts for 8%-10% of pediatric tumors [1,2]. In general, neuroblastoma is thought to originate from undifferentiated neural crest cells (NCC), which can become any of several different cell types, depending on the location within the embryo [3]. NCCs may partially differentiate into neuroblastoma or ganglioneuroblastoma, which are malignant, or may differentiate into benign ganglioneuroma [4]. The primary site of neuroblastoma is typically the adrenal medulla or tissues that originate from the sympathetic nervous system [2]. Although low-and intermediate-risk patients generally have a favorable outcome, high-risk patients show a high mortality rate and fewer than 50% of patients have long-term survival despite recent intensified treatment [5,6].
Several genome-wide analysis studies have been reported which aimed to improve prognosis and develop a new therapeutic strategy for high-risk neuroblastoma [2,7,8]. Segmental chromosomal aberration is common in high-risk neuroblastoma, such as amplification of the MYCN oncogene [9] deletions of chromosomes 1p, 3p, 4p, and 11q, and gains of chromosomes 1q, 2p, and 17q [10][11][12]. These chromosomal aberrations greatly influence the clinical course [13]. In contrast, recent high-throughput genome-wide studies have revealed that there were few recurrent somatic alterations in high-risk neuroblastoma, except MYCN [9] ALK [14][15][16][17] ATRX [18,19] and TERT [20]. These genetic alterations do not account for the entire genetic mechanism leading to high-risk neuroblastoma. Thus, effective targeted therapies for intractable neuroblastoma remain limited.
Based on the super-enhancer landscape of neuroblastoma cell lines, there are two neuroblastoma subtypes: Noradrenergic (ADRN)-type and Mesenchymal (MES)-type, showing distinct expression patterns in core regulatory circuitry (CRC)-related genes [21,22]. However, expression profiling of high-risk neuroblastoma is not yet fully understood. In the present study, we classified 30 cases of International Neuroblastoma Staging System (INSS) [23] Stage 4 neuroblastomas, based on expression profiles of whole transcriptome sequencing (WTS). Then we conducted a combined analysis with data on mutations and copy number alteration (CNA) to explore a new therapeutic strategy for high-risk neuroblastoma.

Patients and materials
This study enrolled 30 specimens collected from pediatric patients (2-140 months) who had been diagnosed with neuroblastoma with INSS Stage 4 neuroblastoma and admitted to Tokyo University Hospital and various hospitals in Japan between January 2003 and December 2015 (S1 Table). Enrolled patients were treated under the various protocols. All patients and/or their parents provided written informed consent, and the Human Genome, Gene Analysis Research Ethics Committee of the University of Tokyo and other participating institutes approved the study protocols. Studies were conducted at Department of Pediatrics, The University of Tokyo, in accordance with the principles of the Declaration of Helsinki. Biopsy samples at the primary tumor site were collected from neuroblastoma patients. No matched normal samples were available. Genomic DNA and RNA of each sample was isolated from frozen biopsy samples by NucleoSpin DNA RapidLyse kit and NucleoSpin RNA kit (Macherey-Nagel Gmbh & Co., Düren, Germany) according to the manufacturer's protocol.
Captured targets were subjected to sequencing using a HiSeq 2000 or 2500 platform (Illumina, San Diego, CA, USA) with a standard 125-bp paired-end read protocol according to the manufacturer's instructions. With mean depths of 411 reads (range = 267-575, S3 Table), sequence alignment and mutation calling were performed using our in-house pipeline "Genomon v.2.5.0," as previously described [25] Reads with a mapping quality score of <25, a base quality score of <30, or five or more mismatched bases were excluded from the analysis. Candidate mutations with a variant allele frequency in tumor samples �0.1 and an EBcall [27] (Empirical Bayesian mutation calling) P � 1 ×10 −3 in coding regions were adopted, and filtered by excluding the following: (i) synonymous mutations and variants without complete ORF information; (ii) known variants listed in the 1000 Genomes Project (Oct 2014 release); NCBI SNP database (dbSNP) build 131, National Heart, Lung, and Blood Institute (NHLBI) Exome Sequencing Project (ESP) 6500, the Human Genome Variation Database (HGVD; October 2016 release), or our in-house SNP database; (iii) variants present only in unidirectional reads; (iv) variants occurring in repetitive genomic regions; (v) variants with <5 supporting reads in tumor samples; and (vi) all variants found in non-paired normal samples (n = 30) showing an allele frequency of >0.0025. Finally, mapping errors were removed by visual inspection on the integrative genomics viewer browser [28] Copy numbers were calculated by allele frequencies and sequence depths of SNPs using our in-house pipeline "CNACS" [29,30] We defined amplification, gain, or loss of genes and segments when calculated signal values were >10, >2.75 or <1.25 (<0.75 on the X chromosome of male subjects), respectively. Significant CNAs were identified using GISTIC 2.0 (q < 0.1).

Whole transcriptome sequencing (WTS)
Total RNA was assessed for integrity and concentration using an Agilent 4200 TapeStation system. All samples had an RNA integrity number higher than 6.5. Libraries for RNA-seq were prepared using the NEBNext Ultra RNA Library Prep kit for Illumina (New England BioLabs, Beverly, MA, USA). Normalized count data, obtained by Genomon v.2.5.0 and the R software package DESeq2, were subjected to cluster analysis. We ascertained cluster stability via consensus clustering by using the top 1,000 most variable genes on the basis of median absolute deviation, with 1,000 iterations, using the R package ConsensusClusterPlus. The R software package Rtsne was used to map the samples to tSNE plot. Heatmaps were generated by the R package pheatmap using count data. Differentially expressed genes were extracted using the R package DESeq2, and pathway analysis was performed using Metascape (http://metascape.org).

Data
The raw data of sequencing are available from the DNA Data Bank of Japan (DDBJ) (accession number hum0035, JGAS00000000246). The results published here are in part based upon data generated by the Therapeutically Applicable Research to Generate Effective Treatments (https://ocg.cancer.gov/programs/target) initiative, phs000467. The data used for this analysis are available at https://portal.gdc.cancer.gov/projects (TARGET-NBL).

Cell proliferation assay
NB-69 and IMR-32 cell lines were obtained from Dr. Yasuhide Hayashi (Jobu University). Mycoplasma infection and short tandem repeat were not tested. To select representative cell lines for each cluster (MYCN-and ATRX-cluster), we used public data of transcriptomic profiling in 39 neuroblastoma cell lines (GSE89413) [31]. Cells were seeded in triplicate at the concentrations of 3,000 cells/well in 96-wells plate. After 24 hours, cells were treated with DMSO (control), Tolcapone (COMT inhibitor, Sigma #SML0150) or GZ-793A (VMAT2 inhibitor, Sigma #SML0851). Cell proliferation was measured by using Cell Counting Kit-8 (Sigma, #96992) at 24, 48, 72, 96 hours after treatment. Cell proliferation rate was calculated by comparing the absorbance before and after the treatments.

Statistical analysis
Statistical analyses were performed using R v3.4.0 software. Adjusted P values of differentially expressed genes were calculated by DESeq2 and were considered statistically significant at adjP < 0.01. P values and q values of pathway analysis were analyzed by using Metascape and were considered statistically significant at q < 0.01. Statistical significance of gene expression level and variant allele frequency were assessed using the Wilcoxon rank-sum test. The resulting values were considered statistically significant at P < 0.05.

Unsupervised consensus clustering based on expression profiles
We performed gene expression profiling of 30 individuals with INSS Stage 4 neuroblastoma, in addition to publicly-available expression data of normal adrenal gland tissues (GSE93499, GSE88682, GSE88668), to characterize neuroblastoma in genomic terms [34] We obtained two stable clusters using unsupervised consensus clustering which validated in tSNE plot (S2 Fig). Importantly, five (17%) cases of neuroblastoma had similar expression profiles to the normal adrenal gland (Fig 2A). Among these five cases, information about the primary tumor site was available in four cases; three were in the adrenal gland, and one in the retroperitoneum (S3 Fig and S1 Table). Although samples of these five cases were obtained from a biopsy of the primary site, results of CNA and variant allele frequency detected by TCS were compatible with the rest of 25 neuroblastoma cases in the other cluster (S4 Fig and S4 Table). Therefore, these five samples seemed to have sufficient numbers of tumor cells. We concluded that the effects of contaminated normal adrenal gland cells were not strong compared to the rest of the samples.
Comparison of expression profiles of these two clusters revealed that several transcription factors involved in CRC, including ETV6 and FOSL1, were highly expressed in the cluster as previously described [22], showing similar expression patterns to normal adrenal glands ( Fig  2B and S5A Fig). Cases in the other cluster also showed a high expression of CRC-related genes in neuroblastoma, such as GATA3 and HAND2 (Fig 2B and S5B Fig). We validated the expression profiles of these two clusters using differentially expressed genes between ADRNand MES-types reported in the classification of super-enhancer-associated transcription factor networks (S5 Table) [21,22] As a result, expression patterns of the cluster showing similar expression profiles to normal adrenal glands corresponded with the MES type, whereas the other cluster corresponded with the ADRN type ( Fig 2C). These two subtypes, MES-and ADRN-clusters were validated by using expression profiles of 94 cases with high-risk neuroblastoma in the Therapeutically Applicable Research to Generate Effective Treatment (TAR-GET) database (S6 Fig and S6 Table). Thus, based on expression profiling, INSS Stage 4 neuroblastoma in the present study was classified into two subtypes, MES-and ADRN-clusters, and cases in the MES-cluster presented similar expression patterns to normal adrenal glands.

Pathway analysis in MES-and ADRN-clusters
For further characterization of MES-and ADRN-clusters, we performed pathway analysis using the 500 top-ranked differentially expressed genes between these clusters (S7 Table). In the MES-cluster, pathways related to secretion and vesicles were enriched significantly, whereas the ADRN-clusters had enriched neuron-and synapse-related pathways (Fig 2D and  2E). These results indicated that the individual cases in MES-and ADRN-clusters harbored features of chromaffin cells and the sympathetic nervous system, respectively. About 80% of chromaffin cells, which are the main components of the adrenal medulla, were reported to be generated from NCC-derived Schwann cell precursors (SCPs) later in embryonic development [35] Furthermore, SCPs and sympathoblasts originate from unique single NCCs traced from E8.5 [35] Thus, the MES-and ADRN-clusters appeared to correspond to differentiation pathways through SCP and sympathoblast cells, respectively, in the development model. These differences might be involved in expression patterns and the development of neuroblastoma of  (Table 1). In fact, all cases that originated from the neck and mediastinum, rather than the adrenal gland or the retroperitoneum, were classified into the ADRN-cluster, while all MES-cluster cases originated from the adrenal gland or the retroperitoneum (S1 Table).

Further classification of ADRN cluster based on expression profiles
We performed unsupervised consensus clustering again only for the ADRN-cluster neuroblastoma cases (n = 25) to further characterize this cluster that accounted for 80% of INSS Stage 4 neuroblastoma cases in the present study. With the top 1000 most variable genes, the ADRNcluster was divided into two groups by 2 nd -unsupervised consensus clustering and tSNE plot (S7 Fig). Based on genetic alterations detected by TCS, these two subtypes were characterized into MYCN-and ATRX-clusters (Fig 3A). The MYCN-cluster contained all the cases with MYCN amplification, as well as frequent alterations of ALK and ARID1A including focal and broad aberrations in the coding region. The ATRX-cluster contained all of the cases with ATRX abnormalities and frequent ATM alterations involving coding region, although no cases with MYCN focal amplification were classified into this cluster. CNA revealed frequent segmental chromosomal aberration in the ATRX-cluster, especially deletions of chromosome 11q, including the ATM locus, and gains of chromosomes 2p and 17q (Fig 3B and S7 Table). In contrast, MYCN-cluster featured deletions of chromosome 1p with a gain of chromosome 2p including the MYCN and ALK loci (Fig 3C and S7 Table). Thus, neuroblastoma cases in the ADRN-cluster were classified into two genetically distinct clusters, MYCN-and ATRX-clusters.

Differentially expressed genes between MYCN-and ATRX-clusters
Pathway analysis using the 500 top-ranked differentially expressed genes between the MYCNand ATRX-clusters (S8 Table) revealed significant enrichment of ribosome-related pathways in the MYCN-cluster (Fig 4). This result was consistent with the fact that MYC family proteins were regulators of ribosome biogenesis and translation [36] In addition to MYCN and ribosome-related pathway genes, ALDH18A1 was expressed significantly in the MYCN-cluster compared to the ATRX-cluster (Fig 4A). ALDH18A1 encodes P5CS, an enzyme involved in the conversion of glutamate to pyrroline-5-carboxylate (P5C). This enzyme is important for producing proline from glutamate, as well as PYCR1, which also highly expressed in the MYCN-cluster (Fig 4B). These results might indicate reprogramming of proline and glutamine metabolism, increasing proline biosynthesis in the MYCN-cluster.
In the ATRX-cluster, SLC18A2 encoding the protein VMAT2 was significantly expressed (Fig 4) compared to the MYCN-cluster. VMAT2 is a vesicular monoamine transporter that accumulates cytosolic monoamines, such as dopamine, into synaptic vesicles. This action was consistent with the pathway analysis results, which enriched vacuole-and vesicle-related pathways in the ATRX-cluster (Fig 4C and 4D). Therefore, metabolic pathways related to dopamine in the ATRX-cluster might be distinct from other neuroblastoma cases.

Targeting the dopamine related metabolic pathways
To explore the vulnerability in the metabolic pathways related to dopamine in the ATRX-cluster, we examined the efficacy of GZ-793A (VMAT2 inhibitor) and/or Tolcapone  degradation of monoamines [37,38]. Consistent with the result from expression profiles, GZ-793A treatment significantly suppressed cell proliferations in NB-69 but not in IMR-32 ( Fig  5). In contrast, Tolcapone treatment attenuated cell proliferation only in higher concentration (20μM) but not 10μM for both NB-69 and IMR-32 (Fig 5). Combination of GZ-793A and Tolcapone represented efficacy even in lower concentration of Tolcapone (10μM) in IMR-32, however, this synergistic effect was not observed in higher concentration (20μM) (Fig 5).

Discussion
In the present study, we classified INSS Stage 4 neuroblastoma into two subtypes, MES-and ADRN-clusters, based on WTS expression profiles. These two clusters were validated with

PLOS ONE
Expression profiles of HR neuroblastoma independent large cohort and consistent with the previous classification of the super-enhancer landscape of neuroblastoma cell lines [21,22]. A small set of core transcription factors forms an interconnected auto-regulatory loop, CRC, which is often driven by super-enhancers [39]. ADRN-type neuroblastoma, showing sympathetic noradrenergic identity, was characterized by CRC modules formed by several transcription factors, such as HAND2, PHOX2B, and GATA3, leading to high expression of these genes. CRC modules, including AP-1 transcription factors, defined the MES-subtype [21,22]. Because the formation of these CRCs is associated with differentiation, the pathway analysis results between MES-and ADRN-clusters might provide information about the origin of neuroblastomas. In the present study, enrichment of secretion-and vesicle-related pathways in the MES-cluster represented the features of chromaffin cells, composing about 80% of the adrenal medulla. In contrast, neuron-and axonrelated pathways were enriched in the ADRN-cluster, which indicated sympathetic nervous system features. Furthermore, primary lesions in MES-cluster cases were limited to the adrenal gland or the retroperitoneum, whereas ADRN-cluster cases occurred in tissues originating from the sympathetic nervous system. These results were consistent with previous reports showing that most chromaffin cells originated from SCPs rather than from sympathoblasts, although both SCPs and sympathoblasts are NCC-derived [35]. Thus, there was a strong correlation between differentiation and classification of super-enhancer or expression profiles in neuroblastoma. Primary neuroblastoma is a mixture of both MES-and ADRN-type cells, with a balance toward the ADRN type, in most samples [21,22]. Because of intratumor heterogeneity and the usage of bulk biopsy samples, the classification of expression profiles in the present study might indicate the main cell type component in each neuroblastoma sample. Therefore, we need to be careful to say that our classification, MES-and ADRN-clusters, is directly associated with the origin of neuroblastoma. There are two possible reasons for the intratumor heterogeneity of neuroblastoma (containing both MES-and ADRN-type cells). First, critical alterations leading to neuroblastoma might occur in NCCs before differentiating into SCPs. In this case, mutated NCCs might differentiate into both MES-and ADRN-types. Second, critical alterations may occur in a later development stage, when NCC derivatives have been already determined to become either MES-or ADRN-type cells. In this case, either of the committed MESor ADRN-type cells might interconvert to the other type. In fact, ADRN cell lines transitioned toward MES-type profiles upon chemotherapy [21,22]. To elucidate the relevance of intratumor heterogeneity and differentiation, further analyses including single-cell analysis might be required.
Cases in the ADRN-cluster were classified into MYCN-and ATRX-clusters in the present study. MYCN-cluster cases showed high expression in ALDH18A1 (encoding P5CS) and PYCR1, important enzymes in converting glutamate to proline. In Burkitt lymphoma, MYC suppressed POX/PRODH expression and increased P5CS and PYCR1, leading to reprogramming of proline and glutamine metabolism [40]. This tumor metabolic reprogramming contributes to tumor cell proliferation despite tumor-specific proline vulnerability as a compensatory mechanism [41]. In invasive breast carcinoma and kidney cancers, extensive proline production is necessary to maintain tumorigenic growth because tumor cell proliferation depends on proline [41]. As a result, high expression of PYCR1 was induced in these cancer cells. Similar to these cancers, neuroblastoma especially in MYCN-cluster might possess a tumor-specific proline vulnerability, which could be used as a target metabolic pathway for treatment.
On the other hand, SLC18A2 encoding VMAT2 was expressed significantly in ATRX-cluster cases and its inhibition with GZ-793A (VMAT2 inhibitor) represented significant attenuation of cell proliferation. VMAT2 is a membrane protein that transports monoamines from the cytosol into synaptic vesicles. Production of monoamines such as dopamine and noradrenaline is a hallmark of neuroblastoma. These monoamines are degraded by enzymes, COMT and monoamine oxidase (MAO), and are finally excreted as monoamine metabolites homovanillic acid and vanillylmandelic acid [42,43] Excess free dopamine in the cytosol undergoes oxidation [38] producing ROS [37]. Since this generation of ROS induces cytotoxicity and neurodegeneration, COMT, MAO, and VMAT2 play an important role in preventing dopamine oxidation in dopaminergic neurons [38]. Therefore, Tolcapone, a potent COMT inhibitor to treat Parkinson's disease, was reported in neuroblastoma cell lines to induce oxidative stress leading to caspase-3-mediated apoptosis and to inhibit tumor proliferation [44]. The efficacy of Tolcapone in neuroblastoma cell lines was also demonstrated in the present study and showed synergy effect with VMAT2 inhibitor in MYCN-cluster representative IMR-32 cell line. However, this synergistic effect was not observed when treated with higher Tolcapone concentration (20μM), suggesting that MYCN-cluster phenotype might use additional pathways for their proliferation and survival because Tolcapone and VMAT2 inhibitor target the same monoamine metabolism pathway. Furthermore, a report from the Children's Oncology Group showed correlation of metaiodobenzylguanidine (MIBG) avidity with high VMAT2 expression in neuroblastoma cases without MYCN amplification [45]. Thus, the combination of VMAT inhibitor with COMT/MAO inhibitor (MYCN-cluster) or MIBG radiation therapy (ATRX-cluster) might be a potential therapeutic strategy for treating neuroblastoma cases.
In the present study, 30 cases with INSS Stage 4 neuroblastoma was classified into MESand ADRN-clusters based on expression profiles. These two clusters showed association with the differentiation process and the origin of neuroblastoma (Fig 6). Further classification of the ADRN-cluster identified MYCN-and ATRX-clusters, characterized by genetic alterations and metabolism with potential therapeutic strategy. Treatment of high-risk neuroblastoma has already been intensified to the maximum limit without exceeding toxicity levels harmful to the patients. Thus, targeting metabolic reprogramming distinct in each cluster might be helpful for the development of a new therapeutic strategy for high-risk neuroblastoma.