Integrative Genomics Identifies Distinct Molecular Classes of Neuroblastoma and Shows That Multiple Genes Are Targeted by Regional Alterations in DNA Copy Number

Neuroblastoma is an important childhood cancer, as it accounts for ∼15% of all pediatric oncology deaths (1). The clinical hallmark of neuroblastoma is heterogeneity, with the likelihood of tumor progression varying widely according to age and disease burden at diagnosis. A large body of data support the hypothesis that the clinical behavior of human neuroblastoma may be reliably predicted based on the analysis of a panel of prognostic variables (reviewed in ref. 2). The Children's Oncology Group (COG) currently stratifies patients into low-, intermediate-, or high-risk categories based on analysis of well-defined prognostic factors, including patient age at diagnosis (3), International Neuroblastoma Staging System (INSS) stage (4), tumor histopathology (5), DNA index (6), and MYCN amplification status (7, 8). This stratification system is currently used to assign therapeutic intensity, mandating that the algorithm be as precise as possible. Clinical experience with this system suggests that the algorithm is useful, but misclassifications almost certainly occur resulting in overtreatment or undertreatment (2). Additional tumor-specific prognostic markers may be required to achieve maximal predictive power.

A large number of genomic aberrations have been defined in neuroblastoma, and the pattern of these somatically acquired changes correlate with tumor behavior (9). Amplification of the MYCN oncogene at 2p24 provides a paradigm for the clinical utility of a tumor-specific DNA rearrangement. MYCN amplification is present in ∼20% of newly diagnosed tumors and remains a powerful biomarker for an aggressive phenotype and poor survival probability even in light of modern dose-intensive therapies (10, 11). Hemizygous deletions and/or loss of heterozygosity (LOH) have been defined at multiple genomic regions, including chromosomal arms 1p, 3p, 4p, 11q, 14q, 16p, and 19q at frequencies between 20% and 45% of cases using a variety of techniques (2, 12). Likewise, unbalanced gain of 17q material may be the most common genomic aberration in neuroblastoma primary tumors, as it is present in 50% to 75% of cases. Perhaps most importantly, these events are nonrandom and are clinically relevant, with LOH at 1p36 and 11q being independently prognostic for disease outcome (13). Despite the clear clinical relevance of these genomic aberrations, no definitive neuroblastoma suppressor gene has been identified.

Genome-scale profiling technology allows for detailed reassessment of the relationship between DNA aberrations with regional gene expression (14, 15). Microarray studies in neuroblastoma to date have mainly focused on defining prognostic signatures or identifying genes associated with oncogenic MYCN (16–18). In addition, McArdle et al. correlated global gene expression profiles with 11q deletions in a limited number of samples but showed the power of this combined approach for gene discovery (19). We therefore sought to determine regional and global gene expression differences between a large number of samples with and without clinically relevant chromosomal aberrations in primary neuroblastomas. These data were used to further our understanding of the genetic basis for neuroblastoma heterogeneity and to prioritize regional candidate genes for further analyses.

Patient samples and molecular characterization. A total of 101 diagnostic neuroblastoma primary tumor samples were selected from the COG (n = 91 prospectively collected) or Children's Hospital of Philadelphia (n = 10) neuroblastoma tumor banks. All specimens were snap frozen in the institutional pathology departments immediately on arrival according to standard COG neuroblastoma biology protocols. Specimens were examined for percent tumor content before RNA extraction and excluded if <70% neuroblastic content were estimated. Each tumor specimen was also assayed for MYCN copy number by fluorescence in situ hybridization, DNA index by flow cytometry, histopathology according to the system of Shimada, and allelic status at chromosome bands 1p36, 11q23, and 17q23 using standard techniques as described previously (refs. 5, 6, 20, 21). We differentiate 11q LOH into two categories based on our prior data (13): balanced 11q LOH = LOH at both 11q and 11p polymorphic loci and unbalanced 11q LOH = loss of 11q loci with no loss at 11p loci. RNA was isolated from ∼50 mg snap-frozen tumor tissue using TRIzol reagent (Invitrogen, Carlsbad, CA) and purified with the RNeasy System (Qiagen, Valencia, CA) according to the manufacturer's recommendations and digested with DNase I (Ambion, Houston, TX). RNA quantity and quality were determined by spectrophotometry, 1% agarose gel, and microfluidics-based electrophoresis (Agilent 2100 Bioanalyzer, Agilent, Palo Alto, CA). Fetal brain RNA (Clontech, Palo Alto, CA) was used as a control sample. The Children's Hospital of Philadelphia Institutional Review Board approved this study and informed consent was obtained before collection of all samples.

An independent set of 72 primary neuroblastomas obtained at diagnosis for patients treated at Memorial Sloan-Kettering Cancer Center (MSKCC) was also made available for this study (16, 22). Specimens were procured and processed in a manner identical to that described above. Samples were analyzed for LOH at 1p36 and amplification of MYCN similar to as described above. Samples showing discordance of MYCN copy number and expression values were excluded from this sample set (22).

Expression profiling. Total RNA (5-15 μg) was used to generate cRNA targets that were hybridized to Affymetrix U95Av2 oligonucleotide arrays containing 12,625 probe sets according to the manufacturer's instructions (Affymetrix, Santa Clara, CA). The quality of labeled cRNA target was assessed using Affymetrix GeneChip Test 3 arrays. Statistical modeling of probe set behavior was conducted using Probe Profiler (Corimbia, Berkeley, CA) as described by the manufacturer. The software implements a model-based approach for data extraction from probe sets and incorporates quality-control features to detect and correct for the contribution of nonspecific cross-hybridization effects. The model weights probe pairs based on consistency of performance. Normalization to remove chip-to-chip variation is accomplished by subtracting the second percentile intensity value and dividing by the interquartile range. A quantitative expression score (e-score) is then calculated for each probe set.

Data from this experiment, including sample annotations, have been deposited in Gene Expression Omnibus database and are available for free at http://www.ncbi.nlm.nih.gov/geo/.

Unsupervised clustering. Unsupervised two-way agglomerative hierarchical clustering with Pearson's correlation coefficient as a distance metric was done using GeneSpring 7.0 (Silicon Genetics, Redwood City, CA). Probe sets with low average expression or low variance were excluded in an effort to minimize noise. Specifically, probe sets exhibiting an average e-score <50 or an estimated SD <75 were excluded, leaving 2,177 probe sets remaining for cluster analysis.

Identification of enriched gene sets. Gene Ontology (GO) annotation was used to characterize the functional relationship of differentially expressed gene classes between the varying sample classes (23). EASE version 2.0 was selected to examine the GO categories of molecular function, biological process, and cellular component (24). EASE functions to determine the overrepresentation of GO terms within a given gene list. For each of the three GO categories, EASE establishes the total number of annotated terms within the input list and population as well as the abundance of a specific term within the input list and population. Each gene class was processed against the population list of 2,177 probe sets and overrepresented GO terms were statistically determined by the one-sided Fisher's exact probability. GO terms with P ≤ 0.05 were selected as significant. In addition, a bootstrap false discovery rate (FDR) was calculated to ensure that the P of the GO terms remained significant over multiple tests. This multiple testing correction creates random gene lists from the population whose size is equal to that of the input list. These random lists allow for the determination of the frequency of false-positive terms in the input list with Ps greater than or equal to the term of interest. In addition to having a P ≤ 0.05, GO terms were required to have a bootstrap FDR ≤0.05.

Differential gene expression. Differential gene expression was measured with both Patterns from Gene Expression (PaGE; ref. 25) and Significance Analysis of Microarray (SAM; ref. 26) algorithms. Samples were first segregated into those with and without a detectable genomic alteration at the region being studied. Binary comparisons were then done for each region and results were visualized with heat maps. For visualization purposes, e-scores were standardized by subtracting the average and dividing by the SD estimated from those samples without a detectable genomic alteration. Standardized e-scores were then binned into five levels based on the number of SD from the mean (>1.5 SD below the mean, 1-1.5 SD below the mean, within 1 SD of the mean, 1-1.5 above the mean, and >1.5 SD above the mean). Heat maps of the standardized and binned e-scores were then generated using Spotfire (Spotfire DecisionSite 8.0, Somerville, MA). Transcriptome maps were generated by placing the Affymetrix U95Av2 GeneChip probe sets located on 1p, 11q, 2p24-25, and 17q in linear order according to chromosomal map position using the May 2004 assembly of the human genome according to the University of California at Santa Cruz Genome Browser (hg 17; http://genome.ucsc.edu; ref. 27).

Reverse transcription-PCR. Real-time quantitative reverse transcription-PCR (RT-PCR) was done as described previously (28). All primer/probe sets spanned exon boundaries to assure specificity for cDNA and primer sequence data are available upon request. Relative expression of target gene was determined by normalization to glyceraldehyde-3-phosphate dehydrogenase (GAPDH) using a standard curve method with 10 serial dilutions according to the manufacturer's instruction. All RT-PCR experiments included a no template control and were done in triplicate. The Wilcoxon's rank-sum test was used to test for statistically significant differences in expression based on class label. Correlations between RT-PCR and microarray data were tested with the Spearman's rank correlation.

Sample selection strategy and quality control. The patient and tumor characteristics for the 101 samples analyzed on the 12,625 oligonucleotide probe set Affymetrix U95Av2 array are shown in Table 1. Samples were selected so that a minimum of 20 cases would be available in four clinically and biologically distinct subsets of neuroblastoma: low-risk, intermediate-risk, high-risk, and high-risk with MYCN amplification. In addition, we attempted to limit biological heterogeneity within each subset by imposing further restrictions on sample selection. For example, although low-risk tumors can be INSS stages I, IIA, IIB, or IVS from patients of a variety of ages, we only selected INSS stage I cases for patients ages <18 months for this cohort because this is the most common and unambiguous representation of low-risk disease. Region-specific allelotyping data are also shown in Table 1 and are consistent with past reports in terms of correlations with high-risk disease features (2, 13, 29). Relationships between expression data and patient survival were not a focus of this report because the median duration of follow-up for the 80 patients who had not died was only 2.3 years.

Quality control was assured by using the standard Affymetrix algorithms as well as metrics available in Probe Profiler. Unscaled brightness was within 2.5 SD of the mean across all experiments. Four primary neuroblastomas (one from each class) were randomly selected from the set of 101 neuroblastomas and assayed twice on the Affymetrix array to assure consistency of results. Pairwise Spearman's rank correlation coefficients for each replicate comparison were 0.97 to 0.99 (data not shown). The standard intensity values from the two housekeeping genes GAPDH and β-actin (ACTB) showed very tight interexperimental consistency and the vast majority of signals were within 2.5 SD of the mean across all experiments.

Class discovery. To characterize the molecular phenotype of human neuroblastoma, we first identified 2,177 probe sets that showed the highest degree of variability across the entire sample set. Unsupervised two-way agglomerative hierarchical clustering results are shown in Fig. 1. There were five main sample classes discovered in this analysis and these were clearly distinguished from the control RNA sample derived from fetal brain. Sample class 1 (S1) was composed almost exclusively of cases with genomic amplification of MYCN (17 of 18) and 1p36 LOH (15 of 18). The only sample in this class that did not have amplification of MYCN showed LOH at 1p36. All of these samples would be classified as high-risk in the current COG classification schema due to MYCN amplification status or, in the one discordant sample, metastatic disease in a patient older than 1 year at diagnosis. Sample class S4, on the other hand, was composed of 28 cases that were almost exclusively without metastatic disease and few or no genomic alterations at the loci assayed. All but one of these samples would be classified as low or intermediate-risk by the current COG system. Nine of the 12 cases in the S2 cluster (n = 12) were older than 1 year with metastatic disease at diagnosis, but all showed normal MYCN and 1p36 status. Sample classes S3 (n = 28) and S5 (n = 15) were the most heterogeneous in terms of known clinical and genomic covariates, as 17 of 28 and 7 of 15 would be classified as high-risk, respectively. Whereas there was only one case with MYCN amplification in each of these sample clusters, the majority showed genomic alterations at 1p36, 11q23, and/or 17q23.

Examination of the gene dendrogram showed three distinct classes, and these correlated strongly with genomic aberrations (Fig. 1). Gene cluster 1 (G1) contained 907 probe sets that show overexpression in the S1 sample group with MYCN amplification. Genes located within or nearby the core region of the MYCN amplicon all were located in G1 (MYCN, NCYM, DDX1, and NAG), and this cluster included multiple bona fide MYCN transcriptional target genes, such as ODC1 and MCM7 (30). This cluster was enriched for genes involved in cellular metabolism, nucleic acid, protein and macromolecule biosynthesis, RNA processing, ribosomal turnover, and RNA processing (Supplementary Table S1).

A second cluster (G2; Supplementary Table S2) contained 437 genes and was enriched for those involved in host immune response and antigen processing. High expression of these genes was associated with sample cluster 2, composed mainly of high-risk cases, but without MYCN amplification. A tight subcluster containing several MHC class II genes showed nearly identical expression patterns across the sample set (Fig. 1).

The final cluster (G3; Supplementary Table S3) of 741 genes was enriched for those highly expressed in the more benign subsets of neuroblastoma. GO terms overrepresented included cellular communication and processing, vesicle-mediated transport, signal transduction, neurogenesis, and neurophysiologic processes. Genes known to be involved in normal sympathetic nervous system development, such as TH, DBH, PHOX2B, GATA3, and NTRK1, all showed higher expression in the nonmetastatic cases within this cluster. This cluster was also enriched for genes mapping to chromosome arm 1p, such as STX12, PTP4A2, GNB1, CLSTN1, PUM1, RERE, CHD5, and BAI2, and these showed lower expression in the metastatic cases with 1p36 LOH. The nonreceptor tyrosine kinase FYN also showed differential expression in a similar pattern.

Region-specific alterations in gene expression. We created transcriptome maps to visualize gene expression pattern based on chromosomal map position (Fig. 2). All genes present on the Affymetrix U95Av2 GeneChip located on 1p, 2p24-25, 11q, and 17q were displayed in linear order according to genome map position. Samples were then segregated into those with and without a detectable genomic aberration at the region being studied. There were a large number of genes showing differential expression at each of the loci under study. For examples, the presence of 1p36 LOH was associated with low or absent expression of hundreds of genes mapping telomeric to cytoband 1p21 (Fig. 2A). Samples with MYCN amplification showed a very tight region of high overexpression noted at 2p24, with MYCN, NCYM, DDX1, and NAG showing striking differential expression based on the presence or absence of the amplification event (Fig. 2B). Samples with unbalanced 11q LOH showed the most marked differential expression of 11q14-q25 mapped transcripts, but many of the samples with balanced 11q LOH also showed a similar pattern (Fig. 2C). A similar pattern in the opposite direction is noted at 17q, but this is less robust (Fig. 2D).

Identification of genes differentially expressed based on genomic alterations. We next determined which genes were differentially expressed throughout the genome in each binary comparison of genomic status at 1p36, 2p24, 11q23, and 17q23 using both PaGE and SAM algorithms. We found no difference between these two techniques (data not shown) and only consider the PaGE data here (full output from this analysis can be visualized at http://stokes.chop.edu/programs/marislab/data). Summary PaGE statistics are shown in Table 2. On a genome scale, the binary comparisons of normal versus aberrant at 1p36, MYCN, 11q, and 17q showed differential expression of 4,386, 6,269, 293, and 3,360 probe sets, respectively. For each analysis, the majority of the most statistically significantly differentially expressed genes in the expected direction (lower in LOH cases for 1p and 11q and higher in gain cases for 2p24/MYCN region and 17q) mapped to the region being studied. Many of the genes showing decreased expression in the 11q analyses also mapped to 3p, a region known to be coordinately deleted in many cases with 11q deletions (31). The much lower number of differentially expressed genes in the chromosome 11 analysis did not change substantively if we considered unbalanced 11q LOH versus all others (305 differentially expressed probe sets) or the comparison of unbalanced 11q LOH versus balanced 11q LOH versus normal 11q status (452 differentially expressed probe sets).

We next focused only on the genes mapping within the genomic region that defined the class labels of aberrant versus normal in the differential expression analyses. Of the 201 genes mapping to 1p35-36, 122 (61%) showed significantly lower expression in the 26 cases with 1p36 LOH compared with only 5 1p35-36 genes that showed significantly higher expression. Similar results were seen at each of the other loci: 27% of the 11q mapped genes showed significantly decreased expression in the 40 cases with 11q LOH, 79% of the 2p24-25 mapped genes showed significantly increased expression in the 20 cases with MYCN amplification, and 40% of the 17q mapped genes showed significantly increased expression in the 38 cases with unbalanced 17q gain and very few genes showed significant differential expression in the opposite direction (Table 2).

To prove that the mRNA copy number estimates from the microarray accurately reflected true mRNA copy number in the tumor cells, we did quantitative RT-PCR for 18 genes shown to be differentially expressed in these analyses. All Spearman's rank correlation coefficient Ps comparing microarray and RT-PCR data were <0.0001 (data not shown). A representative sampling of these data segregated by genomic alterations is shown in Fig. 3.

To prioritize regional candidate neuroblastoma suppressor genes at 1p and 11q, we arranged the probe sets in rank order by their t statistic for differential expression based on allelic status. Figure 4A and B shows the top 50 differentially expressed genes mapping to 1p35-36 and 11q, respectively. This visualization tool allows one to identify samples that show mRNA expression patterns that are discordant with their class assignment. For example, samples 1,009, 1,043, and 1,068 showed no evidence for 1p36 LOH yet showed expression patterns more similar to the samples with 1p36 LOH. Each of these samples also showed MYCN amplification. In contrast, several of the samples coded as having 1p36 LOH but not MYCN amplification showed expression patterns more similar to the cases without 1p36 LOH (e.g., 58, 205, 1,118, and 1,204).

To determine if our methods reproducibly identify transcripts that can be considered as regional candidate tumor suppressor genes or oncogenes, we used an independent data set annotated with 1p36 and MYCN genomic data (Table 1). Of the 122 genes showing significant differential expression in the expected direction at 1p36 in our data set, 71 (58%) also showed significant differential expression in the same direction in the independent validation data set (Table 2). Not surprisingly, concordance was highest for the most significantly differentially expressed genes (Table 3). Table 3 shows the 1p35-36 transcripts concordant in both data sets as ranked by t statistic. Figure 5 compares data for the most differentially expressed genes at the MYCN locus at 2p24-25 in both data sets, again showing a high degree of concordance.

Neuroblastoma has served as a paradigm for the clinical importance of tumor genomic data (2, 13, 29). MYCN oncogene amplification is reliably associated with an aggressive clinical phenotype and remains a powerful prognostic maker even in the era of modern dose-intensive therapy. However, despite a wealth of data on the clinical relevance of other chromosomal aberrations in this disease, no bona fide neuroblastoma suppressor gene has been unequivocally identified. The fact that the majority of hemizygous deletions mapped in neuroblastoma cells are typically dozens of megabases in size led us to speculate that multiple genes are likely targeted by these acquired genomic alterations. This has important implications for gene discovery efforts and for ongoing studies directed toward identifying genetic biomarkers of disease outcome and/or new therapeutic targets. Thus, the primary aim of this study was to combine DNA copy number and mRNA expression data to understand the magnitude of transcriptional alterations at candidate tumor suppressor gene and oncogene loci.

Class discovery methods have been used to differentiate distinct human cancers (32, 33) and more recently have been applied to subclassify within a given histopathology (34–36). It is well known that human neuroblastomas show dramatic clinical and biological heterogeneity. Using unsupervised class discovery methods, we now show that there are distinct expression profiles that correlate with clinically and biologically relevant subsets of neuroblastoma. The most striking mRNA signature is embedded in the cases with MYCN amplification (Fig. 1). These cases almost all cluster together. One case that clustered elsewhere in the sample dendrogram showed significantly lower MYCN mRNA copy number, whereas the other two were the only two infants with MYCN amplified tumors in this study. In general, MYCN-amplified samples showed overexpression not only of genes mapping within or near the typical core MYCN amplicon (37), such as MYCN, NCYM, DDX1, and NAG, but also of known transcriptional targets, such as ODC1 and MCM7 (30, 38). Genes mapping to chromosome 17q were also enriched, and this is likely due to the strong correlation of 17q gain with MYCN amplification (39) and the fact that the NME1 and PPM1D genes likely have an oncogenic role during neuroblastoma progression (40–42). In contrast, these samples showed low or absent expression of a large number of genes mapping to chromosome 1p, a genomic region that is deleted in the vast majority of samples also with MYCN amplification (13, 43).

A novel finding was the highly enriched immunity signature, including much higher expression of HLA molecules in a distinct subset of high-risk neuroblastomas without MYCN amplification (Fig. 1, sample class S2). These were all locally invasive (INSS stage III) or metastatic tumors, and most had evidence for copy number alterations at 11q and/or 17q but not 1p. Although the duration of follow-up is short, it is intriguing that only 2 of 12 patients in this subset have died from tumor progression, although 9 of these cases would be classified as high-risk using the current COG system with expected survival probabilities of <40% (44). Although it is certainly possible that this signature is due to tumor-infiltrating lymphocytes rather than the tumor tissue itself, significant lymphocytic infiltration of neuroblastomas have generally been observed in low-stage cases with the associated opsoclonus-myoclonus syndrome (45–47). We speculate that evidence for immune activation may have an influence on tumor biology, but this very preliminary observation requires careful follow-up. In contrast, low-risk tumors had an expression signature enriched in genes involved in normal sympatheticoadrenal developmental pathways, such as PHOX2B, GATA3, NTRK1, TH, DBH, and FYN, suggesting that loss of a neuronal differentiation/maturation pathway may be critical for developing a more malignant tumor phenotype.

Our gene expression profiling data suggest three distinct classes based on mRNA signatures (S2, S3, and S5) with no clear distinguishing features in terms of known clinical and biological covariates. It will be important to determine, as the data set matures, whether there is a correlation with tumor clinical behavior and patient outcome. It is tantalizing to speculate that the S2 class samples with the strong immune signature might be correlated with improved antitumor control and high survival probability, but the relatively small subset sample sizes and short follow-up duration preclude these analyses at the present time.

Deletions of chromosome arms 1p and 11q are observed in neuroblastomas with a more malignant phenotype and detection of these genomic alterations is being applied in the clinic as prognostic biomarkers (13). Despite enormous effort to define the genes targeted by these deletions, no definitive 1p or 11q neuroblastoma suppressor gene has been identified. Our data are consistent with the hypothesis that perhaps dozens of genes mapping to 1p35-36 and 11q14-25 contribute to the malignant phenotype of high-risk neuroblastomas. By evaluating a completely independent group of samples for which 1p allelic status was known, we could show tight concordance in the 1p genes differentially expressed in both data sets (Table 3). It is also reassuring that the vast majority of 1p genes identified here were also seen by Janoueix-Lerosey et al. as candidate neuroblastoma suppressors using a similar study design in a smaller sample set (48). Specifically, 10 of the 11 1p genes in their study (GNB1, CLSTN1, CDC42, MFN2, SLC35E2, CAPZB, TCEB3, VAMP3, RBAF600, and STMN1) that were also represented on our microarray platform were significantly differentially expressed in both our primary and validation data sets. Caren et al. recently identified several genes at 1p36.22 that show lower expression in neuroblastomas with a poor prognosis, and two of these genes (UBE4B and CORT) were also identified as differentially expressed based on 1p36 allelic status in this study (49). CASP9 was differentially expressed in both data sets, but its role in neuroblastoma biology has been controversial (50, 51). Lastly, CAMTA1 has recently been postulated to function as an oligodendroglioma suppressor gene through haploinsufficiency (52), and its differential expression discovered in this study implies a role in neuroblastoma suppression as well.

Classic positional cloning efforts may be streamlined by considering genome-scale mRNA copy number data. We recently mapped a 729-kb consensus region of deletion at 1p36.3 in a large series of primary neuroblastomas (53). This region overlapped and was partially defined by constitutional deletions in two children with developmental abnormalities who subsequently developed neuroblastoma, suggesting that a gene within this region is involved in neuroblastoma initiation. Of the 24 known transcripts within this interval, CHD5, ICMT, TNFRSF25, KIAA0469, MGC33488, and CAMTA1 each showed highly significant decreased expression in both primary and validation data sets. Likewise, several other positional candidate 1p neuroblastoma suppressor genes have been postulated, and the CDC2L2 (PITSLRE), PRDM2 (RIZ1), CDC42, and CHD5 genes were each highly differentially expressed in both data sets (41, 54–56), supporting a role for each gene in tumor evolution. Taken together, these data support the hypothesis that transcriptional alterations in multiple genes located at the distal short arm of chromosome 1 cooperate during neuroblastoma initiation and/or progression.

There was remarkable concordance between the two data sets regarding transcripts mapping to 2p24-25 that are overexpressed in samples with MYCN amplification (Fig. 5). These data are reassuring for ongoing efforts to map DNA copy number aberrations using mRNA expression profiles (15). Whereas the alterations in mRNA transcription secondary to high-level genomic amplification at a particular locus seem very reproducible despite variations in amplicon size, the downstream consequences of such an event are myriad. Over half of the transcriptomes under study (6,269 probe sets; Table 2) showed significant differential expression in the binary comparison of samples dichotomized on MYCN amplification status. The distribution was approximately equal between genes showing higher or lower expression in the MYCN-amplified cases (3,033 versus 3,236 probe sets, respectively). These data further support the notion that the Myc proteins are promiscuous transcriptional activators and repressors with pleiotropic and diverse effects on cell function.

Only 40% of high-risk neuroblastomas show MYCN amplification, suggesting that alternative pathways to a malignant phenotype exist. Deletion of chromosome arm 11q has recently emerged as a critical genomic event in the evolution of high-risk neuroblastomas that do not have MYCN amplification (13, 57). However, no 11q neuroblastoma suppressor gene has been identified. Like 1p deletions, 11q deletions are often large spanning over 60 Mb from 11q13 to 11qter (58), and the only known constitutional deletions at 11q are likewise very large (59, 60). Whereas a critical region has been mapped to 11q23 (61), several tumors show deletions that map outside of this region (13), and functional studies have suggested other chromosome 11 loci critical for suppression of differentiation (62). Whereas our independent data set was not annotated for 11q allelic status, our primary data allow for prioritization of potential 11q neuroblastoma suppressor genes. Known tumor suppressor genes, such as IGSF4 (63), SDHD (64), MLL (65), and PPP2R1B (66), showed low or absent expression in the majority of cases with 11q LOH. A similar pattern was seen for genes involved in response to DNA damage (ATM and H2AFX; ref. 67), cell adhesion (NCAM1; ref. 68) and nonreplicative senescence (LOH11CR2A; ref. 69). Each could be considered candidate neuroblastoma suppressor genes based on our expression data, but only the SDHD gene has been studied in detail in this disease (70).

Epigenetic modification of DNA is an important somatically acquired mechanism to silence gene expression in human cancer. Several examples exist in neuroblastoma, including dense hypermethylation of promoter sequences for several genes, including CASP8 (71), RASSF1A (72), TSP1 (73), NR1I2 (74), etc. It has also been shown that a CpG island hypermethylation phenotype is associated with more aggressive clinical behavior and poorer survival probability (75, 76). It is certainly possible that many of the genes showing low expression in regions of hemizygous deletion in this study are subjected to epigenetic silencing of the retained allele. Ongoing studies will clarify the frequency of epigenetic versus genetic (mutation) events to result in biallelic inactivation as well as to define if the second event is not required for some genes and that a phenotypic effect can be manifested simply from haploinsufficiency of the deleted allele.

Oncogenic activation of key signaling pathways may provide the most realistic target for pharmacologic manipulation. Although the MYCN oncogene is an obvious drug target, the profound degree of transcriptional deregulation makes complete inhibition challenging, and this abnormality is found in only 20% of all neuroblastoma cases. Unbalanced gain of 17q material is nearly ubiquitous in high-risk neuroblastomas (39, 57, 77), and we postulate that combining genomic and transcriptomic data may help prioritize potential molecular targets for rationale drug design. Our data are consistent with previous reports that NME1 and PPRM1D are overexpressed in neuroblastomas with 17q gain (40–42) but also identify multiple other potential oncogenes, including some that are potentially drugable tyrosine kinases (78), such as MAP3K3 and TLK2. Combining array comparative genomic hybridization and expression profiling data sets at higher resolution genome-wide will provide a much more refined list of molecular targets that may be exploitable.

In conclusion, this study shows that the genomic data embedded in neuroblastoma diagnostic biopsies can be used to subcategorize the disease into molecular subsets. Our data also clearly show that the regional copy number alterations are correlated with a broad array of transcriptional alterations genome-wide. For LOH and low copy number gain events that are generally large enough to be visible at the karyotype level, the number of genes residing within these regions that show significant mRNA copy number alterations is daunting. Our data suggest that multiple genes from several discrete regions of the human genome cooperate to suppress neuroblastoma tumorigenesis and progression. DNA alterations likely target multiple genes; therefore, detection of copy number alterations may provide more prognostically relevant information than mRNA expression signatures. Ongoing studies will test whether DNA or RNA alterations are more useful in a clinical setting or whether the two will provide complementary information.

Grant support: NIH grants R01-CA78545 and R01-CA87847 (J.M. Maris), P01-CA97323 (G.M. Brodeur, S. Shusterman, and J.M. Maris), P01-CA106450 (N-K. Cheung and W. Gerald), R01-CA39771 (G.M. Brodeur), and U10-CA98543 (COG), Alex's Lemonade Stand Foundation (J.M. Maris), Hope Street Kids Foundation (Q. Wang and J.M. Maris), and Abramson Family Cancer Research Institute (J.M. Maris).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.