Identification of key genes in glioma CpG island methylator phenotype via network analysis of gene expression data.

Gene expression data were analysed using bioinformatic tools to demonstrate molecular mechanisms underlying the glioma CpG island methylator phenotype (CIMP). A gene expression data set (accession no. GSE30336) was downloaded from Gene Expression Omnibus, including 36 CIMP+ and 16 CIMP- glioma samples. Differential analysis was performed for CIMP+ vs. CIMP‑ samples using the limma package in R. Functional enrichment analysis was subsequently conducted for differentially expressed genes (DEGs) using Database for Annotation, Visualization and Integration Discovery. Protein‑protein interaction (PPI) networks were constructed for upregulated and downregulated genes with information from STRING. MicroRNAs (miRNAs) targeting DEGs were also predicted using WebGestalt. A total of 439 DEGs were identified, including 214 upregulated and 198 downregulated genes. The upregulated genes were involved in extracellular matrix organisation, defence and immune response, collagen fibril organisation and regulation of cell motion and the downregulated genes in cell adhesion, sensory organ development, regulation of system process, neuron differentiation and membrane organisation. A PPI network containing 134 nodes and 314 edges was constructed from the upregulated genes, whereas a PPI network consisting of 85 nodes and 80 edges was obtained from the downregulated genes. miRNAs regulating upregulated and downregulated genes were predicted, including miRNA‑124a and miRNA‑34a. Numerous key genes associated with glioma CIMP were identified in the present study. These findings may advance the understanding of glioma and facilitate the development of appropriate therapies.

Introduction

Malignant glioma is the most common central nervous system tumour in adults and is associated with significant morbidity and mortality (1). Gliomas are highly invasive and poorly respond to conventional treatments; therefore, further studies to support the development of therapy for them are warranted (2).

Alterations in methylation serve a critical role in the pathogenesis of numerous human malignancies, including gliomas (3). CpG island methylator phenotype (CIMP) has emerged as a distinct molecular subclass of tumours (4). It features extensive, coordinated hypermethylation at specific loci (5,6). Several key genes regulated by methylation have been previously identified. O6-methylguanine-DNA methyltransferase (MGMT), which is responsible for DNA repair, is associated with chemotherapy resistance (7). Previous studies indicated that epigenetic silencing of MGMT via promoter methylation serves an important role in the regulation of MGMT expression in gliomas (8). Bruna et al (9) demonstrated that the methylation of platelet-derived growth factor (PDGF)-B can dictate transforming growth factor-β as an oncogenic factor to promote cell proliferation in human glioma. In addition, Wiencke et al (10) reported that methylation of the phosphatase and tensin homolog promoter defines low-grade gliomas and secondary glioblastoma. Mueller et al (11) also suggested that epigenetic dysregulation of runt-related transcription factor 3 and testin is involved in glioblastoma tumorigenesis. Abnormal DNA methylation of CD133 (12) and tumor protein 53 (13) is also observed in glioma. Additionally, Turcan et al (14) indicated that isocitrate dehydrogenase 1 mutation is sufficient to establish the glioma hypermethylator phenotype. However, the identification of glioma-CIMP (G-CIMP) tumours based on gene expression data has rarely been reported (15). In the present study, gene expression profiles of CIMP-positive (CIMP+) samples were compared with those of CIMP-negative (CIMP−) samples to identify differentially expressed genes (DEGs), which were further subjected to functional enrichment analysis and network analyses. The findings of the present study may extend the understanding of the molecular mechanisms of CIMP+ glioma.

Materials and methods

Gene expression data

A gene expression data set (accession no. GSE30336) was downloaded from Gene Expression Omnibus (14), including 36 CIMP+ glioma and 16 CIMP− samples. Gene expression levels were measured using the GPL571 (HG-U133A_2) Affymetrix Human Genome U133A 2.0 Array (Affymetrix; Thermo Fisher Scientific Inc., Waltham, MA, USA). Probe annotations were also acquired.

Pretreatment and differential analysis

Raw data were pre-treated with the Robust Multichip Average method using the Affy package of R (www.bioconductor.org/packages/release/bioc/html/affy.html). Differential analysis was performed for CIMP+ vs. CIMP− using the limma package (16) of R. |Log (fold change)| >1.0 and P<0.05 were set as cut-offs for significant differential expression.

Functional enrichment analysis

The Gene Ontology (GO; www.geneontology.org/) database is a bioinformatics resource that can provide functional categorization and annotations for gene products via the use of structured, controlled vocabularies (17). The Kyoto Encyclopaedia of Genes and Genome (KEGG; www.genome.jp/kegg) is a database for systematic analysis of the functions of genes or proteins in several specific metabolic and regulatory pathways (18). Functional enrichment analyses of the GO and KEGG databases were conducted using the Database for Annotation, Visualization and Integration Discovery (david.abcc.ncifcrf.gov/) (19). The statistical method for this was based on hypergeometric distribution. P<0.05 was considered to indicate significant functions and pathways.

Construction of protein-protein interaction (PPI) network

Proteins work together to complete certain biological functions. Therefore, revealing PPI is useful in elucidating underlying molecular mechanisms. In the present study, PPI networks were constructed for upregulated and downregulated genes using information from STRING (20). Interactions with the required level of confidence (i.e., score >0.4) were retained in the network. The two networks were visualised using Cytoscape (21).

Proteins in the network were presented as ‘nodes’, and each pairwise protein interaction was represented by an undirected link and the ‘degree’ of a node corresponded to the number of interactions by the protein. ‘Degree’ was calculated for each node.

Prediction of miRNAs and construction of the whole regulatory network

Web-based Gene Set Enrichment Analysis Toolkit (WebGestalt; www.webgestalt.org/option.php) is a comprehensive and powerful analysis toolkit, which can be used for enrichment analysis and microRNA (miRNA)-target prediction by identifying miRNA-binding site motifs. In the present study, miRNAs regulating DEGs were predicted using WebGestalt (22). Count ≥2 was set as the cut-off for predicted miRNAs and the top 10 miRNAs were selected. Following miRNA-target gene network pairs using WebGestalt, PPI networks and miRNA-target gene interactions were integrated. Subsequently, the whole regulatory network was visualised using Cytoscape (21).

Results

DEGs

A total of 41,335 genes were detected and 439 DEGs between CIMP+ and CIMP− samples were identified, including 241 upregulated and 198 downregulated genes in CIMP+ samples.

The GO biological pathway terms enriched for the 241 upregulated genes in CIMP+ samples could be divided into 13 clusters. They were associated with extracellular matrix organisation, defence response, immune response, collagen fibril organisation, and regulation of cell motion. The top 10 terms are listed in Table I.

Table I.

GO biological process terms enriched in the differentially expressed genes.

A, Upregulated genes
GO term	Count	(%)	P-value
GO:0030198 extracellular matrix organization	13	5.676855895	2.62×10−8
GO:0030199 collagen fibril organization	8	3.493449782	1.25×10−7
GO:0002504 antigen processing and presentation of peptide or polysaccharide antigen via major histocompatibility complex class II	8	3.493449782	3.25×10−7
GO:0009611 response to wounding	25	10.91703057	4.95×10−7
GO:0006955 immune response	28	12.22707424	1.58×10−6
GO:0043062 extracellular structure organization	13	5.676855895	3.55×10−6
GO:0006952 defense response	25	10.91703057	6.64×10−6
GO:0016064 immunoglobulin mediated immune response	8	3.493449782	1.05×10−5
GO:0019724 B cell mediated immunity	8	3.493449782	1.34×10−5
GO:0006954 inflammatory response	17	7.423580786	1.59×10−5
B, Downregulated genes
GO term	Count	(%)	P-value
GO:0007423 sensory organ development	11	5.882352941	2.00×10−4
GO:0007155 cell adhesion	20	10.69518717	2.07×10−4
GO:0022610 biological adhesion	20	10.69518717	2.10×10−4
GO:0030182 neuron differentiation	15	8.021390374	2.97×10−4
GO:0048666 neuron development	13	6.951871658	3.27×10−4
GO:0044057 regulation of system process	12	6.417112299	5.55×10−4
GO:0048592 eye morphogenesis	6	3.20855615	9.10×10−4
GO:0051966 regulation of synaptic transmission, glutamatergic	4	2.139037433	1.07×10−3
GO:0031175 neuron projection development	10	5.347593583	1.93×10−3
GO:0015672 monovalent inorganic cation transport	11	5.882352941	2.48×10−3

GO, Gene Ontology.

The GO biological pathway terms enriched for the 198 downregulated genes were divided into 12 clusters. They were associated with cell adhesion, sensory organ development, system process regulation, neuron differentiation and membrane organisation. The top 10 terms are listed in Table I.

KEGG pathway enrichment analysis revealed 16 significant pathways associated with upregulated genes (Table II), including focal adhesion (hsa04510), asthma (hsa05310), ECM-receptor interaction (hsa04512), intestinal immune network for immunoglobulin A production (hsa04672) and allograft rejection (hsa05330). No significant pathway was identified for the downregulated genes.

Table II.

Kyoto Encyclopaedia of Genes and Genome pathways enriched in the upregulated genes.

Term	Count	P-value
hsa04510:Focal adhesion	15	2.96×10−6
hsa05310:Asthma	7	5.16×10−6
hsa04512:Extracellular matrix-receptor interaction	10	6.52×10−6
hsa04672:Intestinal immune network for immunoglobulin A production	8	1.09×10−5
hsa05330:Allograft rejection	7	1.94×10−5
hsa05322:Systemic lupus erythematosus	10	2.52×10−5
hsa05332:Graft-vs.-host disease	7	3.12×10−5
hsa04514:Cell adhesion molecules	11	4.31×10−5
hsa04940:Type I diabetes mellitus	7	4.83×10−5
hsa05416:Viral myocarditis	8	1.27×10−4
hsa05320:Autoimmune thyroid disease	7	1.48×10−4

PPI networks of the DEGs

A PPI network containing 134 nodes and 314 edges was constructed for the upregulated genes (Fig. 1), whereas a PPI network consisting of 85 nodes and 80 edges was obtained for the downregulated genes (Fig. 2).

Figure 1.

PPI network for upregulated genes. The circular nodes represent the names of genes and lines represent the interactions between the genes. PPI, protein-protein interaction.

Figure 2.

PPI network for downregulated genes. The circular nodes represent the name of genes and lines represent the interactions between the genes. PPI, protein-protein interaction.

The top ten nodes with a high degree in the up and downregulated PPI networks are listed in Table III. The top five nodes in the network of upregulated genes were collagen type III α1 (COL3A1), collagen type V α2 (COL5A2), TIMP metallopeptidase inhibitor 1 (TIMP1), collagen type V α1 (COL5A1) and vimentin (VIM). In the network of downregulated genes, the top six nodes were glutamate receptor ionotropic AMPA2 (GRIA2), bone morphogenetic protein 2 (BMP2), protein kinase X-linked (PRKX), v-myc avian myelocytomatosis viral oncogene homolog (MYC), tight junction protein 2 (TJP2) and platelet-derived growth factor receptor α polypeptide (PDGFRA).

Table III.

Top 10 nodes with a high degree in the up and downregulated protein-protein interaction network.

Gene	Degree
Upregulated
COL3A1	22
COL5A2	18
TIMP1	16
COL5A1	16
VIM	15
ANXA2	13
S100A6	12
ANXA1	12
COL4A1	11
CXCL10	11
Downregulated
GRIA2	6
BMP2	5
PRKX	5
MYC	5
TJP2	5
PDGFRA	5
DCX	4
SH3GL2	4
RTN1	4
ID1	4

PPI, protein-protein interaction.

miRNA prediction and regulatory network analysis

The distribution of upregulated and downregulated genes in biological processes was analysed using WebGestalt (Figs. 3 and 4, respectively). The regulatory miRNAs of DEGs were also predicted (Table IV). Among these predicted miRNAs, miRNA-506 and miR-34b targeted the most DEGs in the up- and downregulated regulatory networks. In the upregulated regulatory network, miRNA-506 (miR-506) regulated five upregulated genes: VIM, aryl hydrocarbon receptor, proteolipid protein 2, IQ motif-containing GTPase activating protein 1 (IQGAP1) and syndecan 4. In the downregulated regulatory network, miR-34b regulated four downregulated genes: Sex determining region Y-box 4, PDGFRA, activated leukocyte cell adhesion molecule (ALCAM) and MYC. All predicted miRNAs with their target DEG pairs are presented in the regulatory network (Fig. 5).

Figure 3.

The distribution of upregulated genes in biological processes categories, as assessed by WebGestalt.

Figure 4.

The distribution of downregulated genes in biological processes categories, as assessed by WebGestalt.

Table IV.

Predicted miRs targeting the differentially expressed genes.

A, Upregulated genes
miR	Gene	Statistics
hsa_TTGCACT, miR-130a, miR-301, miR-130b	2	C=52; O=2; E=6.41; R=0.31; raw P=1.0000; adj P=1.0000
hsa_TTTGCAC, miR-19a, miR-19b	2	C=71; O=2; E=8.76; R=0.23; raw P=1.0000; adj P=1.0000
hsa_TGGTGCT, miR-29a, miR-29b, miR-29c	3	C=59; O=3; E=7.28; R=0.41; raw P=1.0000; adj P=1.0000
hsa_TGCCTTA, miR-124a	3	C=84; O=3; E=10.36; R=0.29; raw P=1.0000; adj P=1.0000
hsa_GTGCCTT, miR-506	5	C=105; O=5; E=12.95; R=0.39; raw P=1.0000; adj P=1.0000
hsa_ACATTCC, miR-1, miR-206	3	C=61; O=3; E=7.52; R=0.40; raw P=1.0000; adj P=1.0000
hsa_ACACTCC, miR-122a	2	C=13; O=2; E=1.60; R=1.25; raw P=0.4895; adj P=1.0000
hsa_ATGTTTC, miR-494	2	C=28; O=2; E=3.45; R=0.58; raw P=1.0000; adj P=1.0000
hsa_GGGACCA, miR-133a, miR-133b	2	C=37; O=2; E=4.56; R=0.44; raw P=1.0000; adj P=1.0000
hsa_ATTCTTT, miR-186	2	C=45; O=2; E=5.55; R=0.36; raw P=1.0000; adj P=1.0000
B, Downregulated genes
miR	Gene	Statistics
hsa_ACTGCCT, miR-34b	4	C=41; O=4; E=4.11; R=0.97; raw P=1.0000; adj P=1.0000
hsa_CACTGCC, miR-34a, miR-34c, miR-449	3	C=47; O=3; E=4.71; R=0.64; raw P=1.0000; adj P=1.0000
hsa_AAACCAC, miR-140	2	C=25; O=2; E=2.50; R=0.80; raw P=1.0000; adj P=1.0000
hsa_TAGCTTT, miR-9	2	C=31; O=2; E=3.11; R=0.64; raw P=1.0000; adj P=1.0000
hsa_TGCACTT, miR-519c, miR-519b, miR-519A	3	C=54; O=3; E=5.41; R=0.55; raw P=1.0000; adj P=1.0000
hsa_GTTAAAG, miR-302b	2	C=9; O=2; E=0.90; R=2.22; raw P=0.2255; adj P=1.0000
hsa_AACTGGA, miR-145	2	C=33; O=2; E=3.31; R=0.61; raw P=1.0000; adj P=1.0000
hsa_ACCAAAG, miR-9	2	C=67; O=2; E=6.71; R=0.30; raw P=1.0000; adj P=1.0000
hsa_TGCTGCT, miR-15a, miR-16, miR-15b, miR-195, miR-424, miR-497	2	C=86; O=2; E=8.61; R=0.23; raw P=1.0000; adj P=1.0000
hsa_CTGAGCC, miR-24	3	C=35; O=3; E=3.51; R=0.86; raw P=1.0000; adj P=1.0000

miR, microRNA; C, number of reference genes in the category; O, number of genes in the gene set and also in the category; E, expected number in the category; R, ratio of enrichment; raw P, P-value from hypergeometric test; adj P, adjusted P-value.

Figure 5.

The regulatory network of differentially expressed genes. Upregulated genes are in red, while downregulated ones are in green. Yellow triangles represent microRNAs.

Discussion

In the present study, a total of 439 DEGs were identified, including 241 upregulated and 198 downregulated genes. Functional enrichment analysis predicted that upregulated genes were associated with extracellular matrix organisation, defence response, immune response, collagen fibril organisation and regulation of cell motion, whereas downregulated genes were associated with cell adhesion, sensory organ development, regulation of system process, neuron differentiation and membrane organisation. These findings are consistent with previous reports (23–26). Ulrich et al (27) pointed out that the mechanical rigidity of the extracellular matrix regulates the structure, motility and proliferation of glioma cells. Cell motion and cell adhesion were closely associated with the invasion of glioma cells.

In the present study, a PPI network containing 134 nodes and 314 edges was constructed for upregulated genes, whereas a PPI network consisting of 85 nodes and 80 edges was also obtained for downregulated genes. The top five nodes in the network of upregulated genes were COL3A1, COL5A2, TIMP1, COL5A1 and VIM. TIMP1, as an inhibitor of matrix metalloproteinases, can promote cell proliferation and may have anti-apoptotic function (28). Groft et al (29) reported the differential expression and localisation of TIMP-1 and TIMP-4 in human gliomas and suggested that they may contribute to the pathophysiology of human malignant gliomas. In addition, Aaberg-Jessen et al (30) demonstrated that low expression of tissue inhibitor of TIMP-1 in glioblastoma predicts longer patient survival. Serum TIMP-1 level is also regarded as an independent predictor of survival (31). VIM is a member of the intermediate filament family and functions as an organiser of numerous critical proteins involved in cell adhesion, migration and cell signalling (32). Overexpression of VIM has been reported in central nervous system tumours and it strongly correlates with accelerated tumour growth, invasion and poor prognosis (33). The top six nodes in the network of downregulated genes in the present study were GRIA2, BMP2, PRKX, MYC, TJP2 and PDGFRA. BMP2 is involved in cell differentiation. Deregulation of the BMP developmental pathway in glioblastoma-initiating cells contributes to their tumorigenicity both by desensitising cells to normal differentiation cues and converting otherwise cytostatic signals to pro-proliferative ones (34). Liu et al (35) indicated that BMP2 expression levels may be a potent tool for assessing the clinical prognosis of glioma patients. In addition, Wang et al (36) reported that c-MYC is required for the maintenance of glioma cancer stem cells. Furthermore, Jensen et al (37) demonstrated that astroglial c-MYC overexpression predisposes mice to primary malignant gliomas. Overexpression of PDGFRA has also been reported in gliomas (38,39). Taken these findings together, TIMP-1, VIM, BMP2, MYC and PDGFRA may be associated with the development of glioma.

miRNAs regulating upregulated and downregulated genes identified in the present study were predicted using WebGestalt, including miR-124a and miR-34a. miR-124a can inhibit the proliferation of glioblastoma multiforme cells and induce the differentiation of brain tumour stem cells (40). It is frequently downregulated in glioblastoma and is involved in migration and invasion (41). IQGAP1 is one of its target genes, which has been implicated in the regulation of E-cadherin-mediated cell-cell adhesion (42,43). In addition, miR-34a can inhibit glioblastoma growth by targeting multiple oncogenes (44). In the present study, we predicted that miR-34a regulates PDGFRA and ALCAM, which were downregulated in CIMP+ samples. This is noteworthy because PDGFRA is involved in tumour progression (45) and its suppression by targeting miR-34a could contribute to tumorigenesis in pro-neural malignant gliomas (46), whereas ALCAM is associated with cell adhesion and cell migration (47). Therefore, it is speculated that these miRNAs may be useful for treating gliomas but further confirmation is needed.

Although we identified several DEGs that were important to define a distinct subgroup of glioma and understand the progression of glioma CIMP, there are certain limitations in the present study. The association between DEGs and methylation level in the different CIMPs was not investigated due to the lack of information on DEGs methylation levels in the dataset used. Additionally, experimental or data verification for the DEGs identified in glioma CIMP was not conducted, and in future, samples should be divided into different CIMPs based on methylation analysis to conduct the experimental validation.

In conclusion, several key genes were identified in glioma CIMP, some of which (TIMP1, VIM, BMP2, c-MYC and PDGFRA) may be viewed as potential markers or therapeutic targets for gliomas. In addition, relevant miRNAs, such as miR-124a and miR-34a that regulate genes involved in gliomas were also detected. These findings may provide helpful guidance to reveal molecular mechanisms underlying glioma CIMP.