Neuroblastoma is remarkable for its clinical heterogeneity and is characterized by genomic alterations that are strongly correlated with tumor behavior. The specific genes that influence neuroblastoma biology and are targeted by genomic alterations remain largely unknown. We quantified mRNA expression in a highly annotated series of 101 prospectively collected diagnostic neuroblastoma primary tumors using an oligonucleotide-based microarray. Genomic copy number status at the prognostically relevant loci 1p36, 2p24 (MYCN), 11q23, and 17q23 was determined by PCR and was aberrant in 26, 20, 40, and 38 cases, respectively. In addition, 72 diagnostic neuroblastoma primary tumors assayed in a different laboratory were used as an independent validation set. Unsupervised hierarchical clustering showed that gene expression was highly correlated with genomic alterations and clinical markers of tumor behavior. The vast majority of samples with MYCN amplification and 1p36 loss of heterozygosity (LOH) clustered together on a terminal node of the sample dendrogram, whereas the majority of samples with 11q deletion clustered separately and both of these were largely distinct from the copy number neutral group of tumors. Genes involved in neurodevelopment were broadly overrepresented in the more benign tumors, whereas genes involved in RNA processing and cellular proliferation were highly represented in the most malignant cases. By combining transcriptomic and genomic data, we showed that LOH at 1p and 11q was associated with significantly decreased expression of 122 (61%) and 88 (27%) of the genes mapping to 1p35-36 and all of 11q, respectively, suggesting that multiple genes may be targeted by LOH events. A total of 71 of the 1p35-36 genes were also differentially expressed in the independent validation data set, providing a prioritized list of candidate neuroblastoma suppressor genes. Taken together, these data are consistent with the hypotheses that the neuroblastoma transcriptome is a sensitive marker of underlying tumor biology and that chromosomal deletion events in this cancer likely target multiple genes through alteration in mRNA dosage. Lead positional candidates for neuroblastoma suppressor genes can be inferred from these data, but the potential multiplicity of transcripts involved has significant implications for ongoing gene discovery strategies. (Cancer Res 2006; 66(12): 6050-62)

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 (1618). 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.

Table 1.

Patient and tumor characteristics for 101 test and 72 validation cases

CharacteristicsTest casesValidation cases
n 101 72 
Age (y)   
    <1 41 21 
    ≥1 60 51 
INSS stage   
    I 27 
    II 13 
    III 23 
    IV 50 41 
    IVS 
MYCN   
    Not amplified 81 60 
    Amplified 20 12 
Histology   
    Favorable 48 32 
    Unfavorable 50 33 
    Not available 
DNA index   
    Hyperdiploid 66 26 
    Diploid 26 12 
    Not available 34 
1p36 status   
    No LOH 74 39 
    LOH 26 18 
    Not available 15 
11q23 status   
    No loss 59 ND 
    LOH 40 ND 
    Unbalanced LOH* 17 ND 
    Not available 72 
17q23 status   
    No gain 51 ND 
    Gain 38 ND 
    Not available 12 72 
COG risk group   
    Low 28 22 
    Intermediate 21 
    High 32 30 
    High + MYCN 20 12 
CharacteristicsTest casesValidation cases
n 101 72 
Age (y)   
    <1 41 21 
    ≥1 60 51 
INSS stage   
    I 27 
    II 13 
    III 23 
    IV 50 41 
    IVS 
MYCN   
    Not amplified 81 60 
    Amplified 20 12 
Histology   
    Favorable 48 32 
    Unfavorable 50 33 
    Not available 
DNA index   
    Hyperdiploid 66 26 
    Diploid 26 12 
    Not available 34 
1p36 status   
    No LOH 74 39 
    LOH 26 18 
    Not available 15 
11q23 status   
    No loss 59 ND 
    LOH 40 ND 
    Unbalanced LOH* 17 ND 
    Not available 72 
17q23 status   
    No gain 51 ND 
    Gain 38 ND 
    Not available 12 72 
COG risk group   
    Low 28 22 
    Intermediate 21 
    High 32 30 
    High + MYCN 20 12 
*

See text for definition (unbalanced 11q LOH is a subset of cases with 11q LOH).

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.

Figure 1.

Unsupervised hierarchical clustering of 101 primary neuroblastomas and human fetal brain. Unsupervised two-way agglomerative hierarchical clustering results show five main sample groups, S1 (red), S2 (yellow), S3 (blue), S4 (pink), and S5 (green), in addition to fetal brain control (light blue). Clinical and genomic covariates: Filled box (black), age = older than 1 year at diagnosis; INSS stage = 4; 1p36 = LOH; MYCN = amplified; 11q23 = LOH; 11q23 UNB = LOH; 17q = unbalanced gain. Gray, the condition is not known for that sample. Three main gene clusters are noted as G1 (blue), G2 (yellow), and G3 (red). Genes of interest are noted (with chromosome number in parentheses for those mapping to regions studied for DNA copy number as well), and those shown more than once reflect different probe sets for the same gene on the U95Av2 chip.

Figure 1.

Unsupervised hierarchical clustering of 101 primary neuroblastomas and human fetal brain. Unsupervised two-way agglomerative hierarchical clustering results show five main sample groups, S1 (red), S2 (yellow), S3 (blue), S4 (pink), and S5 (green), in addition to fetal brain control (light blue). Clinical and genomic covariates: Filled box (black), age = older than 1 year at diagnosis; INSS stage = 4; 1p36 = LOH; MYCN = amplified; 11q23 = LOH; 11q23 UNB = LOH; 17q = unbalanced gain. Gray, the condition is not known for that sample. Three main gene clusters are noted as G1 (blue), G2 (yellow), and G3 (red). Genes of interest are noted (with chromosome number in parentheses for those mapping to regions studied for DNA copy number as well), and those shown more than once reflect different probe sets for the same gene on the U95Av2 chip.

Close modal

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).

Figure 2.

Transcriptome maps. All genes present on the Affymetrix U95Av2 GeneChip located on chromosomes cytobands 1p11-1p36 (A), 2p24-25 (B), 11q11-11q25 (C), and 17q11-17q25 (D) are shown in linear order according to map position. Patient samples (columns) were grouped according to their genomic status at the locus of interest. Expression levels were first normalized and then binned into five levels (green, >1.5 SD below the mean; dark green, 1-1.5 SD below the mean; black, within 1 SD of the mean; dark red, 1-1.5 above the mean; red, >1.5 SD above the mean). Genes within or near a typical MYCN amplicon (37, 79) are annotated in (B).

Figure 2.

Transcriptome maps. All genes present on the Affymetrix U95Av2 GeneChip located on chromosomes cytobands 1p11-1p36 (A), 2p24-25 (B), 11q11-11q25 (C), and 17q11-17q25 (D) are shown in linear order according to map position. Patient samples (columns) were grouped according to their genomic status at the locus of interest. Expression levels were first normalized and then binned into five levels (green, >1.5 SD below the mean; dark green, 1-1.5 SD below the mean; black, within 1 SD of the mean; dark red, 1-1.5 above the mean; red, >1.5 SD above the mean). Genes within or near a typical MYCN amplicon (37, 79) are annotated in (B).

Close modal

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).

Table 2.

Region-specific alterations in gene expression

Region studied (Mb)No. differentially expressed genes-genome*No. genes in regionNo. genes on U95Av2No. differentially expressed genes-regionNo. validated differentially expressed genes-region
1p35-36 (0-34) 2,034 ↓ 563 201 122 ↓ 71 ↓ 
 2,352 ↑   5 ↑  
2p24-25 (0-24) 3,033 ↑ 130 47 37 ↑ 22 ↑ 
 3,236 ↓   3 ↓  
11q (53-134) 203 ↓ 1,079 329 88 ↓ Not done 
 90 ↑   3 ↑  
17q (22-80) 1,580 ↑ 931 372 150 ↑ Not done 
 1,780 ↓   58 ↓  
Region studied (Mb)No. differentially expressed genes-genome*No. genes in regionNo. genes on U95Av2No. differentially expressed genes-regionNo. validated differentially expressed genes-region
1p35-36 (0-34) 2,034 ↓ 563 201 122 ↓ 71 ↓ 
 2,352 ↑   5 ↑  
2p24-25 (0-24) 3,033 ↑ 130 47 37 ↑ 22 ↑ 
 3,236 ↓   3 ↓  
11q (53-134) 203 ↓ 1,079 329 88 ↓ Not done 
 90 ↑   3 ↑  
17q (22-80) 1,580 ↑ 931 372 150 ↑ Not done 
 1,780 ↓   58 ↓  
*

Number of differentially expressed genes (probe sets) with PaGE confidence levels ≥0.8 in the entire data set. Arrows indicate direction of change compared with the samples without the genomic alteration being analyzed.

Number of differentially expressed genes with PaGE confidence levels ≥0.8 in the region under consideration.

Number of differentially expressed genes with PaGE confidence levels ≥0.8 in the validation data set that were also present in the study data set for each region under consideration.

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.

Figure 3.

Quantitative RT-PCR validation of oligonucleotide-based gene expression. Representative results for five genes located within genomic regions of interest. Log-transformed data were normalized to GAPDH and represented as box and whisker plots for both microarray (A) and RT-PCR (B) experiments. Statistical significance was determined with the Wilcoxon's test statistic.

Figure 3.

Quantitative RT-PCR validation of oligonucleotide-based gene expression. Representative results for five genes located within genomic regions of interest. Log-transformed data were normalized to GAPDH and represented as box and whisker plots for both microarray (A) and RT-PCR (B) experiments. Statistical significance was determined with the Wilcoxon's test statistic.

Close modal

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).

Figure 4.

Heat map representation of most differentially expressed genes at 1p35-36 (A) and 11q (B). Top 50 genes discovered to be differentially expressed in each binary comparison of LOH versus no LOH for chromosome bands 1p36 and 11q23. Genes are sorted by their PaGE generated t statistic (available at http://stokes.chop.edu/programs/marislab/data).

Figure 4.

Heat map representation of most differentially expressed genes at 1p35-36 (A) and 11q (B). Top 50 genes discovered to be differentially expressed in each binary comparison of LOH versus no LOH for chromosome bands 1p36 and 11q23. Genes are sorted by their PaGE generated t statistic (available at http://stokes.chop.edu/programs/marislab/data).

Close modal

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.

Table 3.

Genes mapping to distal 1p with significantly lower expression in both test and validation data sets

Gene*Locust statistics
Confidence
PennMSKCCPennMSKCC
STX12 1p35.3 −9.09 −2.66 0.996 0.980 
PTP4A2 1p35 −7.62 −1.66 0.996 0.894 
GNB1 1p36.33 −7.21 −1.75 0.996 0.904 
CLSTN1 1p36.22 −7.05 −2.47 0.996 0.973 
PUM1 1p35.2 −6.73 −1.83 0.996 0.915 
RERE 1p36.23 −6.32 −2.48 0.996 0.973 
DNAJC8 1p35.2 −6.22 −2.51 0.996 0.974 
CDC42 1p36.1 −6.22 −1.57 0.996 0.889 
KIAA0962 1p36.1 −5.95 −1.77 0.996 0.909 
MFN2 1p36.21 −5.90 −2.06 0.996 0.942 
KIAA0842 1p36.13 −5.90 −2.39 0.996 0.971 
BAI2 1p35 −5.87 −2.54 0.996 0.974 
ICMT 1p36.21 −5.72 −1.12 0.996 0.835 
PRDM2 1p36 −5.65 −1.47 0.996 0.889 
SLC35E2 1p36.33 −5.64 −2.69 0.996 0.980 
P29 1p36.11 −5.52 −1.90 0.996 0.920 
PIK3CD 1p36.2 −5.51 −1.99 0.996 0.938 
CHD5 1p36.23 −5.45 −1.07 0.996 0.824 
CAMTA1 1p36.23 −5.19 −1.99 0.996 0.938 
CLCN6 1p36 −5.04 −1.44 0.996 0.889 
FRAP1 1p36.2 −5.02 −0.99 0.996 0.805 
SDC3 1p35.2 −4.94 −3.23 0.996 0.983 
PMSCL2 1p36.22 −4.88 −2.22 0.996 0.957 
CTNNBIP1 1p36.22 −4.86 −2.73 0.996 0.980 
CDC2L1 1p36 −4.79 −1.53 0.996 0.889 
RER1 1p36.33 −4.78 −1.57 0.996 0.889 
KIAA0453 1p36.21 −4.77 −2.03 0.996 0.940 
DJ159A19.3 1p36.13 −4.70 −1.77 0.996 0.909 
SRRM1 1p36.11 −4.60 −3.29 0.996 0.985 
SFRS4 1p35.2 −4.51 −2.66 0.996 0.980 
UBE4B 1p36.3 −4.49 −1.69 0.996 0.895 
CAPZB 1p36.1 −4.47 −3.10 0.996 0.980 
TCEB3 1p36.1 −4.47 −1.54 0.996 0.889 
EIF4G3 1p36.12 −4.39 −1.77 0.996 0.909 
VAMP3 1p36.23 −4.35 −2.03 0.996 0.940 
CASP9 1p36.21 −4.35 −0.97 0.996 0.805 
FUCA1 1p36.11 −4.10 −2.02 0.996 0.940 
RER1 1p36.33 −4.07 −1.30 0.996 0.866 
AKR7A2 1p36.13 −3.86 −1.51 0.996 0.889 
ECE1 1p36.1 −3.80 −2.09 0.996 0.945 
CDC2L2 1p36.3 −3.73 −1.99 0.996 0.938 
CDC2L2 1p36.3 −3.73 −2.54 0.996 0.974 
DJ465N24.2.1 1p36.11 −3.60 −3.57 0.995 0.985 
MGC33488 1p36.23 −3.51 −0.97 0.995 0.805 
NUDC 1p36.11 −3.49 −1.44 0.995 0.889 
PTP4A2 1p35 −3.48 −2.41 0.995 0.971 
SHARP 1p36.13 −3.38 −1.53 0.995 0.889 
KIAA0495 1p36.32 −3.29 −1.31 0.995 0.867 
KIAA0469 1p36.23 −3.26 −2.40 0.994 0.971 
VAMP3 1p36.23 −3.13 −1.08 0.994 0.824 
KIAA0450 1p36.32 −3.12 −1.44 0.994 0.889 
TNFRSF25 1p36.2 −3.12 −2.72 0.994 0.980 
PRDM2 1p36 −3.09 −1.88 0.994 0.918 
TNFRSF1B 1p36.22 −2.92 −2.67 0.992 0.980 
RBAF600 1p36.13 −2.91 −1.47 0.992 0.889 
TNFRSF1B 1p36.22 −2.67 −2.69 0.988 0.980 
RUNX3 1p36 −2.61 −1.54 0.987 0.889 
HMGCL 1p36.11 −2.60 −1.07 0.987 0.824 
ARHGEF16 1p36.3 −2.55 −1.31 0.986 0.867 
CASP9 1p36.21 −2.49 −1.54 0.985 0.889 
KIAA1037 1p36.11 −2.47 −2.96 0.984 0.980 
HTR6 1p36.13 −2.44 −1.57 0.983 0.889 
STMN1 1p36.11 −2.34 −1.66 0.98 0.894 
LYPLA2 1p36.11 −2.22 −1.28 0.973 0.866 
TAF12 1p35.2 −2.12 −1.54 0.967 0.889 
ZNF151 1p36.13 −2.07 −2.72 0.963 0.980 
KIAA0478 1p36.12 −1.98 −1.39 0.958 0.881 
LAPTM5 1p35.2 −1.89 −2.90 0.951 0.980 
SCNN1D 1p36.33 −1.84 −1.23 0.947 0.858 
G1P2 1p36.33 −1.73 −1.45 0.935 0.889 
G1P2 1p36.33 −1.66 −1.33 0.926 0.868 
CLIC4 1p36.11 −1.57 −1.24 0.918 0.859 
RAP1GA1 1p36.12 −1.47 −1.37 0.904 0.879 
CDW52 1p36 −1.47 −1.35 0.904 0.873 
CORT 1p36.22 −1.30 −1.63 0.874 0.889 
SDR1 1p36.1 −1.18 −1.50 0.85 0.889 
LCK 1p35.2 −1.15 −1.41 0.845 0.882 
HMGN2 1p36.1 −1.06 −1.84 0.835 0.915 
GJB5 1p35.1 −0.99 −1.52 0.821 0.889 
Gene*Locust statistics
Confidence
PennMSKCCPennMSKCC
STX12 1p35.3 −9.09 −2.66 0.996 0.980 
PTP4A2 1p35 −7.62 −1.66 0.996 0.894 
GNB1 1p36.33 −7.21 −1.75 0.996 0.904 
CLSTN1 1p36.22 −7.05 −2.47 0.996 0.973 
PUM1 1p35.2 −6.73 −1.83 0.996 0.915 
RERE 1p36.23 −6.32 −2.48 0.996 0.973 
DNAJC8 1p35.2 −6.22 −2.51 0.996 0.974 
CDC42 1p36.1 −6.22 −1.57 0.996 0.889 
KIAA0962 1p36.1 −5.95 −1.77 0.996 0.909 
MFN2 1p36.21 −5.90 −2.06 0.996 0.942 
KIAA0842 1p36.13 −5.90 −2.39 0.996 0.971 
BAI2 1p35 −5.87 −2.54 0.996 0.974 
ICMT 1p36.21 −5.72 −1.12 0.996 0.835 
PRDM2 1p36 −5.65 −1.47 0.996 0.889 
SLC35E2 1p36.33 −5.64 −2.69 0.996 0.980 
P29 1p36.11 −5.52 −1.90 0.996 0.920 
PIK3CD 1p36.2 −5.51 −1.99 0.996 0.938 
CHD5 1p36.23 −5.45 −1.07 0.996 0.824 
CAMTA1 1p36.23 −5.19 −1.99 0.996 0.938 
CLCN6 1p36 −5.04 −1.44 0.996 0.889 
FRAP1 1p36.2 −5.02 −0.99 0.996 0.805 
SDC3 1p35.2 −4.94 −3.23 0.996 0.983 
PMSCL2 1p36.22 −4.88 −2.22 0.996 0.957 
CTNNBIP1 1p36.22 −4.86 −2.73 0.996 0.980 
CDC2L1 1p36 −4.79 −1.53 0.996 0.889 
RER1 1p36.33 −4.78 −1.57 0.996 0.889 
KIAA0453 1p36.21 −4.77 −2.03 0.996 0.940 
DJ159A19.3 1p36.13 −4.70 −1.77 0.996 0.909 
SRRM1 1p36.11 −4.60 −3.29 0.996 0.985 
SFRS4 1p35.2 −4.51 −2.66 0.996 0.980 
UBE4B 1p36.3 −4.49 −1.69 0.996 0.895 
CAPZB 1p36.1 −4.47 −3.10 0.996 0.980 
TCEB3 1p36.1 −4.47 −1.54 0.996 0.889 
EIF4G3 1p36.12 −4.39 −1.77 0.996 0.909 
VAMP3 1p36.23 −4.35 −2.03 0.996 0.940 
CASP9 1p36.21 −4.35 −0.97 0.996 0.805 
FUCA1 1p36.11 −4.10 −2.02 0.996 0.940 
RER1 1p36.33 −4.07 −1.30 0.996 0.866 
AKR7A2 1p36.13 −3.86 −1.51 0.996 0.889 
ECE1 1p36.1 −3.80 −2.09 0.996 0.945 
CDC2L2 1p36.3 −3.73 −1.99 0.996 0.938 
CDC2L2 1p36.3 −3.73 −2.54 0.996 0.974 
DJ465N24.2.1 1p36.11 −3.60 −3.57 0.995 0.985 
MGC33488 1p36.23 −3.51 −0.97 0.995 0.805 
NUDC 1p36.11 −3.49 −1.44 0.995 0.889 
PTP4A2 1p35 −3.48 −2.41 0.995 0.971 
SHARP 1p36.13 −3.38 −1.53 0.995 0.889 
KIAA0495 1p36.32 −3.29 −1.31 0.995 0.867 
KIAA0469 1p36.23 −3.26 −2.40 0.994 0.971 
VAMP3 1p36.23 −3.13 −1.08 0.994 0.824 
KIAA0450 1p36.32 −3.12 −1.44 0.994 0.889 
TNFRSF25 1p36.2 −3.12 −2.72 0.994 0.980 
PRDM2 1p36 −3.09 −1.88 0.994 0.918 
TNFRSF1B 1p36.22 −2.92 −2.67 0.992 0.980 
RBAF600 1p36.13 −2.91 −1.47 0.992 0.889 
TNFRSF1B 1p36.22 −2.67 −2.69 0.988 0.980 
RUNX3 1p36 −2.61 −1.54 0.987 0.889 
HMGCL 1p36.11 −2.60 −1.07 0.987 0.824 
ARHGEF16 1p36.3 −2.55 −1.31 0.986 0.867 
CASP9 1p36.21 −2.49 −1.54 0.985 0.889 
KIAA1037 1p36.11 −2.47 −2.96 0.984 0.980 
HTR6 1p36.13 −2.44 −1.57 0.983 0.889 
STMN1 1p36.11 −2.34 −1.66 0.98 0.894 
LYPLA2 1p36.11 −2.22 −1.28 0.973 0.866 
TAF12 1p35.2 −2.12 −1.54 0.967 0.889 
ZNF151 1p36.13 −2.07 −2.72 0.963 0.980 
KIAA0478 1p36.12 −1.98 −1.39 0.958 0.881 
LAPTM5 1p35.2 −1.89 −2.90 0.951 0.980 
SCNN1D 1p36.33 −1.84 −1.23 0.947 0.858 
G1P2 1p36.33 −1.73 −1.45 0.935 0.889 
G1P2 1p36.33 −1.66 −1.33 0.926 0.868 
CLIC4 1p36.11 −1.57 −1.24 0.918 0.859 
RAP1GA1 1p36.12 −1.47 −1.37 0.904 0.879 
CDW52 1p36 −1.47 −1.35 0.904 0.873 
CORT 1p36.22 −1.30 −1.63 0.874 0.889 
SDR1 1p36.1 −1.18 −1.50 0.85 0.889 
LCK 1p35.2 −1.15 −1.41 0.845 0.882 
HMGN2 1p36.1 −1.06 −1.84 0.835 0.915 
GJB5 1p35.1 −0.99 −1.52 0.821 0.889 
*

Of the 146 probe sets representing 122 genes with confidence levels >0.8 in the Penn data set, 111 probe sets representing 98 genes were also shown to be significantly differentially expressed with confidence levels >0.8 in the MSKCC data set. Genes listed more than once in cases where more than one probe set for the given gene on the Affymetrix U95Av2 array showed differential expression. Genes in boldface have been postulated to be neuroblastoma suppressor genes previously (see text) and/or were also found to be differentially expressed based on 1p LOH status by Janoueix-Lerosey et al. (48) or Caren et al. (49).

The t statistic and confidence levels as computed with the PaGE algorithm. Negative numbers denote decreased expression in samples with 1p LOH.

Figure 5.

Heat map representation of most differentially expressed genes at 2p24-25 in the Penn (A) and MSKCC (B) data sets. Genes are sorted by their PaGE generated t statistic (available at http://stokes.chop.edu/programs/marislab/data) and show complete concordance for the genes most significantly differentially expressed between MYCN-amplified and single-copy cases.

Figure 5.

Heat map representation of most differentially expressed genes at 2p24-25 in the Penn (A) and MSKCC (B) data sets. Genes are sorted by their PaGE generated t statistic (available at http://stokes.chop.edu/programs/marislab/data) and show complete concordance for the genes most significantly differentially expressed between MYCN-amplified and single-copy cases.

Close modal

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 (3436). 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 (4042). 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 (4547). 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, 5456), 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 (4042) 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.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

Presented in part at the 96th Annual Meeting of the AACR, Anaheim, CA, April 16-20, 2005 and the 11th Conference of the Advances in Neuroblastoma, Genoa, Italy, June 16-19, 2004.

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.

1
Brodeur GM, Maris JM. Neuroblastoma. In: Pizzo PA, Poplack DG, editors. Principles and practice of pediatric oncology. 4th ed. Philadelphia: J.B. Lippincott Co.; 2002. p.895–938.
2
Maris JM. The biologic basis for neuroblastoma heterogeneity and risk stratification.
Curr Opin Pediatr
2005
;
17
:
7
–13.
3
Breslow N, McCann B. Statistical estimation of prognosis for children with neuroblastoma.
Cancer Res
1971
;
31
:
2098
–103.
4
Brodeur GM, Pritchard J, Berthold F, et al. Revisions of the international criteria for neuroblastoma diagnosis, staging, and response to treatment.
J Clin Oncol
1993
;
11
:
1466
–77.
5
Shimada H, Ambros IM, Dehner LP, et al. The International Neuroblastoma Pathology Classification (the Shimada system).
Cancer
1999
;
86
:
364
–72.
6
Look AT, Hayes FA, Shuster JJ, Douglass EC, Castleberry RP, Brodeur GM. Clinical relevance of tumor cell ploidy and N-myc gene amplification in childhood neuroblastoma: a Pediatric Oncology Group study.
J Clin Oncol
1991
;
9
:
581
–91.
7
Brodeur GM, Seeger RC, Schwab M, Varmus HE, Bishop JM. Amplification of N-myc in untreated human neuroblastomas correlates with advanced disease stage.
Science
1984
;
224
:
1121
–4.
8
Seeger RC, Brodeur GM, Sather H, et al. Association of multiple copies of the N-myc oncogene with rapid progression of neuroblastomas.
N Engl J Med
1985
;
313
:
1111
–6.
9
Brodeur GM. Neuroblastoma: biological insights into a clinical enigma.
Nat Rev Cancer
2003
;
3
:
203
–16.
10
Schmidt ML, Lukens JN, Seeger RC, et al. Biologic factors determine prognosis in infants with stage IV neuroblastoma: a prospective Children's Cancer Group study.
J Clin Oncol
2000
;
18
:
1260
–8.
11
Cohn SL, Tweddle DA. MYCN amplification remains prognostically strong 20 years after its “clinical debut.”
Eur J Cancer
2004
;
40
:
2639
–42.
12
Maris JM, Matthay KK. Molecular biology of neuroblastoma.
J Clin Oncol
1999
;
17
:
2264
–79.
13
Attiyeh EF, London WB, Mosse YP, et al. Chromosome 1p and 11q deletions and outcome in neuroblastoma.
N Engl J Med
2005
;
353
:
2243
–53.
14
Nigro JM, Misra A, Zhang L, et al. Integrated array-comparative genomic hybridization and expression array profiles identify clinically relevant molecular subtypes of glioblastoma.
Cancer Res
2005
;
65
:
1678
–86.
15
Yi Y, Mirosevich J, Shyr Y, Matusik R, George AL, Jr. Coupled analysis of gene expression and chromosomal location.
Genomics
2005
;
85
:
401
–12.
16
Alaminos M, Mora J, Cheung NK, et al. Genome-wide analysis of gene expression associated with MYCN in human neuroblastoma.
Cancer Res
2003
;
63
:
4538
–46.
17
Wei JS, Greer BT, Westermann F, et al. Prediction of clinical outcome using gene expression profiling and artificial neural networks for patients with neuroblastoma.
Cancer Res
2004
;
64
:
6883
–91.
18
Ohira M, Oba S, Nakamura Y, Hirata T, Ishii S, Nakagawara A. A review of DNA microarray analysis of human neuroblastomas.
Cancer Lett
2005
;
228
:
5
–11.
19
McArdle L, McDermott M, Purcell R, et al. Oligonucleotide microarray analysis of gene expression in neuroblastoma displaying loss of chromosome 11q.
Carcinogenesis
2004
;
25
:
1599
–609.
20
Mathew P, Valentine MB, Bowman LC, et al. Detection of MYCN gene amplification in neuroblastoma by fluorescence in situ hybridization: a pediatric oncology group study.
Neoplasia
2001
;
3
:
105
–9.
21
Morowitz M, Shusterman S, Mosse Y, et al. Detection of single-copy chromosome 17q gain in human neuroblastomas using real-time quantitative polymerase chain reaction.
Mod Pathol
2003
;
16
:
1248
–56.
22
Alaminos M, Gerald WL, Cheung NK. Prognostic value of MYCN and ID2 overexpression in neuroblastoma.
Pediatr Blood Cancer
2005
;
45
:
909
–15.
23
Ashburner M, Ball CA, Blake JA, et al. Gene Ontology: tool for the unification of biology. The Gene Ontology Consortium.
Nat Genet
2000
;
25
:
25
–9.
24
Hosack DA, Dennis G, Jr., Sherman BT, Lane HC, Lempicki RA. Identifying biological themes within lists of genes with EASE.
Genome Biol
2003
;
4
:
R70
.
25
Manduchi E, Grant GR, McKenzie SE, Overton GC, Surrey S, Stoeckert CJ, Jr. Generation of patterns from gene expression data by assigning confidence to differentially expressed genes.
Bioinformatics
2000
;
16
:
685
–98.
26
Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response.
Proc Natl Acad Sci U S A
2001
;
98
:
5116
–21.
27
Kent WJ, Sugnet CW, Furey TS, et al. The human genome browser at UCSC.
Genome Res
2002
;
12
:
996
–1006.
28
Wang Q, Hii G, Shusterman S, et al. ID2 expression is not associated with MYCN amplification or expression in human neuroblastomas.
Cancer Res
2003
;
63
:
1631
–5.
29
Vandesompele J, Baudis M, De Preter K, et al. Unequivocal delineation of clinicogenetic subgroups and development of a new model for improved outcome prediction in neuroblastoma.
J Clin Oncol
2005
;
23
:
2280
–99.
30
Shohet JM, Hicks MJ, Plon SE, et al. Minichromosome maintenance protein MCM7 is a direct target of the MYCN transcription factor in neuroblastoma.
Cancer Res
2002
;
62
:
1123
–8.
31
Breen CJ, O'Meara A, McDermott M, Mullarkey M, Stallings RL. Coordinate deletion of chromosome 3p and 11q in neuroblastoma detected by comparative genomic hybridization.
Cancer Genet Cytogenet
2000
;
120
:
44
–9.
32
Golub TR, Slonim DK, Tamayo P, et al. Molecular classification of cancer: class discovery and class prediction by gene expression monitoring.
Science
1999
;
286
:
531
–7.
33
Khan J, Wei JS, Ringner M, et al. Classification and diagnostic prediction of cancers using gene expression profiling and artificial neural networks.
Nat Med
2001
;
7
:
673
–9.
34
Yang XJ, Tan MH, Kim HL, et al. A molecular classification of papillary renal cell carcinoma.
Cancer Res
2005
;
65
:
5628
–37.
35
Haven CJ, Howell VM, Eilers PH, et al. Gene expression of parathyroid tumors: molecular subclassification and identification of the potential malignant phenotype.
Cancer Res
2004
;
64
:
7405
–11.
36
Onken MD, Worley LA, Ehlers JP, Harbour JW. Gene expression profiling in uveal melanoma reveals two molecular classes and predicts metastatic death.
Cancer Res
2004
;
64
:
7205
–9.
37
Reiter JL, Brodeur GM. High-resolution mapping of a 130-kb core region of the MYCN amplicon in neuroblastomas.
Genomics
1996
;
32
:
97
–103.
38
Vasudevan SA, Nuchtern JG. Gene profiling of high risk neuroblastoma.
World J Surg
2005
;
29
:
317
–24.
39
Bown N, Cotterill S, Lastowska M, et al. Gain of chromosome arm 17q and adverse outcome in patients with neuroblastoma.
N Engl J Med
1999
;
340
:
1954
–61.
40
Chang CL, Zhu XX, Thoraval DH, et al. Nm23-1 mutation in neuroblastoma.
Nature
1994
;
370
:
335
–6.
41
Valentijn LJ, Koppen A, van Asperen R, Root HA, Haneveld F, Versteeg R. Inhibition of a new differentiation pathway in neuroblastoma by copy number defects of N-myc, Cdc42, and nm23 genes.
Cancer Res
2005
;
65
:
3136
–45.
42
Saito-Ohara F, Imoto I, Inoue J, et al. PPM1D is a potential target for 17q gain in neuroblastoma.
Cancer Res
2003
;
63
:
1876
–83.
43
Maris JM, Weiss MJ, Guo C, et al. Loss of heterozygosity at 1p36 independently predicts for disease progression, but not decreased overall survival probability, in neuroblastoma patients: a Children's Cancer Group study.
J Clin Oncol
2000
;
18
:
1888
–99.
44
Matthay KK, Villablanca JG, Seeger RC, et al. Treatment of high-risk neuroblastoma with intensive chemotherapy, radiotherapy, autologous bone marrow transplantation, and 13-cis-retinoic acid. Children's Cancer Group.
N Engl J Med
1999
;
341
:
1165
–73.
45
Cooper R, Khakoo Y, Matthay KK, et al. Opsoclonus-myoclonus-ataxia syndrome in neuroblastoma: histopathologic features—a report from the Children's Cancer Group.
Med Pediatr Oncol
2001
;
36
:
623
–9.
46
Gambini C, Conte M, Bernini G, et al. Neuroblastic tumors associated with opsoclonus-myoclonus syndrome: histological, immunohistochemical and molecular features of 15 Italian cases.
Virchows Arch
2003
;
442
:
555
–62.
47
Shimada H. International neuroblastoma pathology classification for prognostic evaluation of patients with peripheral neuroblastic tumors: a report from the Children's Cancer Group.
Cancer
2001
;
92
:
2451
–61.
48
Janoueix-Lerosey I, Novikov E, Monteiro M, et al. Gene expression profiling of 1p35-36 genes in neuroblastoma.
Oncogene
2004
;
23
:
5912
–22.
49
Caren H, Ejeskar K, Fransson S, et al. A cluster of genes located in 1p36 are down-regulated in neuroblastomas with poor prognosis, but not due to CpG island methylation.
Mol Cancer
2005
;
4
:
10
.
50
Teitz T, Wei T, Liu D, et al. Caspase-9 and Apaf-1 are expressed and functionally active in human neuroblastoma tumor cell lines with 1p36 LOH and amplified MYCN.
Oncogene
2002
;
21
:
1848
–58.
51
Abel F, Sjoberg RM, Nilsson S, Kogner P, Martinsson T. Imbalance of the mitochondrial pro- and anti-apoptotic mediators in neuroblastoma tumours with unfavourable biology.
Eur J Cancer
2005
;
41
:
635
–46.
52
Barbashina V, Salazar P, Holland EC, Rosenblum MK, Ladanyi M. Allelic losses at 1p36 and 19q13 in gliomas: correlation with histologic classification, definition of a 150-kb minimal deleted region on 1p36, and evaluation of CAMTA1 as a candidate tumor suppressor gene.
Clin Cancer Res
2005
;
11
:
1119
–28.
53
White PS, Thompson PM, Gotoh T, et al. Definition and characterization of a region of 1p36.3 consistently deleted in neuroblastoma.
Oncogene
2005
;
24
:
2684
–94.
54
Lahti JM, Valentine M, Xiang J, et al. Alterations in the PITSLRE protein kinase gene complex on chromosome 1p36 in childhood neuroblastoma.
Nat Genet
1994
;
7
:
370
–5.
55
He L, Yu JX, Liu L, et al. RIZ1, but not the alternative RIZ2 product of the same gene, is underexpressed in breast cancer, and forced RIZ1 expression causes G2-M cell cycle arrest and/or apoptosis.
Cancer Res
1998
;
58
:
4238
–44.
56
Thompson PM, Gotoh T, Kok M, White PS, Brodeur GM. CHD5, a new member of the chromodomain gene family, is preferentially expressed in the nervous system.
Oncogene
2003
;
22
:
1002
–11.
57
Bilke S, Chen QR, Westerman F, Schwab M, Catchpoole D, Khan J. Inferring a tumor progression model for neuroblastoma from genomic data.
J Clin Oncol
2005
;
23
:
7322
–31.
58
Selzer RR, Richmond TA, Pofahl NJ, et al. Analysis of chromosome breakpoints in neuroblastoma at sub-kilobase resolution using fine-tiling oligonucleotide array CGH.
Genes Chromosomes Cancer
2005
;
44
:
305
–19.
59
Mosse Y, Greshock J, King A, Khazi D, Weber BL, Maris JM. Identification and high-resolution mapping of a constitutional 11q deletion in an infant with multifocal neuroblastoma.
Lancet Oncol
2003
;
4
:
769
–71.
60
Satge D, Moore SW, Stiller CA, et al. Abnormal constitutional karyotypes in patients with neuroblastoma: a report of four new cases and review of 47 others in the literature.
Cancer Genet Cytogenet
2003
;
147
:
89
–98.
61
Guo C, White PS, Weiss MJ, et al. Allelic deletion at 11q23 is common in MYCN single copy neuroblastomas.
Oncogene
1999
;
18
:
4948
–57.
62
De Preter K, Vandesompele J, Menten B, et al. Positional and functional mapping of a neuroblastoma differentiation gene on chromosome 11.
BMC Genomics
2005
;
6
:
97
.
63
Kuramochi M, Fukuhara H, Nobukuni T, et al. TSLC1 is a tumor-suppressor gene in human non-small-cell lung cancer.
Nat Genet
2001
;
27
:
427
–30.
64
Neumann HP, Bausch B, McWhinney SR, et al. Germ-line mutations in nonsyndromic pheochromocytoma.
N Engl J Med
2002
;
346
:
1459
–66.
65
Milne TA, Hughes CM, Lloyd R, et al. Menin and MLL cooperatively regulate expression of cyclin-dependent kinase inhibitors.
Proc Natl Acad Sci U S A
2005
;
102
:
749
–54.
66
Wang SS, Esplin ED, Li JL, et al. Alterations of the PPP2R1B gene in human lung and colon cancer.
Science
1998
;
282
:
284
–7.
67
Bartkova J, Horejsi Z, Koed K, et al. DNA damage response as a candidate anti-cancer barrier in early human tumorigenesis.
Nature
2005
;
434
:
864
–70.
68
Blaheta RA, Hundemer M, Mayer G, et al. Expression level of neural cell adhesion molecule (NCAM) inversely correlates with the ability of neuroblastoma cells to adhere to endothelium in vitro.
Cell Commun Adhes
2002
;
9
:
131
–47.
69
Collado M, Gil J, Efeyan A, et al. Tumour biology: senescence in premalignant tumours.
Nature
2005
;
436
:
642
.
70
De Preter K, Vandesompele J, Hoebeeck J, et al. No evidence for involvement of SDHD in neuroblastoma pathogenesis.
BMC Cancer
2004
;
4
:
55
.
71
Teitz T, Wei T, Valentine MB, et al. Caspase 8 is deleted or silenced preferentially in childhood neuroblastomas with amplification of MYCN.
Nat Med
2000
;
6
:
529
–35.
72
Yang Q, Zage P, Kagan D, et al. Association of epigenetic inactivation of RASSF1A with poor outcome in human neuroblastoma.
Clin Cancer Res
2004
;
10
:
8493
–500.
73
Yang QW, Liu S, Tian Y, et al. Methylation-associated silencing of the thrombospondin-1 gene in human neuroblastoma.
Cancer Res
2003
;
63
:
6299
–310.
74
Misawa A, Inoue J, Sugino Y, et al. Methylation-associated silencing of the nuclear receptor 1I2 gene in advanced-type neuroblastomas, identified by bacterial artificial chromosome array-based methylated CpG island amplification.
Cancer Res
2005
;
65
:
10233
–42.
75
Banelli B, Gelvi I, Di Vinci A, et al. Distinct CpG methylation profiles characterize different clinical groups of neuroblastic tumors.
Oncogene
2005
;
24
:
5619
–28.
76
Abe M, Ohira M, Kaneda A, et al. CpG island methylator phenotype is a strong determinant of poor prognosis in neuroblastomas.
Cancer Res
2005
;
65
:
828
–34.
77
Mosse YP, Greshock J, Margolin A, et al. High-resolution detection and mapping of genomic DNA alterations in neuroblastoma.
Genes Chromosomes Cancer
2005
;
43
:
390
–403.
78
Manning G, Whyte DB, Martinez R, Hunter T, Sudarsanam S. The protein kinase complement of the human genome.
Science
2002
;
298
:
1912
–34.
79
Reiter JL, Brodeur GM. MYCN is the only highly expressed gene from the core amplified domain in human neuroblastomas.
Genes Chromosomes Cancer
1998
;
23
:
134
–40.

Supplementary data