Bioinformatics analysis of the molecular mechanism of diffuse intrinsic pontine glioma.

The present study aimed to elucidate key molecular mechanisms in the progression of diffuse intrinsic pontine glioma (DIPG). The gene expression profile GSE50021, which consisted of 35 pediatric DIPG samples and 10 normal brain samples, was downloaded from the Gene Expression Omnibus database. The differentially-expressed genes (DEGs) in the pediatric DIPG samples were identified. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathways of DEGs were enriched and analyzed. The protein-protein interaction (PPI) network of the DEGs was constructed and functional modules of the PPI network were disclosed using ClusterONE. A total of 679 DEGs (454 up- and 225 downregulated) were identified in the pediatric DIPG samples. DEGs were significantly enriched in various GO terms, and KEGG and Reactome pathways. The PPI network of upregulated (153 nodes and 298 connections) and downregulated (71 nodes and 124 connections) DEGs, and two crucial modules, were obtained. Downregulated genes in module 2, such as cholecystokinin (CCK), gastrin (GAST), adenylate cyclase 2 (brain) (ADCY2) and 5-hydroxytryptamine (serotonin) receptor 7 (HTR7), were significantly enriched in the calcium signaling pathway, the neuroactive ligand-receptor interaction pathway and in GO terms, such as the G-protein coupled receptor (GPCR) signaling pathway, while upregulated genes in module 1 were not enriched in any pathways or GO terms. CCK and GAST associated with the GPCR signaling pathway, HTR7 enriched in the neuroactive ligand-receptor interaction, and ADCY2 and HTR7 involved in the calcium signaling pathway may be key mechanisms playing crucial roles in the development and progression of DIPG.

Introduction

Diffuse intrinsic pontine glioma (DIPG) is a fatal malignant pediatric brain tumor and the leading cause of cancer-related mortality in children (1). Although collaborative efforts have been made to improve treatments, survival times have remained static over the past three decades, and the vast majority of children succumb within 2 years of diagnosis (2). Radiotherapy is the standard treatment offering a significant but transient improvement, while chemotherapy has not shown any benefit thus far (1,3). Due to a lack of therapeutic benefits from biopsies and the lack of available tumor material, molecular studies of DIPG have been greatly limited, thus hampering the development of targeted therapies for the disease (4). Therefore, an improved understanding of the molecular mechanism of DIPG is required to develop more specifically therapeutic approaches for this disease.

The biological properties of pediatric DIPG are not identical to those of cerebral high-grade gliomas of adults, which contribute to the failure of the majority of therapeutic trials (3). Analysis of the molecular mechanism of DIPG is useful for developing an effective therapy. Previous studies uncovered the unique genetic composition of this pediatric brain cancer, with nearly 78% found to possess frequent histone H3 gene (encoding p.Lys27Met) mutations (5). Additionally, Buczkowicz et al performed genomic analysis of DIPG and revealed a novel recurrent activating type I activin A receptor (ACVR1) mutation (6). Mutations in ACVR1 will increase the expression of activin signaling targets, inhibitor of DNA binding 1 (ID1) and ID2, which are known targets of bone morphogenetic protein (BMP) signaling (7–9). In tumors, BMP signaling drives either a differentiation or proliferation response in glioblastoma stem cells, which is to a certain extent associated with the epigenetic state of the cell (10,11). Contrasting effects are associated with BMP signaling depending on the context. These context-dependent consequences may further precipitate the exclusive association between the ACVR1 mutation and DIPG (12,13). Furthermore, whole-genome profiling of pediatric DIPG performed by Zarghooni et al (14) highlighted the roles of platelet-derived growth factor receptor and poly (ADP-ribose) polymerase in DIPG progression. Despite expanded efforts to investigate the genetic basis of DIPG, to date, the molecular mechanisms of its development and progression remain unclear.

In a previous study, GSE50021 was used only as expression data for investigating what drives DIPGs via integration of whole-genome sequencing with methylation and copy number profiling (6). There was no systematic and comprehensive analysis for this expression profile. In the current study, microarray analysis was performed for screening the differentially-expressed genes (DEGs) in DIPG samples, then hub genes, as well as the significant modules and pathways underlying the onset and progression of DIPG, were identified using a comprehensive bioinformatics methods. The study aimed to identify candidate genes and associated pathways of DIPG to clarify its molecular mechanisms.

Materials and methods

Source of data

The GSE50021 microarray expression profile was extracted from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) database based on the GPL13938 Illumina HumanHT-12 WG-DASL V4.0 expression beadchip platform. The gene expression files were deposited by Buczkowicz et al (6). A total of 35 pediatric DIPG samples and 10 normal brain samples were available for the expression array.

After obtaining the raw data, the robust multiarray average (15) method of the Linear Models for Microarray Data (Limma) package (16) in Bioconductor was used to perform background correction, quartile normalization and probe summarization. Finally, the gene expression matrix of the samples was acquired.

Identification of DEGs

As the DEGs were associated with the development of the disease, a t-test in Limma package (16) was used to identify the DEGs in pediatric DIPG samples compared with the normal controls. All genes with a P-value of <0.01 and |log2fold change (FC)|>2.0 were selected as DEGs for further study.

Functional enrichment analysis

Gene Ontology (GO) (17) analysis is applied for functional studies of large-scale genomic or transcriptomic data. The Kyoto Encyclopedia of Genes and Genomes (KEGG) (18) and Reactome (19) are the major recognized pathway-related databases that record information on how genes are networked in various organisms. A functional enrichment analysis was performed by the Database for Annotation, Visualization and Integrated Discovery (DAVID) (20), which is the online tool for systematically extracting biological meaning and functional annotation behind large gene or protein lists.

GO categories, including biological process (BP), molecular function (MF), and cellular component (CC), of the identified DEGs were enriched using DAVID. Also, KEGG and Reactome pathway enrichment analyses of the DEGs were applied for the identification of the main functional and metabolic pathways involved in DIPG using DAVID. A P-value of <0.01 was chosen as the significant threshold.

Functional annotation of DEGs

TSGene [the tumor suppressor gene (TSG) database] (21) is used for further investigation of the molecular mechanisms of TSGs in cancer, and provides detailed annotations for each TSG, such as transcription factor (TF) regulations. The Tumor-Associated Gene (TAG) (22) database consists of five tables, each of which summarizes attributes for a specific entity associated with the TAG.

Based on the data information of TFs, the functional enrichment of the DEGs for transcriptional regulation was assessed and annotated according to the GO and KEGG pathway annotation terms. In addition, the DEGs were mapped into the TSGene and TAG databases to screen the known oncogene and TSGs.

Protein-protein interaction (PPI) network construction and module selection

The PPI network is increasingly used for identifying the associations between proteins in various organisms. The Search Tool for the Retrieval of Interacting Genes (STRING) (23) is an online database providing experimental and predicted PPI information. Cytoscape (24) is a biological graph visualization tool for integrated models of biomolecular interaction networks. In the present study, the STRING online tool was applied to analyze the PPIs of the DEGs with a combined score of >0.4 as the cut-off value, then the PPI networks of the up- and downregulated genes were constructed respectively using Cytoscape software.

In addition, based on connectivity degree, the enriched functional modules of the PPI network were disclosed using the ClusterONE Cytoscape plug-in (25). P<0.01 was chosen as the significance threshold. Furthermore, GO and KEGG pathway enrichment analyses of the DEGs in the significant modules were performed using DAVID to analyze the gene function at the molecular level.

Results

According to t-test analysis of the DEGs in the pediatric DIPG samples compared with the normal controls, 679 DEGs were recorded, including 454 upregulated and 225 downregulated genes. The results showed that the number of upregulated genes was more than the number of downregulated genes.

In the present study, 8 and 6 KEGG pathways corresponding to the up- and downregulated genes, respectively, were enriched using the cutoff criterion. Results of the KEGG pathway analysis showed that the enriched pathways of the upregulated genes included the ribosome, focal adhesion and extracellular matrix-receptor interaction pathways (Table I), while those enriched by the downregulated genes included the calcium signaling, cytokine-cytokine receptor interaction and neuroactive ligand-receptor interaction pathways (Table I).

Table I.

Enriched KEGG pathways of differentially-expressed genes.

KEGG pathways	Gene count	P-value
Upregulated genes
Ribosome	12	1.45×10−5
Focal adhesion	18	2.21×10−5
ECM-receptor interaction	10	1.85×10−4
Shigellosis	8	3.91×10−4
Adherens junction	8	1.32×10−3
Regulation of actin cytoskeleton	14	4.10×10−3
Pathogenic Escherichia coli infection	6	5.92×10−3
Colorectal cancer	6	9.68×10−3
Downregulated genes
Calcium signaling pathway	8	6.10×10−4
Cytokine-cytokine receptor interaction	9	2.11×10−3
Neuroactive ligand-receptor interaction	9	2.52×10−3
Dilated cardiomyopathy	5	2.84×10−3
Arachidonic acid metabolism	4	3.73×10−3
Pancreatic secretion	5	4.66×10−3

Gene count represents the number of genes enriched in the KEGG pathway. KEGG, Kyoto Encylopedia of Genes and Genomes; ECM, extracellular matrix.

In addition, based on the analysis of Reactome pathways, it was found that upregulated genes were enriched in 47 pathways, while downregulated genes were enriched in 16 pathways. The enriched Reactome pathways of upregulated genes were mainly involved in translation, metabolism of proteins and folding of actin by CCT/TriC. Downregulated genes were mainly enriched in G-protein coupled receptor (GPCR) ligand binding, signaling by GPCR and class A/1 (rhodopsin-like receptors).

Furthermore, GO terms in different GO categories (BP, MF and CC) were also enriched by the selected DEGs. Results showed that the upregulated genes were significantly associated with cellular membrane organization, membrane organization and protein targeting to the endoplasmic reticulum (Table II), while the downregulated genes were mainly associated with single-multicellular organismal process, plasma membrane part and neuropeptide hormone activity (Table II).

Table II.

Top 15 significantly enriched GO terms of differentially-expressed genes.

GO ID	Term	Gene count	P-value
Upregulated genes
GO:0016044	Cellular membrane organization	46	4.96×10−9
GO:0061024	Membrane organization	46	5.42×10−9
GO:0045047	Protein targeting to ER	16	1.84×10−8
GO:0072599	Establishment of protein localization to ER	16	2.11×10−8
GO:0016043	Cellular component organization	165	4.61×10−8
GO:0005515	Protein binding	230	1.61×10−7
GO:0003735	Structural constituent of ribosome	15	5.03×10−6
GO:0019900	Kinase binding	24	7.26×10−5
GO:0030276	Clathrin binding	5	9.01×10−5
GO:0032403	Protein complex binding	20	1.77×10−4
GO:0044424	Intracellular part	360	4.02×10−13
GO:0005829	Cytosol	117	1.01×10−12
GO:0005622	Intracellular	361	2.00×10−12
GO:0043226	Organelle	325	8.83×10−12
GO:0005737	Cytoplasm	291	1.22×10−11
Downregulated genes
GO:0044707	Single-multicellular organismal process	98	1.13×10−7
GO:0032501	Multicellular organismal process	99	3.81×10−7
GO:0003008	System process	43	6.41×10−7
GO:0050877	Neurological system process	35	6.61×10−7
GO:0019226	Transmission of nerve impulse	24	5.55×10−6
GO:0044459	Plasma membrane part	46	1.10×10−6
GO:0031226	Intrinsic to plasma membrane	34	3.72×10−6
GO:0044421	Extracellular region part	30	1.45×10−5
GO:0005886	Plasma membrane	76	2.68×10−5
GO:0071944	Cell periphery	77	2.95×10−5
GO:0005184	Neuropeptide hormone activity	5	4.34×10−6
GO:0005179	Hormone activity	8	4.59×10−5
GO:0072509	Divalent inorganic cation transmembrane transporter activity	9	4.67×10−5
GO:0015085	Calcium ion transmembrane transporter activity	8	6.66×10−5
GO:0005262	Calcium channel activity	7	1.77×10−4

Gene count represents the number of genes enriched in the GO term. GO ID, Gene Ontology identification; ER, endoplasmic reticulum.

Based on the TF and TAG data, the results showed that upregulated functional genes contained 20 TFs, including c-abl oncogene 1, non-receptor tyrosine kinase, actin-like 6A and aryl hydrocarbon receptor, and 32 TAGs (8 oncogenes, 15 TSGs and 9 unknown genes associated with cancer progression). Downregulated functional genes contained 7 TFs, including homeobox A7, MEIS homeobox 2 and nescient helix loop helix 2, and 12 TAGs (3 oncogenes, 5 TSGs and 4 unknown genes associated with cancer progression).

PPI network construction and module selection

According to the STRING data, PPI networks of the upregulated (153 nodes and 298 connections) and downregulated (71 nodes and 124 connections) DEGs were obtained (Fig. 1). Upregulated genes with higher node degree mainly were chaperonin-containing TCP1, subunit 2 (β) (CCT2), integrin, β1 (ITGB1), ITGB4, CD44 molecule (Indian blood group), FYN oncogene-related to SRC, FGR, YES (FYN), and cyclin D1 (CCND1) (Fig. 1A). Downregulated genes with higher node degree were as follows: Adenylate cyclase 2 (brain) (ADCY2), prodynorphin (PDYN), cholecystokinin (CCK) and chemokine (C-C motif) receptor 7 (CCR7) (Fig. 1B). The node degrees were >10 for those genes.

Figure 1.

Protein-protein interaction networks of differentially-expressed genes. The red and green nodes indicate (A) upregulated and (B) downregulated genes, respectively. The edges/lines indicate interactions between these genes.

Functional module analysis

In the present study, the upregulated genes were enriched in 23 significant modules, while the downregulated genes were enriched in 15 significant modules by ClusterONE. The most significant modules enriched by the up- and downregulated genes, respectively, were module 1 (Fig. 2A) and module 2 (Fig. 2B). Module 1 had 7 nodes and 16 edges, while module 2 had 14 nodes and 48 edges.

Figure 2.

Functional modules of differentially-expressed genes. (A) Module 1 and (B) module 2. The red nodes indicate upregulated genes and green nodes indicate downregulated genes. The edges/lines indicate interactions between these genes.

Based on the connectivity degree, the hub proteins in module 1 were centrosomal protein 63 kDa (CEP63), CDK5 regulatory subunit associated protein 2 (CDK5RAP2), dynein, cytoplasmic 1, intermediate chain 2 (DYNC1I2), centrin, EF-hand protein, 2 (CETN2) and centrosomal protein 152 kDa (CEP152). The hub proteins in module 2 included ADCY2, CCK, 5-hydroxytryptamine (serotonin) receptor 7 (HTR7), gastrin-releasing peptide, gastrin (GAST) and vasoactive intestinal peptide (VIP).

In addition, the upregulated genes in module 1 were not enriched in KEGG pathways and GO terms, while the downregulated genes in module 2 were significantly enriched in 3 KEGG pathways and 51 GO-BP terms. The pathways enriched by the downregulated genes in module 2 were the calcium signaling, neuroactive ligand-receptor interaction and vascular smooth muscle contraction pathways, and the enriched GO terms mainly were the GPCR signaling pathway, cell surface receptor-linked signal transduction and positive regulation of the homeostatic process.

Discussion

DIPG represents the largest challenge in the field of pediatric neuro-oncology, with no effective treatment and a near 100% fatality rate (6). Additionally, a paucity of DIPG tissue is available for molecular studies (26). Thus, the molecular mechanism underlying DIPG progression remains unclear. In the present study, 679 DEGs (454 up- and 225 downregulated) were identified in pediatric DIPG samples compared with normal controls. Based on the functional modules analysis, the upregulated genes in module 1 were not enriched in any GO terms and pathways, while the downregulated genes in module 2 were significantly enriched in various GO terms and pathways; for example, CCK and GAST were associated with the GPCR signaling pathway, HTR7 was significantly enriched in neuroactive ligand-receptor interaction, and ADCY2 and HTR7 were enriched in the calcium signaling pathway. Thus, the aforementioned downregulated genes and associated pathways merit further consideration.

GAST and CCK were hub proteins with higher node degrees in module 2. GAST and CCK together constitute a family of homologous peptide hormones, which are potent neurotransmitters in the brain and the periphery (27). GAST has been found to function to endogenously modulate the progression of glioma and the invasion to surrounding structures (28). GAST reduces the invasive abilities of gliomas cells in vitro (28), thus, downregulated GAST may enhance the migratory abilities of glioma cells. In addition, CCK is a gut-brain peptide that has been confirmed to be able to induce mitosis in human gliomas (29). A previous study established that CCK-8s and CCK-33 directly stimulate human glioma cell growth and invasion process (28). Also, CCK is downregulated in glioblastoma compared with normal human brain tissues and may function as a TSG in human glioma (30). By contrast, the present study found that GAST and CCK were significantly associated with GPCR signaling pathway. G protein can couple with and activate adenylate cyclase to increase the concentrations of cyclic adenosine monophosphate (cAMP) through the stimulatory GPCR signaling pathways (31). cAMP/protein kinase A (PKA) pathway activation can induce glioma cell differentiation (32). Previous studies have also confirmed that mature amidated GAST peptides act through the G-protein-coupled CCK-2/GAST receptor, and have been shown to play a role in malignancy (33,34). Therefore, the present results are in line with previous findings and suggest that GAST and CCK may play roles in DIPG progression via the GPCR signaling pathway.

HTR7 was a hub protein significantly associated with neuroactive ligand-receptor interactions in the present study. The HTR7 gene has been localized to chromosome 10 (10q21-24), and the loss of 10q underpins the progression of DIPG (35,36). Additionally, a study by Chan et al verified that neuroactive ligand-receptor interactions are strongly associated with DIPG tumorigenesis (37). Thus, the present results are consistent with previous findings and suggest that HTR7 may play roles in DIPG tumorigenesis via neuroactive ligand-receptor interactions.

Furthermore, HTR7 was not only enriched in neuroactive ligand-receptor interactions, but also in the calcium signaling pathway. Additionally, ADCY2 was a hub protein significantly associated with the calcium signaling pathway. Calcium is an intracellular signaling molecule participating in tumor cell growth and metastasis (38). Calcium signaling has significant functions in numerous signaling processes in the proliferation and motility of glioblastoma cells (39). Additionally, ADCY2 encodes a member of the family of adenylate cyclases (40) and the activity of HTR7 receptor mediated by G proteins can stimulate adenylate cyclase (41). Adenylate cyclases catalyze the formation of cAMP (42). cAMP can trigger the calcium signaling pathway and be involved in glioma cell differentiation via the cAMP/PKA pathway (32,43). Although the roles of HTR7 and ADCY2 in the progression of DIPG have not been fully discussed, we hypothesize that HTR7 and ADCY2 could play roles in DIPG tumorigenesis via the calcium signaling pathway mediated by cAMP. Further genetic studies are required to confirm this supposition.

Taken together, the results of the present study indicated that CCK and GAST in association with the GPCR signaling pathway, HTR7 enriched in neuroactive ligand-receptor interaction, and ADCY2 and HTR7 involvement in the calcium signaling pathway may be key mechanisms playing crucial roles in the development and progression of DIPG. The present findings provide novel insights into the molecular mechanism of DIPG progression and suggest directions for future research. A better understanding of the molecular mechanism of the pathogenesis of DIPG will provide a broader perspective for further clinical molecular target therapy experiments. However, as no experiments were performed in the present study, experimental validation is required to confirm these observations.