Purpose: The purpose of this study was to further define the biology of the 11q− neuroblastoma tumor subgroup by the integration of array-based comparative genomic hybridization with microRNA (miRNA) expression profiling data to determine if improved patient stratification is possible.
Experimental Design: A set of primary neuroblastoma (n = 160), which was broadly representative of all genetic subtypes, was analyzed by array-based comparative genomic hybridization and for the expression of 430 miRNAs. A 15-miRNA expression signature previously shown to be predictive of clinical outcome was used to analyze an independent cohort of 11q− tumors (n = 37).
Results: Loss of 4p and gain of 7q occurred at a significantly higher frequency in the 11q− tumors, further defining the genetic characteristics of this subtype. The 11q− tumors could be split into two subgroups using a miRNA expression survival signature that differed significantly in clinical outcome and the overall frequency of large-scale genomic imbalances, with the poor survival subgroup having significantly more imbalances. miRNAs from the expression signature, which were upregulated in unfavorable tumors, were predicted to target downregulated genes from a published mRNA expression classifier of clinical outcome at a higher-than-expected frequency, indicating the miRNAs might contribute to the regulation of genes within the signature.
Conclusion: We show that two distinct biological subtypes of neuroblastoma with loss of 11q occur, which differ in their miRNA expression profiles, frequency of segmental imbalances, and clinical outcome. A miRNA expression signature, combined with an analysis of segmental imbalances, provides greater prediction of event-free survival and overall survival outcomes than 11q status by itself, improving patient stratification. Clin Cancer Res; 16(11); 2971–8. ©2010 AACR.
Loss of chromosome 11q material in neuroblastoma is generally associated with poor clinical outcome. The results presented in this report, however, show that using miRNA expression profiling, the 11q− genetic subgroup can be subdivided into two distinct groups. These subgroups differ significantly in clinical outcome and with respect to the overall frequency of segmental genomic imbalances, indicating the possibility for more refined patient stratification before therapy.
Neuroblastoma, one of the most common solid tumors in children, is derived from precursor cells of the sympathetic nervous system and is responsible for approximately 15% of all childhood cancer deaths (see ref. 1 for review). These tumors are particularly noted for extensive heterogeneity in clinical behavior, ranging from spontaneous regression to aggressive clinical course and death due to disease. Tumors in infants characterized by whole chromosome gains and losses and few structural chromosome abnormalities usually have favorable clinical outcomes and, in some instances, can even spontaneously regress. In contrast, tumors from children >1.5 years of age are often metastatic at diagnosis and become refractory to treatment. These unfavorable tumor subtypes possess a broad spectrum of recurrent chromosomal gains, losses, and amplifications, many of which are highly associated with clinical outcome and form the basis for dividing neuroblastoma into genetic subtypes.
Tumors with amplification of the MYCN transcription factor represent a major genetic subtype of metastatic neuroblastoma with a particularly poor prognosis (2). Loss of chromosome 1p and gain of chromosome 17q regions are common in MYCN-amplified (MNA) tumors, often occurring by a common chromosomal mechanism involving an unbalanced t(1p;17q) (3). Tumors with hemizygous loss of a large segment on chromosome 11q represent another major genetic subtype of metastatic neuroblastoma (4, 5), which can also be associated with unfavorable clinical outcome (6, 7). Loss of chromosome 3p occurs preferentially in 11q− tumors, and interestingly, children with 3p−/11q− tumors have a higher median age at diagnosis than do those with MNA tumors (6, 8, 9). Similar to MNA tumors, the 11q− tumors usually have gain of 17q, but as a consequence of an unbalanced t(11q;17q), which simultaneously results in loss of the 11q material (10, 11). At the molecular level, MNA and 11q− tumors are clearly different, being distinguishable by both mRNA (12, 13) and microRNA (miRNA; ref. 14) genome-wide expression profiling. There are also many additional recurring large-scale chromosome imbalances found in neuroblastomas, including loss of 4p, 9p, and 14q and gain of 1q, 7q, 2p, and 11p (see ref. 15 for review). Approximately 70% to 80% of all metastatic (stage IV) neuroblastomas have neither MNA nor loss of 11q material and form a much more genetically heterogeneous group of tumors (15–17).
Although the 11q− subtype has been shown to be associated with poor clinical outcome (6, 7, 15), exactly what other genetic factors might mediate patient survival within this subgroup has remained poorly defined. Here, we have carried out an integrated analysis of DNA copy number and miRNA expressional profiling on a set of primary neuroblastomas to further define the clinical and genetic characteristics of this important tumor subtype.
Materials and Methods
Neuroblastoma tumor samples were obtained retrospectively from Our Lady's Children's Hospital, Dublin, Ireland (n = 49), and the Children's Oncology Group Neuroblastoma Tumor Bank (n = 111), for a total of 160 tumors [array-based comparative genomic hybridization (aCGH) profiling was carried on the entire set of tumors, whereas miRNA profiling is available on a subset (n = 144)]. Additional 11q− tumors were also obtained retrospectively from University Children's Hospital, Essen, Germany (n = 12), and Ghent University Hospital, Ghent, Belgium (n = 11), for validation purposes. These tumors were also characterized by aCGH and by miRNA expression profiling. Patients were treated under either the European treatment protocol or the U.S. neuroblastoma treatment protocol between 1998 and 2004. Clinical details are summarized in Supplementary Table S1 and Fig. S1. Event-free survival (EFS) time was calculated from the time of enrollment on the front-line or biological study until the time of the first occurrence of relapse, progressive disease, secondary malignancy, or death, or until the time of last contact if no event occurred. Overall survival (OS) time was calculated until the time of death or until last contact. EFS and OS are presented as the estimate ± SE.
miRNA expression profiling
For tumors from Dublin and Children's Oncology Group, DNA was extracted from tumor samples using the DNeasy Kit (Qiagen). DNA was quantified using the Nanodrop spectrophotometer and stored at −80°C. Methods for array hybridization and analysis were as previously reported (20) and summarized in NimbleGen aCGH protocol v5.1, with minor modifications. Briefly, a total of 500 ng of tumor DNA and control DNA were differentially labeled with Cy3 and Cy5 (TriLink Biotechnologies), respectively. Four micrograms of labeled test and reference were hybridized to a 4-plex 72,000 probe array (Roche NimbleGen) for 18 hours in a MAUI hybridization station (BioMicro Systems) at 42°C. Posthybridization washes were done and arrays scanned using the Axon 4000B scanner and GenePix 6.0 (Molecular Devices). Fluorescent intensities were extracted and log 2 ratios calculated and qspline normalized (21) using NimbleScan Software version 2.4. Aberration breakpoints were identified using the CGH-segMNT algorithm in NimbleScan 2.4. Unaveraged log 2 ratios were written out as .gff files and used for data visualization. In general, a significance log 2 ratio threshold of <−0.25 for loss and >0.25 for gain was used to identify DNA copy number imbalances. For tumors from Essen and Ghent, a custom Agilent 44K array was applied for aCGH using standard protocols.
Unsupervised computational analysis
Unsupervised partitional clustering (k-means) was used to optimize a k = 2 split in the 11q− validation data based solely on the expression within the 15-miRNA signature. Pearson's r correlation coefficient was used as a sample object similarity metric. A heatmap was subsequently produced to illustrate this split. These methods were implemented using the k-means and heatmap functions from the stats package of R version 2.10.0.
Patient survival analysis
Kaplan-Meier graphs and log-rank tests for significance were carried out using GraphPad Prism software version 5. Cox proportional hazards regression model was applied to assess the independent effect of each chromosomal aberration on survival. 17q gain was screened for independence of MNA and 11q loss for both EFS and OS. In addition to this model, a screen for independence of 1q loss and 1q gain was also done for both EFS and OS. These methods were implemented using the coxph function in the survival package in R version 2.10.0.
Statistical analysis of miRNA and mRNA signatures
Data from the mirBase miRNA targets database (Microcosm; ref. 22) were used to compile a two-dimensional miRNA/mRNA matrix containing 710 miRNAs, 34,525 transcripts, and 728,289 predicted interactions. All significant (P ≤ 0.05) predicted binding interactions connecting the 15 miRNAs to the various published gene signatures and the average prediction scores were then retrieved. The significance of these interactions was assessed by examining how often these combinations might occur by chance. For each miRNA/mRNA target combination, 100,000 submatrices of equal size were randomly generated and the numbers of connections and scores assessed. The resultant P value was thus used to evaluate the significance of relationship and possible influence the 15-miRNA signature might have on each of the published gene signatures. The above assessment was carried out using in-house generated software implemented in Java version 6 update 15.
Visualization of chromosomal variations and shortest regions of overlap
Global visualization of chromosomal variations for all 160 samples was carried out using an in-house software application implemented using Java and made use of the JFreeChart and Commons packages. Common copy number variants were excluded before visualization. The most commonly detected shortest regions of overlap (SRO) were identified and genomic positions mapped using this application.
DNA copy number analysis was done on 121 primary neuroblastoma tumors using a 72k 4-plex aCGH array. This data set was then combined with a previously reported data set on 39 tumors (23) obtained with a 385k array format for a total analysis of 160 tumors. Figure 1 displays a summary of the aberrations identified across the genome in these tumors, whereas Supplementary Table S1 provides a comprehensive summary of the aCGH results and clinical details for each tumor.
All previously published (4, 6, 8, 9) positive and negative correlations between different genetic abnormalities were detected in our tumor cohort with statistical significance (Table 1). These included inverse correlation between MYCN amplification (MNA) and chromosome 11q loss (P < 0.0001) and positive correlations between loss of 1p and MNA (P = 0.008) and 3p loss with 11q loss (0.0006). Moreover, due to the large size of our tumor set, we were able to identify novel positive correlations between loss of 4p and gain of 7q with loss of 11q (P = 0.002 and P = 0.009, respectively, corrected for multiple comparisons), allowing further definition of the 11q− subtype (Table 1). Loss of 3p and loss 4p were not inversely related within the 11q− subgroup, indicating that these imbalances did not have a complementary effect on tumorigenesis.
Prior studies have indicated that MYCN amplification or 11q loss occurs in approximately 70% to 80% of patients presenting with stage IV disease (15–17). Tumors that lack either of these abnormalities are a much more genetically heterogeneous group that lack predictive biomarkers. The most frequent imbalances found in tumors from our data set lacking either MNA or 11q loss are displayed in Supplementary Table S2 (stage IV disease, n = 20; stages I, II, III, and IV, n = 39). The three most common abnormalities in this subgroup were 17q gain (n = 25; 40%), 16p loss (n = 20; 32%), and 19p loss (n = 19; 31%). Some imbalances, particularly 1q gain, were significantly underrepresented in non-MNA/11q− tumors (n = 1; P = 0.01, corrected for multiple comparisons). There were no discernable positive correlations between genomic imbalances in the non-MNA/11q− tumor group.
Our tumor cohort also exhibited previously published associations between clinical details, genetic abnormalities, and clinical outcome (Supplementary Figs. S1 and S2). As expected, multivariate analysis using the Cox proportional hazards regression model indicated that 1q gain (HR, 2.0; P = 0.04; OS only), 11q loss (HR, 2.1; P = 0.03; EFS only), and MNA (HR, 2.7; P = 0.006; EFS; HR, 2.7; P = 0.008; OS) were independently associated with poor patient survival, although a number of other abnormalities were significantly associated with poor survival in univariate analyses (Supplementary Fig. S2). Unlike the majority of previously published reports (24–26), 17q gain was not an independent predictor of poor survival, as assessed by the Cox proportional hazards regression model. This result, however, is consistent with that of Spitz et al. (27) who reported that 17q gain, in the absence of MNA or 11q−, was not a predictor of EFS or OS.
Shortest regions of common gain or loss of genetic material in tumors
Regions affected by DNA copy number alterations were further analyzed to identify regions commonly deleted or gained in the tumor series. Supplementary Table S3 summarizes these SROs with the approximate genomic position of each breakpoint, the observed frequency of each aberration, and comparison of SROs previously reported. Regions that were classified as copy number variant in the Database of Genomic Variants (http://projects.tcag.ca/variation/) were excluded from this analysis. Detected SROs ranged in size from 0.9 to 10.4 Mb for losses and 0.3 to 6.6 Mb for DNA copy number gains. The most frequent SRO identified was detected on 17q, encompassing a region of 5.2 Mb. This region of common gain contains 18 annotated genes including 5 members of the ATP-binding cassette, subfamily A transporter family (ABCA5, ABCA6, ABCA8, ABCA9, and ABCA10). Interestingly, overexpression of transcripts in this gene family has been previously reported in chronic lymphocytic leukemia (28), high-grade prostatic intraepithelial neoplasia (29), and breast cancer (30) and have also been shown to confer drug resistance on tumor cells (31). Other cancer-related genes located within this 17q SRO include SOX9 (32) and SSTR2 (33).
Two SROs (gains) were identified on chromosome 1q encompassing regions of 1.2 and 3 Mb, respectively. The 3-Mb SRO imbalance was further refined through integration of a previously reported common region of gain on 1q (17). This 1q SRO has been mapped to 206.4-207.7 Mb and contains the genes PLXNA2, mir-205, and LOC642587. PLXNA2 is a member of the plexin family of proteins, which function in association with semaphorin proteins during neural development (34). Although mir-205 has been previously implicated in cancer development (35–37), a screen of 430 miRNAs in 145 neuroblastoma tumors did not identify this miRNA as significant in the pathogenesis of neuroblastoma (38).
Dichotomization of the 11q− genetic subtype by miRNA expression profiling
Fischer et al. (39) recently reported that a 144-gene expression signature originally developed by Oberthuer et al. (40) could stratify 11q− tumors into two subgroups: one with favorable clinical outcome and a second group with unfavorable outcome. To determine if our cohort of 11q− tumors exhibited a similar dichotomous pattern of clinical outcome, we analyzed 11q− tumors (n = 36) using a 15-miRNA expression profiling signature that is predictive of survival, as reported by Bray et al. (38). The signature was trained to identify poor survival on a set of neuroblastomas that included all genetic subtypes of the disease, and then validated on a similar set of tumors that was also broadly representative. The set of 11q− tumors used in the current study is completely independent of the tumors used by Bray et al. (38) to derive the miRNA survival signature. The unsupervised partitional clustering algorithm, k-means, was used to divide the 11q− tumors into two subgroups on the basis of miRNA expression. The two 11q− subgroups differed significantly with respect to EFS (P = 0.04) and OS (P = 0.03; Fig. 2). For group A (n = 16), the probability for 4-year EFS was 14.5% ± 9.3 SE, whereas this increased to 54.7% ± 12.8 SE for group B. The probability of 4-year OS for group A was 25.5% ± 12.5 SE, whereas for group B the rate of survival increased to 65.0% ± 12.8 SE.
Examination of the genomic imbalances found in the favorable (group B) versus unfavorable (group A) 11q− tumor subgroups indicated that no single imbalance was overrepresented with statistical significance. Gain of 17q was present at nearly equal frequencies in the two groups and could not account for the difference in survival for the two different groups of patients with 11q− tumors, consistent with our observation that 17q gain is not independently predictive of survival. Nevertheless, the overall frequency of genomic imbalances was significantly higher in the unfavorable tumor subgroup A (mean 5.5 versus 3.8; P = 0.019; Fig. 2), indicating the possibility that an overall increase in imbalances is associated with poor survival. Interestingly, two of the patients that were identified as having favorable tumors on the basis of miRNA profiling (i.e., false negative results) had a higher number of segmental imbalances (n = 7 in each tumor) than the mean number (5.5) found in the unfavorable tumor subgroup. Thus, an analysis of segmental imbalances following miRNA expression profiling could potentially improve OS prediction. This combined approach led to a sensitivity of 82.4%, a specificity of 55.0%, a positive predictive value of 77.7%, and a negative predictive value of 75.0% on our data set. We conclude that 11q− tumors are composed of two subgroups with very divergent clinical behavior and that these subgroups can be distinguished on the basis of miRNA expression profiling and analysis of DNA copy number alterations.
There was no relationship between patient age at diagnosis and survival for the 11q− tumor data set, as only 6 patients with 11q− tumors were <1.5 years, with 3 of these cases occurring in each group. The relatively small number of patients <1.5 years of age at diagnosis, which is generally associated with better survival, was not surprising given that 11q loss is significantly correlated with older patients (6, 15). In addition, there was also no statistically significant difference in tumor stage between the unfavorable and the favorable subgroups.
Next, we determined whether any of the 15 miRNAs in our survival signature target the genes in the 144-gene signature used by Fischer et al. (39). For the downregulated miRNAs associated with poor survival, 10 potential targets with upregulated genes were identified. However, the number and significance of predicted interactions were not significantly overrepresented (P = 0.20) compared with background (see Materials and Methods). For the upregulated miRNAs associated with poor survival, 13 potential targets with downregulated genes were identified, which approached statistical significance (P = 0.06). More recently, Vermeulen et al. (41) have described a gene signature predictive of clinical outcome for neuroblastoma. A statistically significant number of interactions between downregulated miRs from our signature and upregulated genes from this signature was not evident. However, 18 putative interactions between upregulated miRNAs and downregulated genes from the gene signature were identified, which, according to our simulation studies, is significantly higher than background (P = 0.0007). Some miRNAs from our signature that are upregulated in unfavorable tumors, such as miR-572 (Supplementary Table S4), are predicted to target multiple downregulated genes from the Vermeulen signature (DDC, PRKCZ, PIK3R1, and PTPRF), further indicating that the miRNA signature regulates the gene survival signature to at least some extent. Caution is required in the interpretation of these results, however, as these are not experimentally validated interactions.
Neuroblastoma tumors with loss of chromosome 11q material represent a major genetic subgroup of tumors with adverse clinical outcome. Loss of chromosome 3p was previously shown to occur preferentially within this tumor subgroup (6, 8, 9), and we now further show that loss of chromosome 4p and gain of chromosome 7q are significantly overrepresented in the 11q− tumor subset. Although our results further refine the shortest region of loss or gain on each of these chromosomes, the genes contributing to tumorigenesis on each of these chromosomes remain unidentified. To a lesser extent, 2p gain, 11p gain, and 14q loss are also preferentially associated with the 11q− tumor subtype.
We also show by miRNA expression profiling that 11q− tumors are actually composed of two clinical subtypes with significant differences in both EFS and OS. This finding is completely consistent with a recent report by Fischer et al. (39) showing that 11q− tumors can be split into two subgroups with distinctly different clinical outcomes on the basis of mRNA expression profiling. Moreover, we show by computational methods that miRNAs from our expression survival signature show an enrichment for predicted target sites in the 3′ untranslated region of mRNAs from two published gene survival signatures for neuroblastoma: Obertheur et al. (P = 0.06; ref. 40) and Vermeulen et al. (P < 0.0007; ref. 41). Therefore, miRNAs from the survival signature (38) could potentially influence the gene expression signatures predictive of aggressive disease pathogenesis, at least to a limited extent. It is noteworthy that Mestdagh et al. (42) have also recently shown a significant correspondence between miRNAs that are upregulated either directly or indirectly by MYCN/cMYC and genes that are significantly downregulated in unfavorable neuroblastoma (42). Their miRNA signature was based primarily on MYCN/cMYC regulated miRNAs and therefore differs significantly from the signature used by Bray et al. (38), which was trained solely on the basis of miRNAs associated with poor patient survival. This indicates that differences in miRNA profiles within the 11q− subgroups are not solely MYC driven.
Indeed, the mechanisms regulating the miRNAs in the survival signature are likely to be varied and complex. Bray et al. (38) have shown that expression levels of many miRNA loci are affected by MYCN levels and by large-scale losses and gains of genomic regions in neuroblastoma tumors. The unfavorable 11q− tumors identified by miRNA expression profiling have a significantly higher number of genomic imbalances, relative to the favorable 11q− subgroup, indicating that genomic imbalances might be contributing to the differences in miRNA expression in the two different 11q− subgroups.
Unbalanced gain of 17q, which has been reported to be predictive of poor patient survival by a number of groups (24–26), occurs at nearly equal frequencies in the favorable and unfavorable 11q− subgroups, leading us to conclude that 17q gain does not contribute to the differences in EFS and OS between these groups in a large way. This is consistent with our multivariate analysis of the entire tumor cohort indicating that 17q gain is not independently predictive of poor clinical outcome. Our result is consistent with a previous report by Spitz et al. indicating that 17q gain is not predictive of survival in tumors lacking MNA or 11q loss (27). Thus, our results do not support a role for 17q gain in contributing to an aggressive tumor phenotype, although the overall high frequency of 17q gain in neuroblastoma (including tumors that lack MNA or 11q loss) indicates that it contributes to tumorigenesis in some manner, such as early tumor initiation.
Our finding of an increased frequency of segmental imbalances in the unfavorable 11q− tumor subgroup is completely consistent with reports that the overall pattern of genomic imbalances is highly predictive of clinical outcome in neuroblastoma (43, 44). Janoueix-Lerosey et al. (43) concluded that the presence of segmental alterations is a powerful predictor of poor clinical outcome and can be used to improve patient stratification. Our demonstration that an increased frequency of segmental alterations in 11q− tumors is also correlated with lower EFS and OS provides further confirmation that the two subgroups identified by the miRNA profiling are biologically distinct.
Attiyeh et al. (7) reported following the analysis of a large neuroblastoma tumor cohort that 11q loss by itself is independently predictive of EFS and progression-free survival but is not an independent predictor of OS. Loss of 11q was also significantly associated with EFS, but not OS, in our data set. The specific application of a 15-miRNA expression signature, combined with an analysis of segmental imbalances, however, allows for the identification of poor OS for patients with 11q− tumors with a sensitivity of 82%. The signature seems to be better at identifying truly aggressive 11q− tumors than genuinely favorable tumors, as specificity was only 55%.
In conclusion, we show that two 11q− neuroblastoma tumor subgroups occur, which differ significantly in EFS and OS. These subgroups differ significantly in their miRNA expression profiles and, as shown by Fischer et al. (39), in their mRNA expression profiles. The tumor subtype with poor EFS and OS has significantly more segmental genomic imbalances than the favorable subgroup, indicating that a combination of expression profiling (miRNAs and mRNAs) with analysis of DNA copy number alterations will lead to improved prognostication of this often fatal tumor subtype.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Grant Support: Grants from Science Foundation Ireland (07/IN.1/B1776), Children's Medical and Research Foundation, the NIH (5R01CA127496), the Belgian Foundation Against Cancer, Foundation of public interest, the Fund for Scientific Research (G.0198.08), GOA (01G01910), and the VLK (Vlaamse Liga Tegen Kanker). P.G. Buckley is supported by a postdoctoral fellowship from the Irish Research Council for Science, Engineering and Technology. K. De Preter is a postdoctoral researcher with the Fund for Scientific Research-Flanders. P. Mestdagh was sponsored by BOF (01D31406).
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.
- Received December 8, 2009.
- Revision received March 31, 2010.
- Accepted April 15, 2010.