Meta-mining of Neuroblastoma and Neuroblast Gene Expression Profiles Reveals Candidate Therapeutic Compounds

Purpose: Neuroblastoma is a heterogeneous childhood tumor with poor survival outcome for the aggressive type despite intensive multimodal therapies. In this study, we aimed to identify new treatment options for neuroblastoma based on integrative genomic analysis.

Experimental Design: The Connectivity Map is a database comprising expression profiles in response to known therapeutic compounds. This renders it a useful tool in the search for potential therapeutic compounds based on comparison of gene expression profiles of diseased cells and a database of profiles in response to known therapeutic compounds. We have used this strategy in the search for new therapeutic molecules for neuroblastoma based on data of an integrative meta-analysis of gene copy number and expression profiles from 146 primary neuroblastoma tumors and normal fetal neuroblasts.

Results: In a first step, a 132-gene classifier was established that discriminates three major genomic neuroblastoma subgroups, reflecting inherent differences in gene expression between these subgroups. Subsequently, we screened the Connectivity Map database using gene lists generated by comparing expression profiles of fetal adrenal neuroblasts and the genomic subgroups of neuroblastomas. A putative therapeutic effect was predicted for several compounds of which six were empirically tested. A significant reduction in cell viability was shown for five of these molecules: 17-allylamino-geldanamycin, monorden, fluphenazine, trichostatin, and rapamycin.

Conclusions: This proof-of-principle study indicates that an integrative genomic meta-analysis approach with inclusion of neuroblast data enables the identification of promising compounds for treatment of children with neuroblastoma. Further studies are warranted to explore in detail the therapeutic potential of these compounds.

Neuroblastoma is the most common extracranial solid childhood tumor with a remarkable variation in clinical presentation ranging from favorable localized tumors that can spontaneously regress to highly metastatic disease with unfavorable outcome. Despite intensive multimodal therapies, survival rates for aggressive neuroblastomas are still disappointingly low. One possible approach toward development of more effective and less toxic therapies is to gain insight into the biochemical pathways that are deregulated in neuroblastoma and to use this information for the design of molecular therapies. However, at present, only three genes, MYCN, PHOX2B, and ALK, have been directly linked to neuroblastoma pathogenesis (i.e., MYCN amplification in a subgroup of aggressive metastasizing tumors, PHOX2B germ-line mutations in a subset of familial cases and in a small percentage of sporadic cases, and ALK mutations in ∼12% of high-stage sporadic cases; refs. 1–5). One of the challenges in unraveling the altered signaling in neuroblastoma cells is the genetic heterogeneity of the tumor. Genome-wide DNA copy number and transcriptome analysis are now offering the possibility to investigate this genetic heterogeneity in much more detail. Array comparative genomic hybridization (CGH) studies identified three major genetic subgroups within sporadic neuroblastoma cases (i.e., triploid tumors without segmental alterations, MYCN-amplified neuroblastoma with 17q gain and 1p deletion, and MYCN single-copy neuroblastoma with 11q deletion and 17q gain; refs. 6, 7). The relationship between these subgroups and gene expression profiles has thus far not been investigated thoroughly and most expression studies aimed particularly at the identification of prognostic genes (8–14).

Disturbance of normal differentiation pathways is an important event in oncogenesis of multiple tumor types, as exemplified by disruption or mutation of hematopoietic differentiation genes in leukemia and implication of the patched-hedgehog signaling pathway in brain tumors (15, 16). The proven role of PHOX2B and ALK, both involved in nervous tissue differentiation, in the oncogenesis of familial neuroblastoma and the occurrence of somatic ALK and PHOX2B mutations in sporadic neuroblastoma tumors indicates that future data-mining efforts for neuroblastoma gene expression profiles should ideally also include data from the normal fetal counterpart of neuroblastoma cells, such as the adrenal neuroblasts. The limited accessibility to these fetal human neuroblasts has, however, precluded such studies in the past. Using dedicated protocols for isolation and amplification of mRNA, we have been the first to determine the expression profiles of these normal neuroblasts, present as small cell clusters in fetal adrenal glands (17). These data strongly supported the hypothesis that the neuroblasts are the “neuroblastoma cells of origin” and now open the way for more powerful data mining of the available tumor expression profiles. In a first comparative expression analysis of the normal neuroblasts versus the malignant neuroblastomas, several neuroblastoma candidate genes, including ALK, were identified.

In this study, we investigated whether the genetic heterogeneity, reflected in the three major genetic subgroups, of neuroblastoma tumors is also reflected in their expression profile. For each subgroup, we could identify a signature with genes differentially expressed between normal progenitor neuroblasts and malignant neuroblastomas. Instead of looking for signaling pathways with therapeutic potential that are enriched in these gene lists, we directly submitted the lists to the Connectivity Map database (18). This public resource contains gene signatures obtained after treatment of cells with several chemical compounds and has pattern-matching tools to detect similarities between these signatures and gene lists uploaded by the user. Application of this tool on the subgroup-specific signatures directly allowed the identification of new candidate therapeutic small molecules for neuroblastoma.

Data preprocessing. Raw data (.Cel files) of four published neuroblastoma expression microarray studies obtained from two different platforms were collected and reanalyzed [i.e., hgu95av2 Affymetrix (11, 19) and hgu133a with (17) and without (20) preamplification]. The Bioconductor repository of the statistical language R was used to do gcRMA normalization (21). In the next step, probe IDs of both platforms were matched using the array comparison spreadsheet.³

Significance analysis of microarrays (22) was used to identify probe IDs that are differentially expressed between hgu133a with or without preamplification as well as probe IDs that are differentially expressed between hgu133a and hgu95av2 arrays. These probe IDs were omitted for further meta-analysis (Supplementary Data 2).

Statistical analysis. For the construction of a prediction analysis of microarrays (PAM) classifier, the MCRestimate Bioconductor package (9) was used and rank products (RP) analysis was done with the RankProd Bioconductor package (23).

Neuroblastoma gene server analysis and positional gene enrichment analysis were done using in-house developed tools (24).⁴

^,5

For the Connectivity Map analysis, RP lists were loaded in the Connectivity database (build 1). In the field for up-regulated genes, the lists that are higher expressed in the neuroblasts compared with the neuroblastomas were uploaded and vice versa for the down-regulated gene field. Using the gene set enrichment analysis algorithm (25), cmap calculates enrichment scores for the different compounds. Compounds with significant P values (<0.1) after permutation analysis were listed in Table 3 (in this way, only compounds from the database that were tested more than twice were included).

Cell viability assay. Cells were seeded in duplicate in 96-well plates (10,000 per well), incubated for 6 h to allow adherence to the surface, and then treated in duplicate for 24 and 48 h with six-point dilution series of the compounds dissolved in either ethanol [trichostatin (Sigma-Aldrich) and rapamycin (LC Laboratories)] or DMSO [17-allylamino-geldanamycin (17-AAG; Sigma-Aldrich), monorden (Calbiochem), fluphenazine (Sigma-Aldrich), and colforsin (Sigma-Aldrich)]. Concentration ranges were according to the literature for trichostatin A (0-8 μmol/L; ref. 26), rapamycin (0-25 μmol/L; ref. 27), 17-AAG (0-8.1 μmol/L; ref. 28), and monorden (0-5 μmol/L; ref. 29). For fluphenazine, we used a concentration range of 0 to 20 μmol/L, as this compound was not yet used for treatment of cancer cells. Cell viability was determined using the CellTiter-Glo Luminescent Cell Viability Assay (Promega) with untreated cells as a reference (average signals were plotted).

Quantitative PCR validation. Cell lines IMR-32, NGP, CLB-GA, and SK-N-SH were treated with 2.5 μmol/L 17-AAG and harvested at 24, 48, and 72 h after treatment. In parallel, untreated cells were harvested at each time point.

RNA was isolated from the cells using the miRNeasy kit (Qiagen). Subsequently, 2 μg RNA from each sample was treated with RQ1 DNase I (Promega) and desalted using a Microcon-100 spin column (Millipore). cDNA synthesis was done on the eluate with the iScript cDNA synthesis kit (Bio-Rad). All manipulations were conducted according to the manufacturer's instructions. Reverse transcription-PCR amplification reactions were carried out in a 384-well plate with a total reaction volume of 8 μL, containing 10 ng of template cDNA, 4 μL of 2× SYBR Green I reaction mix (Eurogentec), 4 μL nuclease-free water (Sigma-Aldrich), and 250 nmol/L of each primer. Cycling conditions were as follows: 10 min at 95°C followed by 45 cycles of denaturation (10 s at 95°C) and annealing (45 s at 60°C) and elongation (1 s at 72°C). All reactions were done on LC480 (Roche). Primers for CLCN6 (RTPrimerDB-ID7804), CYP1B1 (ID7805), DDIT4 (ID7806), EGR3 (ID7807), and HSPA6 (ID7808) were designed and validated using RTPimerDB's batch primer design tool (30–32). LC480 Cq export files were imported into qBasePlus (Biogazelle; ref. 33) and normalized using five stably expressed reference genes: AluSq, HMBS (ID4), HPRT1 (ID5), SDHA (ID7), and UBC (ID8; ref. 34).

Genomic subgroup classification of the neuroblastoma tumors. For the integrative meta-analysis, mRNA expression microarray data of four published studies were merged (11, 17, 19, 20). Integration of genomic and transcriptomic data of the tumors was conducted through differential expression profiling of the neuroblastoma tumors that belonged to the different genomic subclasses. Based on (array)CGH and fluorescence in situ hybridization data, neuroblastoma samples of these studies were classified in subtype 1 (no MYCN amplification, whole 17 gain, no structural aberrations), subtype 2A (11q deletion, no MYCN amplification, 17q gain or normal 17), and subtype 2B (MYCN amplification, no 11q deletion, no 3p deletion, 17q gain or normal 17) tumors. One hundred forty-six of 207 tumors (70%) could be unequivocally classified in one of the three prototype subgroups (Table 1). The remainder could not be classified as either the genomic information was incomplete or missing (6 from study 2 and 20 from study 3) or a genomic aberration pattern was present characteristic for more than one subgroup, confirming previous studies (6, 35). For this analysis, we included only tumors unequivocally assigned to the subgroups defined above.

Gene expression classification of genomic subgroups. We first examined whether the genomic aberration pattern of neuroblastoma tumors is reflected in their mRNA expression profile. Therefore, we established a classifier aiming to distinguish neuroblastoma tumors of the three different genomic subgroups. A complete 10-times repeated 10-fold cross-validation (9, 36) was done on the merged data sets to estimate the PAM classification accuracy of the patients in one of the three subgroups. This estimation reflects the validity of the expression data of the different studies for prediction of the genomic subgroup using a PAM classifier. Three patients of each subgroup from each study (total = 36) were used in this training step and showed that the estimated performance of classification in one of the three subgroups using PAM is very good (accuracy of 92% = 33 of 36; Table 2A).

We selected 136 probes (132 unique genes) that were included in at least 95 of the 100 predictive gene lists generated in the cross-validation process, as these genes are likely to be the ones with the highest predictive value (9). We used an in-house developed neuroblastoma gene server⁴ to compare this gene list with genes that have been reported as differentially expressed in 29 published gene expression profiling studies on neuroblastoma (total of 2,864 genes after exclusion of the published gene lists from the four data sets used in this study). We found that as much as 67 of the 132 PAM genes (50%) are included in neuroblastoma gene server. The finding of these genes in other neuroblastoma expression profiling studies strongly suggests that these genes are genuinely deregulated in neuroblastoma. Nine were included in the Cancer Gene Census list (ALK, ATIC, CBFB, ETV1, GMPS, IGL@, MYCN, NTRK1, and RET), 15 are nervous system development genes (Gene Ontology: 0007399; ALCAM, ALK, ASCL1, CDH4, DRD2, GAL, GAP43, NRCAM, NRG1, NTRK1, PMP22, PTN, SH3GL3, SLIT3, and SPOCK2), 11 are cell cycle genes (Gene Ontology: 0007049; BCAT1, BUB1, CADM1, CCNB1, CHEK1, CKS2, DLG7, DUSP4, E2F3, KIAA0367, and PTN), and 5 are apoptosis genes (Gene Ontology: 0006915; CADM1, KIAA0367, IL7, PMAIP1, and SCG2).

Using positional gene enrichment analysis (24), we searched for chromosomal loci that are overrepresented in the classifier list. Significant enrichment was found for a region on 1p36 that is frequently deleted in high-stage neuroblastoma tumors (ATP6V0B, CDC42, CLCN6, CLSTN1, FUCA1, GNB1, PINK1, PRDM2, PTPRF, STX12, and VPS13D; adjusted P < 0.002).

In the next step, the performance of the established PAM classifier was tested on the remaining 110 samples. This resulted in an accuracy of 91% (Table 2B), which is comparable with the estimated accuracy of 92%. The results showed that the gene classifier most easily identified subtype 2B tumors, whereas discrimination of subtype 1 and 2A tumors was more difficult. These data clearly show that the genomic differences in neuroblastoma tumors are reflected in their expression profiles and that the classification of neuroblastoma tumors in these three different genomic subgroups might be relevant toward the biology of neuroblastoma and must be taken into account in downstream gene expression analysis.

Genes differentially expressed between neuroblasts and neuroblastomas from different genomic subgroups. To gain further insights in the genes that play a role in the genetic subgroups of neuroblastoma, genes were identified that are differentially expressed in the normal neuroblasts (17) versus the neuroblastoma samples of three genomic subgroups using the RP algorithm (P < 0.05; ref. 23). More detailed information on the gene lists can be found in Supplementary Data 1. Comparison of the different gene lists showed a high overlap between the lists but also indicated that there are many subgroup-specific genes (Fig. 1). This analysis further confirms that tumors of subgroup 1 and 2A share more genes (86 + 11) differentially expressed compared with normal neuroblasts than subgroup 1 and subgroup 2B tumors (26 + 9) or subgroup 2A and subgroup 2B tumors (27 + 14).

Connectivity Map analysis of differentially expressed genes and cell viability screening of candidate therapeutics. The RP gene lists were analyzed using the Connectivity Map (cmap) database (build 1; ref. 18). This database contains genome-wide transcriptional data from cultured human cells treated with bioactive small molecules. Submission of the gene lists to cmap allows to find similarities between the transition of the malignant neuroblastoma phenotype to the normal neuroblast phenotype and the change of cancer cells after treatment with a certain compound. It is hypothesized that a similar transition profile indicates a potential therapeutic effect of the compound on neuroblastoma growth. Table 3 summarizes the significant results (permuted results with P < 0.1) for the different gene lists (neuroblasts versus neuroblastoma type 1, neuroblasts versus neuroblastoma type 2A, and neuroblasts versus neuroblastoma type 2B) and ranks the compounds according to the mean enrichment score across the different signatures. In the top list, there was a positive enrichment for four HSP90 inhibitors, three histone deacetylase inhibitor compounds, and three dopamine blocking agents. For functional screening, we selected representative compounds of each of these classes (with the most significant P values; monorden and 17-AAG, trichostatin A, and fluphenazine) as well as the compound with the highest overall enrichment score (forskolin) and rapamycin that showed a significant enrichment for all three signatures. The added value of our strategy in which we contrasted tumoral expression profiles characteristic of the different genomic subgroups with the normal neuroblast profiles is clear when doing cmap analysis on the list of genes that are differentially expressed between the neuroblasts and all neuroblastoma samples, irrespective of their genomic subgroup. This analysis provided us with only six potential therapeutic compounds (data not shown), of which four were in common with the compounds listed in Table 3. This is in sharp contrast to the 9, 14, and 13 compounds that are significantly enriched in the 1, 2A, and 2B signatures, respectively.

Viability assays on five neuroblastoma cells [CLB-GA, IMR-32, SKNBE (2c), NGP, and SK-N-SH] showed that five of six selected compounds have indeed an effect on neuroblastoma cell viability (Fig. 2). As a negative control, we observed that two compounds (iloprost and estradiol) that were not significant according to the cmap analysis had no effect on neuroblastoma cell growth (Supplementary Data 3).

To check whether the compounds are indeed acting via activation or repression of genes from the gene lists (used for the cmap analysis), we selected the top ranking genes from the cmap results, that is, the highest median ranking score for the different cmap treatments for one compound (i.e., 17-AAG). Expression of five presumed differentially regulated genes (three up-regulated and two down-regulated genes) was tested using real-time quantitative reverse transcription-PCR on four cell lines and three time points. It was observed that four of the five genes (CLCN6, DDIT4, HSPA6, and CYP1B1) show the expected effect (i.e., up-regulated or down-regulated for all three time points) in at least one of the four cell lines (Supplementary Data 4).

Despite recent advances in chemotherapeutic and transplantation options, the overall prognosis for advanced neuroblastomas remains very poor. For those surviving, the intensive treatment schemes may lead to both short-term and long-term side effects, including risk for secondary tumors. Therefore, more effective and less toxic drugs are urgently needed. In this study, we aimed to identify new therapeutic compounds for treatment of neuroblastoma patients through an integrative genomic meta-analysis and comparative expression analysis with normal neuroblast cells. The success of this study is based on three important and innovative data-mining steps. First, we used publicly available expression data to allow a more robust statistical analysis. Four published expression microarray data sets generated on two different Affymetrix platforms were merged and reanalyzed (11, 17, 19, 20). In this way, more than 200 neuroblastoma tumors were analyzed. Second, in one of the studies, the transcriptome of fetal adrenal neuroblasts, the normal counterpart cells of the malignant neuroblastomas, was profiled (17). Inclusion of these data is of particular interest as comparison with the expression profiles from the malignant tumors may lead to the identification of developmental pathways that are disturbed and can be therapeutically targeted. Third, we did a unique integrative genomics approach. In previous integrative genomics efforts, tumors with and without a given genomic aberrations were compared (19, 20, 37, 38). Here, we focused on the global genomic landscape of the tumors. As previously described (6, 7), three major genomic subgroups can be recognized (i.e., a group of tumors with only numerical aberrations; a group of tumors with MYCN amplification, 1p deletion, and 17q gain; and a group characterized by 11q and 3p deletion and 17q gain). In fact, also a fourth subgroup of patients with few or no genomic aberrations can be recognized. This last subgroup might represent tumors with few or no genomic alterations, such as the acute myeloid leukemias with normal karyotypes that carry nucleophosmin mutations (39). These neuroblastomas might also be characterized by specific mutation patterns. The possibility of normal cell contamination should also be considered but is less likely, as in most studies samples with >60% of tumor cell content were used. Based on (array)CGH and fluorescence in situ hybridization results, the tumors were classified in one of the three major genomic subgroups. In a first step, we were able to establish a 132-gene classifier that allows to distinguish tumors that belong to the different genomic subgroups based on mRNA expression data, in keeping with an assumed association of specific differential gene expression profiles for each of the three major genomic neuroblastoma entities. Not surprisingly, this classifier contains numerous genes that have been reported in the context of neuroblastoma (e.g., MYCN, NTRK1, ALK, CADM1, and ASCL1). In the next step, the neuroblast expression profile was compared with the expression profiles of neuroblastomas from the different subgroups. Instead of looking at genes and pathways that are disrupted and might be therapeutically targeted, we directly submitted the differentially expressed gene lists to the Connectivity Map database (cmap; build 1). In this way, we were able to identify agents that may have therapeutic potential in neuroblastoma. We showed that only by considering the genomic subgroup in which the tumors belong, we could identify many more potential therapeutic compounds. The finding of different compounds with the same mode of action further supported the validity of our approach. Moreover, for several of these compounds, a therapeutic activity has been reported previously. Histone deacetylase inhibitors (e.g., trichostatin A) were shown to have promising therapeutic activity for different cancer entities (40). Several studies have also shown the growth-suppressive activity of these inhibitors on neuroblastoma (41, 42). In addition, two HSP90 inhibitors came up from the data-mining screen, of which the 17-AAG molecule was already shown to have therapeutic potential in neuroblastoma (43, 44). Moreover, this compound was evaluated in a phase I study on pediatric tumors (45). Another HSP90 inhibitor and macrocyclic antifungal antibiotic, monorden (alias radicicol), with the ability to suppress transformation by diverse oncogenes such as Src, Ras, and Mos (46), was not yet tested in neuroblastoma cells. In addition, the antiproliferative or antiapoptotic activity of rapamycin, an inhibitor of mammalian target of rapamycin, and phenotiazines, such as fluphenazine, has been reported for neuroblastoma (27, 47, 48). To evaluate in vitro therapeutic potential of the selected compounds, six of the strongest candidate compounds were tested by assessing the viability of neuroblastoma cells on treatment. This analysis revealed that five of these compounds effectively reduced neuroblastoma cell viability using concentrations in the low micromolar range, suggesting that our approach provides a rational means for prioritization of small-molecule compounds for preclinical evaluation.

Recently, cmap build 2 was released containing 6,100 expression profiles from 1,309 bioactive small molecules (only 164 compounds in the first build of cmap). Reanalysis of our gene lists in the new build leads to the identification of other promising compounds with therapeutic potential in neuroblastoma (data not shown). Using the same criteria, 160 potentially interesting compounds were identified, including 17 of 20 from the list in Table 3. This small discrepancy might be due to the different statistical algorithm that is used in the new version to identify significantly enriched compounds. One of the top five compounds from this analysis is quinostatin, which is, like rapamycin, an inhibitor of the mammalian target of rapamycin pathway and might be an interesting therapeutic (49, 50). Further analysis is warranted to test this compound in neuroblastoma.

The presented compound screening approach is of particular interest for pediatric tumors such as neuroblastoma because numerous obstacles exist in the development of therapies for these entities not in the least the rare occurrence of this disease. Typically little or no efforts are undertaken for drug development for such rare disorders due to the limited expected return. With so few people affected by neuroblastoma, there is a reduced market incentive for industry-based drug development. Therefore, it is of interest to screen a list of small-molecule compounds of which the majority is already approved by the Food and Drugs Administration and which do not need further preclinical trials to evaluate toxicity.

No potential conflicts of interest were disclosed.

We thank Nurten Yigit for technical assistance with cell culturing and reverse transcription-PCR analysis.

European Community under the FP6 (project: STREP: EET-pipeline, number: 037260), “Kinderkankerfonds,” “Stichting tegen Kanker,” Fund for Scientific Research, and BOF. K. De Preter and T. Van Maerken are supported by the Fund for Scientific Research.