This article is mentioned in:

Introduction

Neuroblastoma (NB) is a cancer of the peripheral sympathetic neurons, which occurs during childhood, and the percentage of NB of all pediatric oncology-associated mortality is ~10% worldwide (1). NB predominately occurs in young children, of which the average age is 17 months (2) NB derives from undifferentiated neural crest cells and metastasizes to other organs (3). Patients with NB have poor prognosis with a 5-year survival rate of no more than 30% (4).

With the development of the field of biology, increasing numbers of markers correlated with NB have been identified. MYCN is a well-known prognostic marker of NB (5), and MYCN amplification is the most significant molecular marker of risk in NB (6). Tang et al demonstrated that patients with NB exhibiting low expression levels of MYCN and TrkA have a 5-year survival rate of 63.7%, whereas the 5-year survival rate of patients with high expression levels of the two is 88.1% (7). In addition, the alternative TrkAIII splice variant is involved in the pathogenesis of advanced stage human NB through microtubules, which are involved in the promotion and maintenance of NB cells (8). In addition to these reports, NCYM has been observed to act as a cis-antisense gene of MYCN, and contributes to the stabilization of MYCN in human NB by inhibiting GSK3β (9). NCYM, co-amplified with MYCN, is involved in the pathogenesis of NB (10). Although several studies have reported that MYCN is key in the pathogenesis of NB, the genes correlated with MYCN amplification remain to be fully elucidated. Therefore, the present study aimed to examine the molecular mechanism of MYCN amplification in NB.

Wu et al analyzed expression profiles using Ingenuity Pathway Analysis, and found that the aryl hydrocarbon receptor downregulates the expression of MYCN through the mediation of E2F1 in the protein-protein interaction (PPI) network, and further confirmed that the aryl hydrocarbon receptor regulates the activity of the MYCN promoter and results in the downregulated expression of MYCN (11). To identify more genes correlated with MYCN in the present study, the determination of differentially expressed genes (DEGs) between NB cell lines exhibiting MYCN amplification (MYCN amplification group) and NB cell lines with a normal MYCN copy number (control group) were identified. Hierarchical clustering and Gene Ontology (GO) analysis were performed for these DEGs, and a protein-protein interaction (PPI) network was constructed for the DEGs to identify the key genes associated with MYCN. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was also performed for the DEGs in PPI network. The correlation of MYCN and the key gene associated with to MYCN was determined using Pearson's correlation coefficient, in order to investigate whether their association was synergistic or antagonistic.

Materials and methods

Microarray data

Microarray data (accession. no. GSE53371) were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) (11). The microarray platform of GSE53371 was GPL887 Agilent-012097 Human 1A Microarray (V2) G4110B (Feature Number version; Agilent Technologies, Inc., Santa Clara, CA, USA). A total of 20 samples were available, which included 10 NB cell lines with MYCN amplification and 10 NB cell lines with a normal MYCN copy number.

Data preprocessing and identification of DEGs

Normalization of the microarray data were performed using a median normalization method (12). The DEGs were identified using the Linear Models for Microarray Data (Limma) package (13). The raw P-value was adjusted into a false discovery rate (FDR) using the Benjamini & Hochberg method (14,15). The genes with an FDR <0.05 and |log2FC (fold change)|>1 were considered significantly different between the MYCN amplification group and the control group.

Hierarchical clustering analysis of DEGs

Based on the Euclidean distance, which is the actual distance between two points that may be calculated using the Pythagorean formula, hierarchical clustering analysis (16) was performed to evaluate the sample specificity. In the hierarchical clustering analysis, the expression value of the DEGs were extracted, and subsequent analysis was performed of the DEGs using the Pheatmap package in R language (http://cran.r-project.org/web/packages/pheatmap/index.html) (16,17).

Functional enrichment analysis of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID) is a comprehensive functional annotation tool (18). Gene Ontology (GO) analysis was performed for the DEGs using the DAVID database (http://david.abcc.ncifcrf.gov/). The raw P-value was adjusted into the FDR using Benjamini & Hochberg's method and an FDR <0.05 was selected as the cut-off criterion.

Construction of the PPI network

The Search Tool for the Retrieval of Interacting Genes (STRING) database is used to collect and predict PPI associations (19). In the present study, a PPI network was constructed for the DEGs using STRING, and a confidence score >0.4 was selected as the cut-off criterion. Subsequent visualization of the PPI network was performed using Cytoscape (http://cytoscape.org/) (20).

Pathway analysis of DEGs in the PPI network

The WEB-based GEne SeT AnaLysis Toolkit (WebGestalt) aims to examine large sets of genes, and is used for functional enrichment analysis, including GO enrichment, pathway enrichment and transcription factor analysis, by integrating functional categories (21,22). KEGG (http://www.genome.ad.jp/kegg) pathway enrichment analysis was performed for the DEGs in the PPI network using WebGestalt, and P<0.05 was selected as the cut-off criterion.

Correlation analysis

Based on the PPI network, the expression values of the genes correlated directly with MYCN were extracted, and correlation analysis was performed using the expression values of the genes. For this correlation analysis, the correlation of MYCN and the key gene associated with MYCN was determined using Pearson's correlation coefficient (23) with SPSS 13.0 (SPSS Inc., Chicago, IL, USA).

Results

Identification of DEGs

In total, 172 DEGs between the MYCN amplification group and the control group were identified, which included 137 downregulated DEGs and 35 upregulated DEGs.

Hierarchical clustering analysis of the DEGs

Following hierarchical clustering analysis, there were two distinct gene clusters identified, including the MYCN amplification cluster and control cluster (Fig. 1). These results suggested that there were two different gene expression patterns with marked color differences, which were available for use to distinguish between the MYCN amplification samples and the control samples.

Figure 1 Hierarchical clustering analysis of differentially expressed genes. Red represents the higher levels of gene expression; blue represents lower levels of gene expression.

GO enrichment analysis of the DEGs

GO enrichment analysis revealed that the five enriched terms were as follows: Regulation of action potential in neuron (P=1.24E-05), myelination (P=2.95E-05), regulation of action potential (P=4.71E-05), ensheathment of neurons (P=4.85E-05) and axon ensheathment (P=4.85E-05), as listed in Table I. The KCNMB4, PLP1, CLDN1, LGI4, MAL, PMP22 and GAL3ST1 genes were involved in the regulation of action potential in neuron biological process, and the percentage of the genes enriched in this term was 21.88% (Fig. 2).

Figure 2 Pie chart of Gene Ontology enrichment analysis for the differentially expressed genes. The percentage was proportional to the differentially expressed genes enriched in each function.

Table I Top five most enriched GO terms for the differentially expressed genes.

Table I

Top five most enriched GO terms for the differentially expressed genes.

GO ID	GO term	Gene count	P-value	FDR	Gene
GO:0019228	Regulation of action potential in neuron	7	1.24E-05	0.015898	KCNMB4, PLP1, CLDN1, LGI4, MAL, PMP22, GAL3ST1
GO:0042552	Myelination	6	2.95E-05	0.018878	PLP1, CLDN1, LGI4, MAL, PMP22, GAL3ST1
GO:0001508	Regulation of action potential	7	4.71E-05	0.020080	KCNMB4, PLP1, CLDN1, LGI4, MAL, PMP22, GAL3ST1
GO:0007272	Ensheathment of neurons	6	4.85E-05	0.015539	PLP1, CLDN1, LGI4, MAL, PMP22, GAL3ST1
GO:0008366	Axon ensheathment	6	4.85E-05	0.015539	PLP1, CLDN1, LGI4, MAL, PMP22, GAL3ST1

[i] GO, Gene Ontology; FDR, false discovery rate.

PPI network analysis of DEGs and pathway analysis

Analysis of the PPI network revealed that there were 85 nodes and 142 gene links, and SOX2 was identified as the key gene correlated directly with MYCN (Fig. 3).

Figure 3 Protein-protein interaction network of differentially expressed genes. The blue square node represents MYCN; the green square node represents SOX2. The line between two genes represents their interaction.

KEGG pathway enrichment analysis of the DEGs in the PPI network indicated that the ASPA, MAOA and ALDH1A3 genes were involved in the histidine metabolism pathway (P=0.012376), and the LAMA2, FOS, ITGA6, GLI3, GLI1 and RAD51 genes were involved in the extracellular matrix (ECM)-receptor interaction pathway (P=0.040716; Table II).

Table II Top two most enriched KEGG pathways of the differentially expressed genes in the protein-protein interaction network.

Table II

Top two most enriched KEGG pathways of the differentially expressed genes in the protein-protein interaction network.

ID	KEGG pathway	Gene count	P-value	Gene
hsa00340	Histidine metabolism	3	0.012376	ASPA, MAOA, ALDH1A3
hsa04512	ECM-receptor interaction	6	0.040716	LAMA2, FOS, ITGA6, GLI3, GLI1, RAD51

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix.

Correlation analysis of genes

Correlation analysis indicated that the expression of SOX2 was inversely correlated with the expression of MYCN in the NB cell lines, and the correlation coefficient of SOX2 and MYCN was −0.83 (Fig. 4).

Figure 4 Correlation analysis of the expression levels of SOX2 and MYCN. Correlation coefficient= −0.83. Values were obtained relative to the control.

Discussion

In the present study, 172 DEGs between the MYCN amplification group and control group were identified, which included 137 downregulated and 35 upregulated DEGs. Functional enrichment analysis indicated that the KCNMB4 gene was involved in the regulation of action potential in neuron pethway, and the FOS, GLI3 and GLI1 genes were involved in the ECM-receptor interaction pathway. Correlation analysis demonstrated that the expression of SOX2 was inversely correlated with the expression of MYCN in the NB cell lines.

Neurons encode and convey information by generating action potentials (24). According to a previous report by Hall et al, the expression of specific Kv3.1 glycoprotein, which is a voltage-gated potassium (Kv) channel, has an impact on the wave form of action potentials and N-glycans associated with the Kv3.1 protein affect the migratory rate of The NB cell (25). In addition, Leung et al demonstrated that the voltage-gated Kv channel has an effect on cAMP-stimulated neuritogenesis in mouse NB N2A cells (26). These findings suggest that the voltage-gated Kv channels are important in the development of NB through the modulation of action potentials. Large-conductance Ca2+-activated K+ (BK) channels contribute to the excitability of neurons (27). The two subunits of BKβ are BKβ3 (KCNMB3) and BKβ4 (KCNMB4) (28). Therefore, KCNMB4 may be involved in the pathogenesis of NB through the voltage-gated Kv channel. In the present study, KCNMB4 was involved in the regulation of action potential in neuron biological process. Based on these findings, it was hypothesized that KCNMB4 is involved in the development of NB by the regulation of action potentials in neurons.

SOX2, located in chromosome 3q26.33, is a transcription factor which is involved in the regulation of stem cell properties (29,30). Previously, it has been reported that neural stem cell-like cells are isolated from certain types of neural cancer, including glioblastoma, medulloblastoma and NB (31,32). In addition, Das et al demonstrated that SOX2 is involved in the maintenance of undifferentiated stem cells by the regulation of all-trans retinoic acid, which induces the differentiation of NB cells. Furthermore, miR-340 has been found to be involved in the pathogenesis of NB by mediating the expression of SOX2 (33). Wang et al found that FoxM1 is involved in the tumorigenicity of aggressive NB cells by activating the expression of SOX2 (34). According to a report by Yang et al, the expression of OCT4 and SOX2 are involved in the progression of NB (35). These results indicate that SOX2 is important in the pathogenesis of NB. The results of the present study, which revealed that the expression of SOX2 was directly correlated with the expression of MYCN, were consistent with previous findings. Therefore, SOX2 may be involved in the pathology of NB through the regulation of MYCN.

In the present study, PPI network analysis indicated that GLI1 regulated SOX2 by the mediation of GLI3, and FOX also regulated the expression of SOX2 by the mediation of DCN. GLI1 and GLI3 are mediators of the SHH pathway, which has an effect on the early development of the central nervous system (36). In addition, SHH-GLI1 is involved in the self-renewal of cancer stem cells (37). Furthermore, Shahi et al reported that the expression of SMO and GLI3 are partially correlated in NB (38). These results suggest that GLI1 and GLI3 may be involved in the development of NB by the mediation of SOX2. In addition, FOX proteins are a conserved transcriptional regulator superfamily, and FOX family transcription factors are key in the progression of cancer (39). In NB, FOXR1, fused to MLL or PAFAH1B due to interstitial deletions, function as oncogenes (40). Santo et al demonstrated that events of intra-chromosomal deletion/fusion at 11q23 activate the expression of FOXR1 in NB (41). This suggested that FOX may also be involved in the pathogenesis of NB by the regulation of SOX2. According to these findings, the present study hypothesized that FOS, GLI3 and GLI1 may regulate the expression of MYCN by the mediation of SOX2. Although the findings of the present study were consistent with previous observations that the FOS, GLI3 and GLI1 genes are involved in the pathogenesis of NB, the present study revealed that the FOS, GLI3 and GLI1 genes were involved in the ECM-receptor interaction pathway. Therefore, FOS, GLI3 and GLI1 may be involved in the ECM-receptor interaction pathway to contribute to the development of NB.

In conclusion, the present study identified 172 DEGs between the MYCN amplification group and the control group, including 137 downregulated and 35 upregulated DEGs. KCNMB4 may contribute to the development of NB by the regulation of action potentials in neurons. FOS, GLI3 and GLI1 may be involved in the pathogenesis of NB by the regulation of SOX2, the expression of which was inversely correlated to the expression of MYCN in the NB cell lines. These findings have provided a novel opportunity to investigate the pathogenesis of NB, however, these results require confirmation by further investigations.

References