System-Level Analysis of Neuroblastoma Tumor–Initiating Cells Implicates AURKB as a Novel Drug Target for Neuroblastoma

Neuroblastoma (NB) is the most common and fatal extracranial solid tumor of children. Despite aggressive chemotherapy and radiation therapy, <40% of patients with high-risk disease achieve long-term survival (1). The failure of current therapies is due to a high rate of tumor relapse that is lethal in most cases (2). New therapies that can prevent or treat disease recurrences are therefore necessary to achieve better prognosis for children with high-risk NB.

It has been suggested that failure of chemotherapy, radiation therapy, and targeted approaches to permanently cure aggressive cancers such as high-risk NB is due to inherent drug resistance of cancer stem cells or tumor-initiating cells (TICs) that survive the treatment and subsequently regenerate the tumor (3). Cancer stem cells or TICs have been described in a variety of hematopoietic and solid malignancies, including those of the breast, brain, pancreas, liver, skin, and colon (4). Primary TIC lines have also been isolated from NB tumors and metastases from patients, and these cells recapitulated metastatic NB in an orthotopic mouse model with as few as 10 cells (5). NB TICs and cancer stem cells share several properties, including the ability to self-renew and differentiate into cell types observed in the bulk tumor, express stem cell markers, and exhibit enhanced tumorigenic potential as compared with established cell lines (5). A recent study using chronic myeloid leukemia stem cells provided proof of principle that targeting a cancer stem cell–enriched gene could lead to the eradication of such cells and a potential disease cure (6). Therefore, NB TICs, which are primary cell lines with very high tumorigenic potential in immunosuppressed mice provide a model for the development of improved therapies for recurrent and metastatic NB.

We describe here an RNA sequencing approach (7, 8) we used to characterize a panel of NB TICs and normal neural crest-like skin-derived precursor cells (SKP) isolated from human foreskins (9), aiming to identify transcripts preferentially abundant in NB TICs. SKPs are multipotent precursors isolated from human foreskin that are able to self-renew and differentiate into various neural crest derivatives, including peripheral neurons and Schwann cells (9). Because NB is a tumor of neural crest precursors (10), SKPs provide a normal reference transcriptome for the identification of candidate gene expression changes associated with the TIC phenotype (5). To increase the specificity of the identified gene expression changes, we also used RNA sequencing data derived from breast, skin, brain, lymph node, ovary, cervix, and lung tumors, and cell lines to identify transcripts enriched in abundance in NB TICs compared with other tumor types.

To identify existing therapeutics that could be applied to the treatment of recurrent NB, we used Ingenuity Pathway Analysis software (Ingenuity Systems) to analyze the functional significance of identified genes and match them against a database of available drugs. Our work implicates an existing Aurora B kinase inhibitor, currently in clinical trials for acute myeloid leukemia, as a drug repositioning candidate for recurrent NB. This study provides the first analysis of NB TIC transcriptomes and illustrates the principle that TIC-specific expression profiles are useful for pinpointing candidate novel therapeutic options for NB.

NB TICs and SKPs were cultured as previously described (5,9). Briefly, the cells were cultured in DMEM-F12 medium, 3:1 (Invitrogen), containing 2% B27 supplement (Gibco), 40 ng/mL basic fibroblast growth factor 2, and 20 ng/mL epidermal growth factor (both from Collaborative Research; proliferation media) in 75 cm² flasks in a 37°C and 5% CO₂ tissue-culture incubator. The cell growth conditions were normalized such that NB TICs were cultured for 7 days and SKPs for 14 days postplating prior to harvesting in exponential growth phase and RNA isolation for transcriptome analysis. Details of the NB TIC and SKP samples used in this analysis are outlined in Table 1.

RNA sequencing libraries from NB TICs and SKPs were constructed from DNase I treated mRNA as previously described (11). The libraries were sequenced on an Illumina Genome Analyzer. The read length and number of reads generated for each library is provided in Supplementary Table S1. The reads were aligned to the human reference genome (National Center for Biotechnology Information Build 36.1) and a database of known exon junctions (7) using MAQ software version 0.7.1 in paired-end mode, and the duplicate read pairs were removed (12). The number of bases sequenced per the number of exonic bases mapped was used as a measure of gene expression level for each gene (11). The expression values were square-root transformed and used in the lmFit function to estimate fold changes between the compared groups by fitting a linear model (13). NB TICs versus SKPs and NB TICs versus other cancers (Supplementary Table S1) were compared this way. The F-statistical with Benjamini-Hochberg multiple testing correction implemented in the eBayes function was used to assess the significance of differential expression. Those genes with Benjamini-Hochberg–corrected P < 0.05 were considered statistically significant.

Cells were collected and lysed in Trizol, and RNA was purified using RNeasy mini kit (Qiagen). RNA samples (Table 1) were analyzed on Affymetrix GeneChip Human Exon 1.0 ST Arrays. The data were background corrected and normalized using the Robust Multichip Average procedure implemented in the Affymetrix Expression Console software. Gene-level expression summaries were computed based on all core probes. Differential gene expression was assessed using the lmFit function of the Linear Models for Microarray Data (LIMMA) Bioconductor package (13). The significance of differential expression was assessed as described above for the sequencing data.

A list of significantly differentially expressed genes from each analysis (NB TICs versus SKPs RNA sequencing, NB TIC versus tissue pool RNA sequencing, and NB TIC versus SKPs microarray) was compared against the lists from the other two analyses to derive sets of common upregulated and downregulated genes. Ingenuity Pathway Analysis software (Ingenuity Systems) was then used on these sets to select canonical pathways significantly enriched among each set (P < 0.05).

A crude membrane fraction was prepared as follows. NB88R2 cells were swollen in hypotonic buffer (20 mmol/L Tris, pH 7.4; 10 mmol/L KCl; 5 mmol/L sodium vanadate; 1 mmol/L phenylmethylsulfonylfluoride) and lysed by dounce homogenization. The cleared cell lysate was centrifuged for 15 minutes at 6,000 × g to collect the crude membrane fraction. The protein fraction was resuspended in urea buffer (8 mol/L urea, 2 mmol/L HEPES, 2.5 mmol/L sodium pyrophosphate, 1 mmol/L β-glycerophosphate, and 1 mmol/L vanadate; Cell Signaling Technology) and was reduced and alkylated with 4.5 mmol/L dithiothreitol (DTT) and 10 mmol/L iodoacetamide, respectively.

Whole-cellular fraction was prepared as follows. NB88R2 cells were lysed in urea lysis buffer (8 mol/L urea, 2 mmol/L HEPES, 2.5 mmol/L sodium pyrophosphate, 1 mmol/L β-glycerophosphate, and 1 mmol/L vanadate) and sonicated (3 bursts of 4 W for 10 s). The cell lysate was cleared by centrifugation (20,000 × g for 15 min at 4°C) and was reduced and alkylated with 4.5 mmol/L DTT and 10 mmol/L iodoacetamide, respectively. Proteins were digested with trypsin and purified using C18 reverse phase resin prior to mass spectrometry.

The gel-free two-dimensional liquid chromatography coupled to shotgun tandem mass spectrometry (MudPIT) analysis was done as described (14) with the following modifications: approximately 60 μg (membrane fraction) or 40 μg (whole-cell fraction) of digested protein was analyzed on a linear ion-trap LTQ-Orbitrap mass spectrometer (ThermoFisher). Samples were loaded using a Proxeon HPLC system (Thermo Fisher Scientific) and subjected to an eight-cycle MudPIT. All data was analyzed using Sequest (ThermoFinnigan; version SRF v. 5) and X! Tandem (http://www.thegpm.org/; version 2007.01.01.2 for membrane fraction or version TORNADO 2009.04.01.3 for whole-cell fraction) search algorithms using the Human International Protein Index database (version 3.41 with 72,155 entries or version 3.66 with 86,845 entries for membrane and whole-cell fractions, respectively). Sequest and X! Tandem were searched with a fragment ion mass tolerance of 0.50 or 0.40 Da for membrane and whole-cell fraction, respectively, and a parent ion tolerance of 2.0 or 5.0 ppm for membrane or whole-cell fraction, respectively. The iodoacetamide derivative of cysteine was specified as a fixed modification in Sequest and X! Tandem. The oxidation of methionine was specified as a variable modification. Proteins were accepted based on the following criteria. At least two peptides per protein were identified with a probability threshold of 95% or greater as derived by the Peptide Prophet algorithm (15) and an overall protein identity of >95.0% using the Protein Prophet algorithm (16).

NB12, NB88R2, and FS283 were dissociated into single cells and seeded in triplicates at 3,000 cells per well in 50 μL medium containing 30% SKPs conditioned media in non–tissue culture–treated 96-well plates (Corning Life Sciences). AZD1152 (Selleck Chemicals LLC) was dissolved in dimethyl sulfoxide (DMSO) to a stock concentration of 50 mmol/L, from which 1:3 fold sequential dilutions were prepared. Intermediate dilutions of the compound were made in medium and immediately added to the cells in a volume of 50 μL. Cells treated with 0.05% DMSO in the absence of the drug were used as a control for optimal cellular proliferation, whereas wells containing media only were used to determine the background fluorescence; alamarBlue (10 μL) was added to each well after 72 hours, followed by incubation for an additional 24 hours. Fluorescence intensity was measured using PHERAstar SpectraMax Plus384 microplate reader (BMG Labtech) with an excitation filter of 535 nm and an emission filter of 590 nm. Percentage reduction of alamarBlue was calculated as ((mean fluorescence of treated wells - background fluorescence)/(mean fluorescence of DMSO-treated wells - background fluorescence)) * 100. Half maximal effective concentration (EC50) curves were generated using GraphPad Prism 5 software (GraphPad Software, Inc.).

Cells were harvested, washed with cold HBSS, and lysed with NP40 lysis buffer containing 10 mmol/L Tris (pH 8.0), 150 mmol/L NaCl, 10% glycerol, 1% Nonidet P-40, 1 mmol/L phenylmethylsulfonylfluoride, 1 mmol/L orthovanadate, and proteinase inhibitor cocktail tablet (Complete, Mini, EDTA-free, Roche). Cells were lysed for 10 to 20 minutes on ice and centrifuged for 10 minutes at 12,000 rpm at 4°C. Protein amounts were determined by BCA Assay (Pierce), and 40 μg of protein was loaded per lane. Western blots were probed with rabbit polyclonal anti-Aurora B antibody (Abcam; ab2254) and anti–glyceraldehydes-3-phosphate dehydrogenase (Santa Cruz; sc-47724) antibody in 5% w/v nonfat dry milk in TBS/0.1% Tween-20 over night at 4°C. Blots were developed using ECL or ECL-plus reagent (GE Healthcare Life Sciences).

Cell lines were infected with either a mock treatment or lentivirus-encoding shRNAs of interest at a multiplicity of infection of 1.0. Seventy-two hours postinfection, the virus was removed, and cells were seeded in triplicate at a density of 10,000 per well in 24-well plates. The remaining cells were used for RNA isolation to determine the efficiency of knockdown by quantitative reverse transcriptase (qRT) PCR. Viable cell numbers were determined on days 1, 3, 5, and 7 postplating by removing cells from wells and counting via hemocytometer.

TICs have been identified and characterized in NB primary tumors and metastases, and have been shown to be associated with tumor relapse (5). Genes whose expression is altered in NB TICs compared with normal equivalent cells would thus provide potential targets for refractory and recurrent NB. To identify such genes, we sequenced transcriptomes from 10 NB TIC lines isolated from tumors and metastases of six high-risk NB patients (Table 1) using Illumina RNA sequencing (7, 8). The NB TIC lines used in this study included those isolated from patients during disease relapse and one from a patient in remission. Because we previously showed that line NB67 isolated from the bone marrow of a patient in clinical remission who subsequently relapsed was tumorigenic, we included this NB TIC line in the analysis (5). To generate reference normal expression profiles, we sought a primary human nontransformed or immortalized cell that would as closely as possible resemble the putative origin of NB cells, the human embryonic neural crest. SKPs are primary (nonimmortalized) human neural crest-derived precursor cells that are resident in children and give increase in culture and in vivo to the same differentiated cell types as embryonic neural crest stem cells, including peripheral neurons, bone and cartilage, smooth muscle, and Schwann cells (9). SKPs have a normal karyotype and are cultured under the same conditions as NB TICs; their use as reference normal cells for the analysis of NB TICs has been reported (5). Therefore, for this comparative transcriptome study, we also sequenced the transcriptomes of three SKP lines from three children without cancer diagnosis (9).

We used the LIMMA package to identify genes with evidence of differential expression between NB TICs and SKPs. This analysis revealed 817 and 1,913 genes either increased or decreased in abundance in NB TICs versus SKPs.

We considered it likely that, within the list of differentially expressed genes, there would be candidate NB TIC markers and also transcripts generally associated with a proliferative phenotype. Targeting gene products that are nonspecifically expressed in proliferating cell types would potentially result in increased toxicity, particularly in children whose organ systems are undergoing growth and development. To identify gene expression differences specific to NB TICs, we compared our NB TIC RNA sequences to RNA sequencing data from 30 cancer samples available in house. These samples were derived from seven tissue types, including ovary, lymph nodes, lung, blood, brain, skin, and cervix (Supplementary Table S1) and were included as an additional reference set for the identification of transcripts enriched specifically in NB TICs. Our analysis revealed 449 and 1,059 genes were upregulated or downregulated in the comparison of NB TICs versus SKPs and NB TICs versus other tissues, respectively.

To confirm the differential expression of candidate NB TIC-enriched transcripts identified using RNA sequencing, we analyzed eight NB TIC lines from five patients and five SKP lines from four cancer-free children (Table 1) using Affymetrix Human Exon 1.0 ST Arrays. This platform provides independent confirmation of gene expression at the level of exons (17, 18). Exon array analyses confirmed the differential expression of 321 (71%) and 819 (77%) genes in NB TICs compared to other tissues, respectively (Fig. 1; Supplementary Table S2). These genes represented robust sets of NB TIC-enriched and depleted transcripts that we analyzed further to identify the pathways disrupted in NB TICs.

A robust prognostic gene expression signature of 59 genes was recently developed using primary NB microarray data sets and literature review; high expression levels of 17 of the genes were associated with high-risk disease, whereas high expression levels of 42 other genes were associated with low-risk disease (19). We determined the overlap of these gene lists with the lists of NB TIC-enriched and -depleted transcripts we obtained as described above. Of the 321 NB TIC-enriched genes, six (ODC1, MCM2, TYMS, PAICS, CDCA5, BIRC5) were part of the 17-gene high-risk signature. On the other hand, none of the NB TIC-enriched genes were part of the 42-gene low-risk signature, supporting an association between NB TICs and high-risk rather than low-risk disease. Moreover, three of the genes linked to low-risk disease (AKR1C1, PMP22, and PLAT) were detected among our list of 819 NB TIC-depleted genes. Our analysis suggested that NB TIC expression profiles are at least partially consistent with known molecular features of primary NB tumors and that NB TICs therefore have utility as models of NB. On the other hand, because 9 of 10 NB TIC lines we examined were derived from metastases rather than primary tumors, NB TIC gene expression signature may differ substantially from that of primary NB tumors.

To assess the functional significance of transcripts differentially abundant in NB TICs, we conducted a pathway enrichment analysis using the Ingenuity software (Ingenuity Systems). The analysis revealed several signaling pathways significantly upregulated in NB TICs (P < 0.05), including the BRCA1 DNA damage response pathway, the mitotic role of pololike kinase pathway, the role of CHK proteins in cell cycle checkpoint control pathway, the cell cycle G₂/M DNA damage checkpoint regulation pathway, and the ATM signaling pathway (Fig. 2A). In contrast, the axonal guidance pathway, the CXCR4 signaling pathway, the integrin-linked kinase (ILK) signaling pathway, the germ cell-sertoli cell junction signaling pathway and the transforming growth factor β signaling pathway were significantly downregulated in NB TICs compared with SKPs and the tissue pool (Supplementary Fig. S2). Of the 321 genes significantly upregulated in NB TICs, 13 were known members of the BRCA1 DNA damage response pathway (Fig. 2B). In addition, eight genes were associated with pololike kinase and cell cycle checkpoint control pathways that are direct downstream targets of BRCA1 signaling (Fig. 2B). The elevated ATM expression level in the NB TICs could be associated with the stem cell characteristics of these cells because a recent study has found that ATM expression is required for the proliferation of neural stem cells and DNA repair (20).

Supporting the notion that the BRCA1 pathway is a key signaling cascade in NB, the BRCA1 pathway involves the protein product of the BARD1 locus, common variations in which have been linked to susceptibility to high-risk NB by a single nucleotide polymorphism–based genomewide association study (21). In addition, a recent study reported that BRCA1 in complex with BARD1 possesses ubiquitin ligase activity that can modify histones H2A and H2B, thereby affecting nucleosome dynamics (22). Our result is therefore consistent with the notion that misregulation of the BRCA1 pathway may contribute to genetic and epigenetic aberrations observed in high-risk NB (23).

To assess the contribution of the NB TIC-enriched genes to the NB TIC proteome, we conducted a MudPIT analysis (14) of whole-cell lysate and a membrane-enriched fraction of NB TIC line NB88R2 generated from a bone marrow metastasis of a high-risk patient. The MudPIT approach effectively identifies hundreds to thousands of concurrently expressed proteins for global or subcellular fraction–specific proteomic profile analyses (24). The MudPIT analysis of the whole-cell lysate isolated from line NB88R2 cells revealed 819 proteins in which each protein was identified by at least two peptides. A similar analysis identified 1,530 proteins in the membrane-enriched fraction isolated from the same line.

Of the 321 TIC-enriched genes, we identified in our transcriptome analysis, 75 were detected by MudPIT in either whole-cell or membrane-enriched lysate of line NB88R2 or both (Supplementary Table S3). Forty-five of the detected proteins were products of genes that were expressed in the 75% to 100% expression percentile in the NB88R2 line, whereas only two protein products were detected for genes expressed in the 0% to 24% expression percentile, showing a correlation between transcript abundance and MudPIT analysis (Supplementary Fig. S2). Twenty-one percent (16 of 75) of the detected proteins were associated with DNA replication, recombination, and repair functional category (Ingenuity Systems), including PARP1, PCNA, UBE2N, FEN1, HMGB2, and RFC, which forms a major complex interacting with BRCA1 (25). This result suggests a general correlation between the gene expression and protein level in at least one NB TIC line.

Because the most direct pharmacologic intervention is inhibition of a target protein (26), we focused further functional analyses on genes upregulated rather than downregulated in NB TICs with respect to SKPs and other tissues. Drug repositioning, in which existing drugs are used for novel indications, is a powerful approach to novel therapy development because it greatly reduces the cost and time required to clinically develop a new therapeutic option (27). We therefore aimed to use NB TIC-enriched genes to identify targets of existing therapeutics with the concept that such drugs could in the future be potentially effective against recurrent NB. We applied the Ingenuity Pathways Analysis tool to map the 321 NB TIC-enriched transcripts, as well as their interacting partners, to known drugs. This analysis revealed thirty known drug targets among the NB TIC-enriched genes and their interacting partners defined by Ingenuity Knowledge Base (bold type in Table 2 indicates the NB TIC-enriched genes). Many of the predicted drug targets have been explored preclinically or clinically for the treatment of NB (Table 2). Drugs predicted by our analysis included both general chemotherapeutics, such as etoposide, becatecarin, doxorubicin, flavopiridol, and vincristine, all of which, we coincidentally learned, are currently approved or in trials for NB, as well as targeted agents such as BCL2 inhibitors, evaluated for the treatment of NB (28). Several agents predicted by our analysis, such as HDAC inhibitors and PARP inhibitors, have already shown promise in the management of chemotherapy-resistant NB (29, 30). In addition, our analysis predicted genes and gene products targeted by existing drugs that previously have not been implicated clinically as therapeutic targets for high-risk NB. These molecules include AURKB, ADORA2A, CXCL10, SLC1A4, COL14A1, TNFRSF10B, ITGA2b, and IL6.

The Aurora kinase family includes three serine/threonine kinases involved in the control of the cell cycle. Aurora A and B kinases have shown promise as anticancer agents for the treatment of solid tumors and leukemias (31). Although an Aurora A kinase inhibitor is in an ongoing phase I/II clinical trial for NB (NCT00739427), Aurora B kinase inhibitors have not been investigated in relation to NB. A recent report suggested a direct link between Aurora B kinase and BARD1, a key component of the BRCA1 signaling pathway that is also associated with susceptibility to NB (32). That report, together with the aberrant expression of the BRCA1/BARD1 pathway observed in our study, provided a rationale for exploring the antiproliferative potential of Aurora B kinase inhibitors in NB TICs.

To assess whether elevated mRNA levels at the AURKB locus in NB TICs corresponded to increased levels of AURKB protein, we did Western blot analysis using whole-cell lysates from three NB TIC lines (NB12, NB88R2, and NB122R) and two SKP lines (FS274 and FS227). This analysis revealed a strong presence of the AURKB protein in NB TICs but detected no protein in SKPs, supporting the gene expression result (Fig. 3A). To gain further insight into the role of AURKB in controling NB TIC proliferation, we performed shRNA knockdown experiments in NB TIC lines NB12 and NB88R2. NB TICs infected with lentiviruses encoding two separate shRNAs to AURKB showed 77% to 80% growth reduction compared with NB TICs infected with lentiviruses carrying mock shRNAs to green fluorescent protein or β-galactosidase (Fig. 3B and C). The observed reduction in proliferation following AURKB knockdowns supports the premise that AURKB signaling is important for the viability of NB TICs.

To assess whether pharmacologic inhibition of AURKB would have the same effect on NB TIC proliferation as the AURKB knockdowns done above, we used AZD1152, a selective AURKB inhibitor that is currently undergoing phase I/II testing in patients with acute myelogenous leukemia (NCT00497991). NB TIC lines (NB12 and NB88R2), as well as FS283 SKP line, were treated with a range of AZD1152 concentrations, and cell growth was assessed 96 hours later using alamarBlue reduction (33) as a read-out of cellular metabolic activity. As shown in Fig. 3D, proliferation of NB TICs is reduced following inhibition of AURKB, showing low micromolar EC50 values (1.5-4.6 μmol/L). In contrast to this, SKPs were less sensitive to AZD1152, exhibiting higher EC50 values (12.4 μmol/L). The selective activity of AZD1152 in NB TICs, which is likely due to the differential protein abundance in NB TICs compared with SKPs, provides a foundation for further exploring AURKB as a drug target for pediatric NB. Our work provides the first high-resolution system-level analysis of NB TICs and a proof of principle that inhibiting an NB TIC-enriched pathway is selectively cytotoxic to these cells while sparing their normal stem cell counterparts. The apparent selectivity of AURKB inhibition is particularly important in pediatric oncology as a pool of normal stem cells must be maintained for proper development.

No potential conflicts of interest were disclosed.

We thank Robyn Roscoe and Diane Miller for project management support; Mark Robinson for statistical advice on the use of the LIMMA package for the analysis of RNA sequencing data; and the Genome Sciences Centre Library Core, Sequencing, and Bioinformatics teams for their technical assistance with this work.

Grant Support: Canadian Institutes of Health Research, the National Cancer Institute of Canada, The James Birrell Fund for Neuroblastoma Research, Lilah's Fund for Neuroblastoma Research, McLaughlin Centre for Molecular Medicine, Solving Kid's Cancer, the Stem Cell Network, The British Columbia Childhood Cancer Parents' Association, The Jordan Hopkins Foundation for Cancer Research, Genome Canada, Genome British Columbia, The British Columbia Cancer Foundation, The Terry Fox Research Institute, The Ontario Institute for Cancer Research, The Michael Smith Foundation for Health Research training (O. Morozova and R.D. Morin) and scholarship (M.A. Marra), Natural Sciences and Engineering Research Council of Canada fellowship (O. Morozova), Vanier Scholar of the Canadian Institutes of Health Research (R.D. Morin), and Canada Research Chairs (D.R. Kaplan F. Miller, and M.F. Moran).

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.