Genome-Wide Microarray Expression and Genomic Alterations by Array-CGH Analysis in Neuroblastoma Stem-Like Cells

Neuroblastoma has a very diverse clinical behaviour: from spontaneous regression to a very aggressive malignant progression and resistance to chemotherapy. This heterogeneous clinical behaviour might be due to the existence of Cancer Stem Cells (CSC), a subpopulation within the tumor with stem-like cell properties: a significant proliferation capacity, a unique self-renewal capacity, and therefore, a higher ability to form new tumors. We enriched the CSC-like cell population content of two commercial neuroblastoma cell lines by the use of conditioned cell culture media for neurospheres, and compared genomic gains and losses and genome expression by array-CGH and microarray analysis, respectively (in CSC-like versus standard tumor cells culture). Despite the array-CGH did not show significant differences between standard and CSC-like in both analyzed cell lines, the microarray expression analysis highlighted some of the most relevant biological processes and molecular functions that might be responsible for the CSC-like phenotype. Some signalling pathways detected seem to be involved in self-renewal of normal tissues (Wnt, Notch, Hh and TGF-β) and contribute to CSC phenotype. We focused on the aberrant activation of TGF-β and Hh signalling pathways, confirming the inhibition of repressors of TGF-β pathway, as SMAD6 and SMAD7 by RT-qPCR. The analysis of the Sonic Hedgehog pathway showed overexpression of PTCH1, GLI1 and SMO. We found overexpression of CD133 and CD15 in SIMA neurospheres, confirming that this cell line was particularly enriched in stem-like cells. This work shows a cross-talk among different pathways in neuroblastoma and its importance in CSC-like cells.


Introduction
Neuroblastoma is a neuroendocrine tumor of unknown etiology, derived from primordial neural crest cells which afterwards develop into adrenal medulla and sympathetic ganglia [1]. It is one the most prevalent cancers in childhood and nearly 50% of the cases take place in children younger than two years. The estimated 5-years survival rate is 90% and 50% for patients with non-high-risk and high-risk neuroblastoma, respectively. Most neuroblastomas occur sporadically and develop with very diverse clinical behavior, from spontaneous regression to aggressive malignant progression and resistance to chemotherapy [2,3].
The actual treatment of neuroblastoma depends on the clinic stage of the tumor, but it is commonly based on radiopharmaceutical therapy in combination with surgery [4]. Current studies and clinical trials are combining conventional chemotherapy with monoclonal antibodies, stem cell transplants and retinoids, but because of the complexity of this pathology, progress remains extremely slow [5].
Some authors propose that this diverse clinical behavior of neuroblastoma might be due to molecular differences in cell subpopulations [6]. The cancer stem cells (CSC) model might be an explanation for this heterogeneous behavior [7]. This model proposes that only a small subpopulation with characteristics of stem cells within the tumor has the ability to proliferate and maintain its growth. Even if a tumoral mass shows a substantial decrease in size in response to therapy, if the CSC are spared, it will regrow leading to a relapse [8,9]. Some studies propose that CSC operate with the machinery and developmental programs expressed in normal stem cells [10,11]. There is growing evidence of some signalling pathways involved in self-renewal of both normal and tumor tissues as Wnt [12,13], Sonic Hedgehog (Hh) [14,15], Notch [16,17] and Transforming Growth Factor Beta (TGF-b) [18][19][20] signalling pathways, that might contribute to tumorigenesis when deregulated [14,21].
In this line, the development of new therapies based in molecular targets may be of great value for the treatment of neuroblastoma. The observation of some pathways acting on multiple levels to promote the development of neuroblastoma and CSC subpopulation, has prompted new therapeutic strategies to treat not only this neoplasm but other brain and nervous system tumors [2,18,22]. Therefore, this project performs a genomic analysis of CSC-like by array-CGH and expression array, with special focus on altered signalling pathways that might explain the stem cell phenotype of the CSC subpopulation. The aim of this study is to provide a powerful tool to open up new targets for therapy or redirect current cancer treatments towards CSC in order to achieve total elimination of tumor cell population and improve treatment effectiveness.

Cell lines culture
Two commercial neuroblastoma-derived cell lines were used: SK-N-DZ cell line (ATCC NuCRL-2149) provided by the American Type Culture Collection (ATCC, Manassas, VA, USA) and SIMA cell line (DSMZ Nu ACC 164) provided by the German Collection of Microorganisms and Cell Cultures (DSMZ, Braunschweig, Germany). Both were grown as monolayer with DMEM/GlutaMAX medium supplemented with 10% fetal bovine serum, 5% non-essential aminoacids, 1% penicillin/ streptomycin and 0.1% amphotericin B. To obtain the neurospheres cultures in order to enrich them in CSC-like cells, after chemical dissociation of SK-N-DZ and SIMA cell lines, 5610 5 cells were transferred to 25 cm 2 flasks (positioned vertically) and grown in suspension with 5 ml of selective medium: DMEM-F12/ GlutaMAX plus 1% penicillin/streptomycin and 0.1% amphotericin B, supplemented with Epidermal Growth Factor (EGF) (20 ng/ml), Fibroblast Growth Factor Basic (FGFb) (20 ng/ml) and B27 Supplement (1X). All cultures were maintained at 37uC in a humidified atmosphere of 5% CO 2 /95% air.

DNA/RNA isolation and reverse transcription
DNA from 1610 6 cells pellet was extracted using Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA) and total RNA from 2610 6 cells pellet, employing Illustra RNAspin Mini Kit (GE Healthcare, Buckinghamshire, UK) following manufacturer's instructions in both cases. DNA and RNA samples were stored at 220uC and 280uC until utilization, respectively. DNA and RNA concentrations were measured with a SmartSpec Spectophotometer (Bio-Rad, Hercules, CA, USA). Reverse transcription was carried out from 2 mg of total RNA (previously denatured at 65uC for 5 min) with a mixture of 0.5 mg of random primers, 0.5 mM dNTPs, 10 mM DTT, First-Strand buffer 1X and 200U SuperScript II Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA) in a total volume of 40 ml. The mixture was incubated for 10 min at 25uC, 50 min at 42uC and 15 min at 70uC, followed by chilling on ice. The cDNA was diluted 1/5 and stored at 280uC until utilization.

Array-CGH
DNA samples were analyzed using a 244K microarray (Agilent Technologies, Santa Clara, CA, USA). Oligonucleotide aCGH processing was performed as detailed in the manufacturer's protocol (version 4.0; http://www.agilent.com). Data were extracted from scanned images using feature extraction software (version A. 8.5.3,Agilent). Raw data text files from the latter were then imported for analysis into CGH Analytics 3.4.40. Aberrations were detected with the ADM2 algorithm and filtering options of a minimum of 5 probes and abs(log2Ratio).0.3. Aberration segments were individually reviewed using build 35, hg17 of UCSC. Anomalies that were localized to regions with high-copy repetitive or GC-rich DNA sequences including telomeric regions were excluded. Gains and losses for the oligonucleotide dataset were defined as a linear ratio $1.2 or #0.8, respectively. High and low-level amplification events were defined as a linear ratio $4 or 2,ratio,4, respectively. The data are described in accordance with MIAME guidelines and have been deposited in ArrayExpress under E-MTAB-2866 accession number. All the analysis and statistics related with the array-CGH was performed at Institut Gustave Roussy by the bioinformatics team.

Microarray experiment, data normalization and analysis
SurePrint G3 human Gene Expression Microarray Kit (Agilent, ID:028004) was used to analyze the transcriptional profiles of all samples. This array contains 42,404 60-mer oligonucleotides targeting 27,988 Entrez Gene RNAs and 7,419 lincRNAs (long intergenic non coding RNAs) in an 8660K format slide (each slide contains 8 arrays of 60,000 features with 1440 spots reserved for internal quality control). The transcripts used can be checked in the EMBL-EBI Database (http://www.ebi.ac.uk/arrayexpress/ arrays/A-GEOD-14550).
Array performance and analysis were carried out at the Institut Gustave Roussy, using the Agilent Feature Extraction software. Data sets were normalized using quantile normalization, and fold change between cell line and neurospheres was calculated for each gene. The whole dataset was filtered by intensity of the probes. Spots with too high intensities (saturating intensities) and irregular spots (spots inside whose neighbouring pixel showed very different intensities) were not further considered. Likewise, genes with a signal-to-noise ratio lower than 6 fold the average standard deviation observed for all negative control probes were excluded. Finally, .1 fold-change between normalized gene expression measured in neurospheres (or CSC-like cells) versus standard cell lines was used as the criteria to select differentially expressed genes.
The list of differentially expressed genes in both cell lines was imported and classified by DAVID (Database for Annotation, Visualization and Integrated Discovery, http://david.abcc.ncifcrf. gov/), version 6.7 [23][24][25]. Enriched categories in biological processes and molecular function were defined by a p value ,0.05 and more than 2 genes differentially expressed in each category. After the enrichment analysis, a classification by PANTHER (Protein ANalysis THrough Evolutionary Relationships, www. pantherdb.org/) Data Base [26][27][28] was performed. Finally, a deeper pathway analysis was carried out by comparing the same set of genes with the NCBI BioSystems Database (http://www. ncbi.nlm.nih.gov/biosystems). Microarray expression data are available at ArrayExpress with accession number E-MTAB-2867.

Real-Time polymerase chain reaction (RT-qPCR)
A first screening of our primer library was performed according to the microarray data and pathway analysis. 44 pairs of primers were selected and aligned with the mRNA and DNA sequences obtained from RefSeq (http://www.ncbi.nlm.nih.gov/RefSeq/) and UCSC Genome Browser (http://genome.ucsc.edu/). A deeper analysis of each pair of primers was performed with Primer3 (http://primer3.sourceforge.net/) and Primer Blast (http://www.ncbi.nlm.nih.gov/tools/primer-blast/) in order to obtain product length, %GC, melting temperature (Tm), selfcomplementarity or possible products on unintended templates. After the computational analysis, 17 out of 44 primers (Table 1) were selected for further studies. RT-qPCRs were performed in order to define the optimal Tm for each pair of primers. The melting curves obtained were analyzed to check amplification efficiency in each temperature, and detect additional peaks displaced from the desired amplicon peak (determined by agarose gels checking size and identifying possible intron amplifications). For the quantitative PCR reaction, 2.5 ml of the generated cDNA was added to a PCR mix containing 12.5 ml SYBR Green Supermix (Bio-Rad, Hercules, CA, USA) and 12.5 pM forward and reverse primers, in a total volume of 25 ml. The reaction was performed according to the following protocol: 10 min at 95uC, followed by 40 cycles of 15 s at 95uC, 30 s at the specific Tm for each primer (Table 1) and 20 s at 72uC. A melting curve was added at the end of the protocol as a quality control. All samples were run in triplicate and three independent experiments were carried out. Gene expression levels between samples were normalized using the expression levels of the HPRT1 gene. Relative gene expression was analyzed according to the 2 2DDCt method [29].

Statistical analysis
The statistical analysis of arrays are detailed in the description for each technique. The data of gene expression measured by RT-qPCR was graph as the mean 6 standard deviation. The 2 2DDCt values were analysed for normality distribution and homogeneity of variance. Since the data did not show normality distribution, the statistical significance between groups was compared by pairs (the expression level of each gene in standard cultures vs neurospheres cultures) with the non-parametric Mann-Whitney U test. The significance level was set as p,0.05. The statistical program GraphPad Prism for Windows version 5.04 was used.

Genomic profile by array-CGH analysis
To identify the potential common chromosomal alterations indicating possible zones implicated in the CSC-like phenotype, an array-CGH analysis was performed. In the case of SK-N-DZ, we only found six chromosomal regions amplified and no deleted region. SIMA cell line showed more regions of gain and loss, specifically 10 and 30 regions, respectively (data summarized in Table 2). Notwithstanding we did not find any region that was lost or gained in both cell lines and the significance was found only in some points of the genome, we identified numerous chromosomal sequences with copy number variants (CNV) both in SK-N-DZ and SIMA cell lines (Table 2). Whole genome profiles are presented in Figure 1. Among the genes that were included in these regions identified by the array-CGH, the analysis by PANTHER did not show any pattern or group of genes directly related with signalling pathways involved in the CSC-like phenotype. Only a few genes were associated with pathways already described in tumorigenesis and stemness (as p53) but were not principal effectors and its position was downstream, including genes that were related with different regulation points and pathways (data not shown) making the analysis difficult. Despite these results, we noticed that several amplified areas included CpG islands ( Table 2). Using the miRBase of the University of Manchester (www.mirbase.org), we were able to identify the microRNAs included in the gained areas of the chromosomes 1, 2 and 12 (Table 3).

Gene expression profiles by microarray analysis
In order to characterize the differentially expressed genes in CSC-like compared to the standard tumor cell line, mRNA expression was analyzed by microarray analysis. After the first intensity filter, 25,368 genes out of the whole dataset of 42,405 (,60%), reached the minimum quality threshold and were selected for the differential expression analysis. When setting a threshold value of .1 fold-change, 4,831 genes (,11.4%) were identified as differentially expressed in SK-N-DZ neurospheres and 6,613 genes (,15.6%) in SIMA neurospheres. Figure 2A shows the number of genes differentially up and downregulated in each cell line. Among these gene sets, only 757 genes have significantly changed expression in both cell lines ( Figure 2B). The sets of genes up and dowregulated for each cell line are summarized in Table S1 and the genes shared by both cell lines are represented in Tables S2, S3 and S4.
Differentially expressed genes with .1 fold-change in microarray analysis were categorized according to the biological processes they are involved in, and to the molecular functions they code for, by using the DAVID database. Among all categories, those that appeared significantly altered or enriched (p,0.05 and more than 2 genes differentially expressed in both cell lines by category) were selected for further studies (Table 4).
Further analysis of the significantly enriched categories and signalling pathways involved was performed with PANTHER database and information from the NCBI BioSystems database. As a result, we observed significant alteration in representative signalling pathways known to be also involved in stem cell selfrenewal, as Wnt, Notch, Hh and TGF-b ( Figure 3).

Biological processes and signalling pathways expression in CSC-like subpopulation
The expression microarray experiment showed up a high percentage of genes altered in neurospheres (Figures 2 and 3). A complete analysis of the data highlighted some of the most relevant biological processes and molecular functions affected ( Table 4) that might be responsible for the CSC-like phenotype of this subpopulation. First, we observed an apparent alteration of ion transporters. In our analysis, several genes of Ca 2+ and K + channels showed up as differentially expressed, as KCNMA1 and CACN1AG. Secondly, we detected alterations in some of the most commonly studied signalling pathways in brain tumors and stem cell development: Wnt, Notch, Hh and TGF-b.
Analyzing Wnt signalling pathway, we observed downregulation of some of the upstream regulators of the pathway: FZD8 (Frizzled receptor), PRKCH (protein kinase C) and FRAT1, inhibitor of GSK3-mediated phosphorylation of b-catenin, and therefore a positive regulator of the pathway that could act as a proto oncogene [30,31]. The upregulation of POU5F1 as a downstream target of the pathway was also noticed.
In Notch signalling pathway, overexpression of the activator ligand of the pathway, JAG1 was identified. We also detected upregulation of HES7, one of the final targets of the pathway [32,33]. This evidence showed an apparent activation of the pathway in CSC-like comparing with the standard cell line.
Expression microarray data showed upregulation of PTHC1, the main inhibitor of Hh signalling pathway, but also one of the gene targets of the pathway which transcription would be increased if the pathway is active, as we confirmed with the RT-qPCR data ( Figure 4A).
Likewise, when analyzing TGF-b signalling pathway, the expression microarray data suggested an activation of the pathway ( Figure 4B). Ligands like TGFB3 and BMP8B were upregulated, as well as the signal transducer SHC. Quite the contrary, we found downregulation of inhibitory SMADs (like SMAD6 and SMAD7) and the pseudoreceptor BAMBI, which could sequester ligands and inhibit signal transduction by different mechanisms [34].
Parallel behaviour in SK-N-DZ and SIMA stem-like cells RT-qPCR was performed in order to corroborate the expression microarray results. In both cell lines, we confirmed the decrease in the expression of SMAD6, SMAD7 and BAMBI, as the increase observed in JAG1. The results suggested that the CSC-like subpopulation of both cell lines could have a similar expression profile after enrichment with the neurosphere assay ( Figure 5A).

Hh and TGF-b analysis in SIMA stem-like cells
Based on our previous experience in both Hh and TGF-b signalling pathways, further expression analysis by RT-qPCR was performed on some of the most relevant components of both signalling networks in SIMA cell line ( Figure 4). Results confirmed Hh aberrant activation by overexpression of its key components as PTCH1, GLI1 and SMO. Likewise, TGF-b seemed to be activated by overexpression of its receptors (TGFBR1 and TGFBR2) and ligands (TGFB2), (despite TGFB3 is downregulated even though it is also a ligand of the pathway), and downregulation of its repressors (SMAD6, SMAD7 and BAMBI).
We confirmed by RT-qPCR the increase of the levels of mRNA of JAG1, a Notch activating ligand that is related with TGF-b pathway [35]. Its upregulation supports Notch activation as aforementioned.

Expression of stem-cell markers
To confirm the enrichment of the isolated population in CSClike cells, analysis of CSC-specific cell surface markers was performed. RT-qPCR on CD133 and CD15 was carried out in both SK-N-DZ and SIMA neurospheres ( Figure 5B). On the basis of these results, SIMA culture was highly enriched in CSC-like cells as shown by the significant increase in CSC surface markers expression. An interesting point is that SK-N-DZ subpopulation did not show any overexpression of the stem-cell markers analyzed.

Discussion
The heterogeneous clinical behavior of neuroblastoma tumors has been widely studied [36][37][38][39]. It is already described that genomic rearrangements cannot explain per se this heterogeneity neither the CSC-like subpopulation [40]. In this work, we hypothesized that this heterogeneity might be due to genomic rearrangements together with changes in the expression profile of tumor cell subpopulations, particularly in CSC-like cells. In connection with this model, the genomic profile of CSC-like cells and the regulatory mechanisms or signalling pathways involved in CSC promotion could explain the heterogeneity and biology of this subpopulation and help to its identification in neuroblastoma.
Chromosomal aberrations are common in neuroblastoma, including numerical whole chromosomal gains, segmental chromosomal gains and losses and somatic mutations [2]. However, array-CGH did not show large modified areas in the genome of CSC-like subpopulations neither individually nor shared by both cell lines (Table 2). Despite this result, we identified other genetic elements involved in the biology of neuroblastoma and in the regulation of several tumors as CNV [41][42][43] and CpG islands [44][45][46] (Table 2). Both elements are being investigated as possible therapeutic targets in neuroblastoma [5] and were already identified in other tumors confirming their therapeutic potential [47][48][49][50][51]. These results suggest the fact that these two elements could be involved in tumorigenesis and CSC-like cells generation in neuroblastoma, and, as in other tumors, this might be a starting point to develop or improve current therapies. This is the case of decitabine (5-aza-29-deoxycytidine). This agent interferes with DNA methylation but the preclinical assays show that it must be administrated at doses that are not tolerable in order to produce a biological effect [52]. Our findings propose a genetic basis for the effectivity of this kind of compounds in neuroblastoma highlighting the importance for further research in demethylating agents.
In addition, some microRNAs were identified in lost and gained areas (Table 3). In the last years, these elements have become promising therapeutic targets in neuroblastoma [53]. For example, in amplified loci of SIMA neurospheres, hsa-mir-4254 and hsamir-4262 were indentified (Table 3). These two microRNAs were already related with the stem cell phenotype [54], supporting the hypothesis that the CSC might be using the same regulatory machinery as normal stem cells [10,11]. Other microRNA identified were hsa-mir-4420, hsa-mir-4429 and hsa-mir-4497, all described as possible new microRNAs involved in malignant human B cells regulation [55], and hsa-mir-3605, a novel microRNA linked to human cervical cancer [56].
In contrast, the expression array results indicated a modified expression pattern in CSC-like cells including changes in different processes and cell functions (Figures 2, 3 and Table 4). Different studies indicated some of the most important pathways involved in the development and maintenance of CSC subpopulation. Among them, the most remarkable are Wnt, Notch, Hh and TGF-b,   which have been already described for different tumors and especially in brain tumors [14]. In neuroblastoma, Wnt and Notch genes have been identified as possible CSC-like markers. The difficulty in the identification of specific genes and accurate markers of CSC-like lies in the heterogeneity of this tumor and the poor knowledge related with this subpopulation, since no specific and definitive gene expression profile has been proposed for neuroblastoma CSC-like cells [57]. Interestingly, the expression pattern in CSC-like cells involves not only Wnt and Notch but also Hh and TGF-b signalling pathways. For example, deregulation of several genes of Wnt pathway was detected. The results indicated a decrease in the expression of upstream regulators as FZD8, PRKCH and FRAT1 and an increase in the downstream target POU5F1, suggesting the inactivation of Wnt signalling pathway. POU5F1 was described in several tumors (urothelial, prostatic, cervical, breast and lung) and is implicated in different functions as chemoresistance, prognosis marker, proliferation, apoptosis, migration and invasion [58][59][60][61]. In addition to all these functions, this gene plays an important role in stemness and self-renewal of stem cells in normal and tumor tissues [60,62,63] and was found to be overexpressed and correlated with the progression of neuroblastoma [64]. This gene could interact with members from other pathways as TGF-b, Notch and Wnt [65,66] working together in the regulation of stem cell pluripotency. Our results show alterations not only in POU5F1 but also in members of TGF-b, Notch and Wnt signalling pathways indicating the potential of this gene to be used as a new therapeutical target. Our findings are supported by other works that described its potential in other tumors [62,64,67] and as a possible new CSC marker in neuroblastoma [67][68][69].
In this context, an apparent activation of Notch pathway was detected in neurospheres. It has been demonstrated that inhibition of Notch in cancer cells has the potential to slow down cell proliferation and induce apoptosis, despite the extensive crosstalk of this pathway with other major cancer pathways as Ras, Akt, NF-kB, Wnt, Hh and TGF-b [70][71][72]. This might happen because in neuroblastoma as in other tumors, the CSC subpopulation seems to be especially sensitive to inhibition of stem cell pathways as Notch [73].
Additionally, Hh was activated in CSC-like cells by upregulation of transcription targets and activators (PTCH1, SMO and GLI1) regardless downregulation of GLI2. It is known that GLI2 is a bi-functional transcription factor, with activator and repressor regions. Therefore the pathway could still be activated if GLI2 is silenced depending on the context [74,75]. In neuroblastoma CSC-like the Gli factors act through cooperative functional interactions in target gene regulation [76] although in our case it seems that GLI2 does not play an important role. Alike, the TGF-b signalling pathway is activated in neurospheres cultures ( Figure 4B) despite the downregulation of TGFB3, a ligand of the pathway. A possible explanation could be that the downregulation of BAMBI ( Figure 4B) would let a greater proportion of ligands free, increasing the concentration of extracellular ligands and consequently facilitating the activation of the pathway or acting as a negative regulator of the TGF-b pathway [34]. Interestingly, JAG1, one of the targets of TGF-b [77] and also the main activator of Notch signalling pathway is upregulated, supporting Notch signalling pathway activation. Kurpinski et al described JAG1 as a crosstalk point between these two pathways in the regulation of muscular stem cells [78] which probably suggests an important role of this gene in the generation of CSC-like cells in neuroblastoma. These findings also give us the possibility of studying targets as JAG1 that participate in several genetic pathways involved in CSC-like cells phenotype. Table 2. Cont. It is clear that the expression profile that we observed in the CSC-like cells is similar to other profiles described for this subpopulation in different tumors including neuroblastoma. The identification of CSC is based on the expression of different markers, but in the case of neuroblastoma it remains unclear which markers are accurate to identify this subpopulation, due to its heterogeneity. Our findings show differences between the two cell lines used ( Figure 5B) in two CSC markers widely described: CD133 and CD15. This result could highlight the limitations of neurospheres culture [79] since it was proved that differences in the culture conditions could affect the enrichment and isolation of CSC-like cells [80]. However, since the genes related with CSC followed a similar pattern in both cell lines, the more plausible explanation is that the phenotype of CSC-like cells is not yet fully characterized in neuroblastoma [57]. Coulon et al [57] questioned the expression of CD133 to identify CSC in neuroblastoma, highlighting the need of characterization of other CSC markers or expression profiles that allow the recognition of this subpopulation.  The choice of CSC markers is still controversial and the search of specific genes that help to identify it is the real challenge for the characterization and research of new therapeutical targets in neuroblastoma [81]. Regardless the growing evidence in new approaches and in the identification of new therapeutic targets, the treatment for neuroblastoma is still evasive. Nowadays, it includes surgery, radiation and/or chemotherapy depending on the patient's stage and risk stratification of the disease [5]. Only some new compounds are included in the schema as 13-cis-retinoic acid together with anti-GD2 antibodies and interleukin-2, a combination that has increased the progression-free survival [82]. The literature includes an increasing number of new compounds with great potential to be introduced in the treatment. This is the case of drugs that targeted MYCN and ALK inhibitors that currently are in preclinical evaluation [83]. The evidence of possible targets involved in the development and regulation of other tumors could be an approach to develop new strategies in neuroblastoma. A good example are the different pathways involved in tumorigenesis and CSC regulation, as Wnt [62,64,[67][68][69], Notch [84], Hh [85][86][87] and TGF-b [18,88,89]. In this connection, an interesting case are the ion transporters, described as essential for cell proliferation and which appear to have a role in cancer development [90][91][92]. Evidence is particularly extensive for Ca 2+ and K + channels [93] and Romania et al already described their importance in neuroblastoma [94]. In our analysis, several genes coding for Ca 2+ and K + channels were altered in CSC-like subpopulation like KCNMA1, previously described in prostate cancer [95] and in breast cancer [96,97], and CACN1AG also altered in gastric cancers, colorectal cancers and acute myeloid leukemia (AML) [98]. Several reports propose that ion transporters have clinical potential not only as therapeutical targets but also as prognosis  (**)  FGFR2, NAMPT, ACVRL1, FGFR3, ERBB3, PPARG, MITF, TGFB3,  TNFSF12, JAG1, FLT3LG, EDNRA, CD9, S1PR3, INS, TGM2, TESC,  LYN, PTPRF, SF1, PTPRU, DHRS2, TNFSF13B, ID2, DBP, CHRM1   Neuroblastoma Stem-Like Cells PLOS ONE | www.plosone.org markers or to improve actual treatments [99][100][101][102][103], as blocking channel activity seems to impair the growth of some tumors, including neuroblastoma [99,104]. The presence of altered genes of Ca 2+ and K + channels in neuroblastoma CSC-like subpopulations confirms and opens the possibility to develop new compounds that could act at different levels and subpopulations within neuroblastoma. This is only one example of the importance of redirecting treatments towards more effective CSC molecular therapies to achieve total elimination of the tumor cell population and avoid relapses [40,105].
In conclusion, we found that the expression profile of neuroblastoma CSC-like cells differs significantly from the expression profile in standard neuroblastoma cell lines, which is seen as altered signalling pathways involved in stem cell proliferation pathways as Wnt, Notch, Hh and TGF-b, suggesting a cross talk among them and with other pathways. The CSC markers analysis reveals no correlation between the cell lines, displaying the heterogeneity of neuroblastoma tumors and CSC-like population but opening the possibility to identify other potential markers. Taking together, these results confirm the importance of different pathways in CSC regulation and the identification of possible candidate targets for molecular CSC therapies in neuroblastoma. The results presented highlight new information for the CSC-like phenotype in neuroblastoma and indicate the importance of redirecting current cancer treatments towards CSC molecular therapies to achieve total elimination of the tumor cell population and improve treatment effectiveness.