Ependymomas (EPNs) are one of the most commonly diagnosed malignant neuroepithelial tumors in children and adults, accounting for ~1.8% of all primary central nervous system (CNS) tumors (1). It has been demonstrated that ependymomas in children <5 years of age account for >50% of all cases (2). Following comprehensive treatment with gross total resection and adjuvant radiotherapy, only 70% of patients with intracranial ependymoma <5 years of age were cured (3). In 2016, the latest World Health Organization classification for ependymoma was published as follows: Subependymoma; myxopapillary ependymoma; ependymoma; anaplastic ependymoma; and ependymoma, RELA proto-oncogene, NF-κB subunit fusion-positive (4). However, the pathological criterion is limited for clinical application; it was demonstrated that 7% of cases were misdiagnosed and were subsequently reclassified as ependymoma, which is indicative of the requirement for a more precise molecular classification system (5). Therefore, an increased number of studies are required to improve our understanding of the molecular mechanisms involved in EPN.microRNAs (miRNAs/miR) are a class of non-protein-coding small single-stranded RNAs, which bind to the 3′-untranslated region (UTR) of target mRNAs, inhibiting gene expression at the post-transcriptional level (6,7). It was suggested that high expression levels of miR-124-3p significantly decreased the progression-free survival time of patients with EPN by negatively regulating tumor protein p53 nuclear protein 1 (TP53INP1), thereby indicating a potential role of miR-124-3p as a therapeutic biomarker (8). Other studies involving microarray analyses revealed that the expression level of cyclin D1, which is involved in DNA repair, was increased in recurrent EPNs compared with the primary EPNs (9). With considerable advances in high-throughput technologies, increasing numbers of miRNAs and genes have been confirmed to serve important roles in the diagnosis, treatment and prognosis of EPNs, which may lead to further investigations into the molecular mechanisms of its pathogenesis.In the present study, five gene expression profiles datasets (GSE25604, GSE50161, GSE66354, GSE74195 and GSE86574) and one miRNA expression profile dataset (GSE42657) were downloaded from the Gene Expression Omnibus (GEO) to identify differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs) between EPN samples and normal brain tissue samples. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, in addition to protein-protein interaction (PPI) network analyses of these DEGs, were performed to identify associated hub genes and modules. Furthermore, an miRNA-mRNA network analysis was conducted to identify the key miRNAs (and their targeted genes) that may be utilized as biomarkers of EPN. Taken together, the present study may provide an improved understanding of the genetic mechanisms of EPN.
Materials and methods
Microarray data
mRNA (GSE25604, GSE50161, GSE66354, GSE74195 and GSE86574) and miRNA (GSE42657) expression profiles were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/); data from the GSE25604 dataset, including 15 EPN and 7 normal samples, generated from GPL571 platform (Affymetrix Human Genome U133A 2.0 Array) were submitted on Nov 24, 2010 and updated on April 05, 2017. The GSE50161 (including 46 EPN and 13 normal samples), GSE66354 (including 64 EPN and 13 normal samples), GSE74195 (including 13 EPN and 5 normal samples) and GSE86574 (including 29 EPN and 10 normal samples) datasets were generated using the same microarray platform. miRNA expression data from the GSE42657 dataset were submitted on November 30, 2012 and updated on February 12, 2016. GSE42657, from the GPL8179 platform (illumina Human v2 MicroRNA expression beadchip) contained 14 EPN and 7 normal samples. Survival data of EPN from GEO and The Cancer Genome Atlas databases were not available.Furthermore, 10 EPN and 10 normal brain tissue samples were obtained from the Department of Neurosurgery, Shanghai Ninth People's Hospital Affiliated to Shanghai Jiao Tong University School of Medicine (Shanghai, China). The present study was approved by the Ethics Committees of the Shanghai Ninth People's Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, and informed patient consent was obtained prior to clinical research. All tissue samples were immediately frozen in liquid nitrogen and subsequently stored at −80°C.
Identification of DEGs and DEMs
Following transformation into gene symbols according to probe annotation files of each platform, the original expression data were preprocessed using Linear models for microarray data using Limma package in R (v.3.4.1; http://www.r-project.org/). This included background correction, quantile normalization, summarization and probe ID-to-gene symbol transformation.GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) is an online tool used to identify differentially expressed molecules among ≥ two sample groups in a GEO dataset. In the present study, GEO2R was used to identify DEGs and DEMs between EPN and normal samples, and a threshold of adjusted P<0.05 and |log2 fold-change (FC)|>1 was applied. To control for Type I errors and false discovery rate, P-values were adjusted using the Benjamini & Hochberg false discovery rate method provided by Limma package (10).
Functional and pathway enrichment analysis
Following data preprocessing and DEG acquisition, aberrantly expressed DEGs were uploaded to the Database for Annotation, Visualization and Integrated Discovery (DAVID; v.6.8; http://david.ncifcrf.gov/) for functional and pathway enrichment analysis (11). KEGG pathway database records networks of cellular molecular interactions, and variants of these interactions specific to particular organisms (12). GO enrichment analysis includes ‘Biological process’ (BP), ‘Molecular function’ (MF) and ‘Cellular component’ (CC). In addition, the criterion P<0.05 was considered to indicate statistically significance. Reactome pathway analysis was also conducted (https://www.reactome.org/) (13), the results of which were visualized using Cytoscape 3.5.1 (14).
PPI network construction and module selection
Following uploading of the DEGs to the Search Tool for Retrieval of Interacting Genes (STRING) database (version 10.5; http://www.string-db.org/), a PPI network was constructed and visualized in Cytoscape (v.3.5.1) (14,15). Hub genes were identified using cytoHubba plug-in and DEG modules were distinguished using Molecular Complex Detection (MCODE) plug-in (16,17). The criteria were set as scores >3 and nodes >4. Furthermore, in order to conduct GO Biological process (BP) and KEGG pathway enrichment analysis, DEGs in the selected modules were further analyzed using DAVID. Terms with P<0.05 were considered to be significant. The transcription factor (TF)-gene regulatory network was performed using the iRegulion plug-in (18).
Identification of potential miRNA target genes
miRWalk2.0 (http://zmf.umm.uni-heidelberg.de/mirwalk2) comprises a collection of predicted and validated miRNA-target interactions. The predicted interactions were assessed using the following 12 miRNA target prediction programs: DIANA-microT v4.0 (19), DIANA-microT-CDS (20), miRanda-rel2010 (21), mirBridge (22), miRDB 4.0 (23), miRmap (24), miRNAMap (25), PicTar 2 (26), PITA (27), RNA22 v2 (28), RNAhybrid 2.1 (29) and Targetscan 6.2 (30). In order to improve prediction reliability, the genes predicted by ≥8 of the databases among 12 databases were selected as the DEM targets. Validated databases were then used to search for verified targeted genes. Finally, the sums of predicted and validated genes, which were repeated with the obtained DEGs, were recognized as the eventual target genes of the identified DEMs.
Construction of a negative miRNA-mRNA regulatory network
miRNAs promote mRNA degradation or inhibit the translation of negatively regulated targeted genes (7). Therefore, a negative miRNA-mRNA regulatory network was selected and visualized in Cytoscape. The top 3 significant DEMs and their target hub genes were also identified.
Association analyses and efficacy evaluation
To verify the clinical value of the hub genes, association analyses and efficacy evaluation were conducting using mRNA (GSE25604, GSE50161, GSE66354, GSE74195 and GSE86574) and miRNA (GSE42657) expression profiles, respectively. The associations between hub genes and the corresponding clinical features were illustrated using by boxplots, and statistical significance was determined using an independent sample t-test. Efficacy evaluation was conducted by receiver operating characteristic (ROC) curve analysis with the ‘pROC’ package in R language. Genes with an area under the curve (AUC) value >0.7 were considered to distinguish tumor tissues from normal tissues.
Reverse transcription-quantitative (RT-q) PCR
Total RNA was extracted from tissues using TRIzol® reagent (Takara Biotechnology Co., Ltd.,) according to the manufacturer's protocol. RT-qPCR was performed using the FastStart essential DNA Green Master mix (Roche Diagnostics) and a miRNA qPCR Assay kit (CWBio Technology Co., Ltd.), and performed using the iCycler iQ Real-Time Detection System (Bio-Rad Laboratories, Inc.). snRNA U6 and GAPDH were used as internal controls for miR and mRNA, respectively. The expression levels were quantified using the 2−ΔΔCq method (31) and the fold change for target genes was normalized to the appropriate internal control. The primer sequences are listed in Table SI.
Statistical analysis
SPSS software (v.22.0; IBM Corp.), GraphPad Prism (v.7.0; GraphPad Software, Inc.) and R software (v.3.4.1) were used for statistical analysis. The statistical significance between EPN and normal brain tissues was determined using an independent sample t-test. P<0.05 was considered to indicate a statistically significant difference. The bar plots were generated by GraphPad Prism or R software.
Results
Identification of DEGs
Following preprocessing, boxplots of the samples from five selected datasets are demonstrated in Fig. S1. Based on the cut-off criteria of P<0.05 and |log2 FC|>1, a total of 12,115, 6,030, 4,493, 3,026 and 7,751 DEGs between EPN and normal brain samples were identified from the GSE25604, GSE50161, GSE66354, GSE74195 and GSE86574 datasets, respectively (Fig. S2). Among these datasets, a total of 948 overlapping DEGs were identified (Fig. 1), including 609 upregulated and 339 downregulated genes, were screened out among the five datasets. In addition, the expression levels of 129 DEMs were identified as being significantly different between EPN and normal brain tissues, of which 39 were upregulated and 80 were downregulated miRNAs in EPN tissues (Table SII).
Figure 1.
Venn diagram of differentially expressed genes among the five datasets.
Following uploading of the 609 upregulated and 339 downregulated DEGs to DAVID, GO categories and KEGG pathways with P<0.05 were identified, and the top 5 significant terms were selected from corresponding GO categories on up- and downregulated DEGs (Table SIII). GO analysis results revealed that upregulated DEGs were significantly associated with ‘extracellular exosome’ and ‘extracellular matrix’ in CC, ‘extracellular matrix organization’ and ‘cell adhesion’ in BP, and ‘extracellular matrix structural constituent’ and ‘protein binding’ in MF. The downregulated DEGs were significantly enriched in ‘cell junction’ and ‘plasma membrane’ in CC, ‘chemical synaptic transmission’ and ‘neurotransmitter secretion’ in BP, and ‘calcium ion binding’ and ‘GABA-A receptor activity’ in MF.In addition, the top 10 significantly enriched pathways of the up- and downregulated DEGs were selected from significant KEGG pathways (Table SIV). The upregulated DEGs were significantly enriched in ‘ECM-receptor interaction’, ‘focal adhesion’ and ‘graft-versus-host disease’, while downregulated DEGs were associated with ‘GABAergic synapse’, ‘retrograde endocannabinoid signaling’ and ‘morphine addiction’. Reactome pathway analysis results indicated that DEGs were associated with ‘G1/S-specific transcription’ and ‘NOTCH1 intracellular domain regulates transcription’ (Fig. 2).
Figure 2.
Significantly enriched Reactome pathway analysis of differentially expressed genes.
A PPI network of 948 DEGs was produced using the STRING online tool, with a cutoff of score >0.7 (Fig. 3; Table SV). Using the cytoHubba plug-in, DEGs with the top 20 betweenness scores and with degrees ≥24 (top 21) were identified as hub genes. A total of 6 nodes were selected as hub genes, including cyclin dependent kinase 1 (CDK1), CD44 molecule (Indian blood group) (CD44), proliferating cell nuclear antigen (PCNA), MYC, synaptotagmin 1 (SYT1) and kinesin family member 4A (KIF4A) (Table SVI). Among these 6 genes, betweenness of MYC was the greatest, and at 48, CDK1 had the highest number of node degrees. The TF-gene regulatory network is demonstrated in Fig. S3.
Figure 3.
Protein-protein interaction regulatory network of DEGs in ependymoma. Nodes represent DEGs. Red nodes indicate upregulated DEGs and green nodes indicated downregulated DEGs. Nodes with higher degree values are depicted with larger shapes. Edges/lines stand for the regulatory association between any 2 nodes. DEGs, differentially expressed genes.
Furthermore, MCODE analysis revealed a total of 19 available modules, from which the top 5 significant modules were selected (Fig. 4); BP and KEGG enrichment analysis of DEGs between these modules were subsequently performed (Table SVII). Enrichment analysis demonstrated that the DEGs in module 1 were primarily enriched in ‘morphine addiction’ and ‘GABAergic synapse’ (KEGG analysis), and ‘chemical synaptic transmission’ in the BP category of the GO analysis. In the KEGG pathway analysis, the DEGs in modules 2–5 were mostly associated with ‘protein digestion and absorption’, ‘pathogenic Escherichia coli infection’, ‘spliceosome’ and ‘GABAergic synapse’, respectively.
Figure 4.
Top 5 modules in the protein-protein interaction network for DEGs. (A) Module 1. (B) Module 2. (C) Module 3. (D) Module 4. (E) Module 5. Nodes represent DEGs. Edges/lines stand for the regulation association between any 2 nodes. Red and green nodes represent upregulated and downregulated genes, respectively. DEGs, differentially expressed genes.
Construction of a negative miRNA-gene regulatory network
Among the 39 upregulated DEMs and 339 downregulated DEGs, 19,319 predicted and 16,116 validated miRNA-target pairs were identified. A total of 26,033 predicted and 15,986 validated miRNA-target pairs were also identified from 80 downregulated DEMs and 609 upregulated DEGs, respectively. Subsequently, 812 miRNA-DEG pairs from upregulated DEMs and downregulated DEGs were selected for the construction of a miRNA-DEG network, which was visualized using Cytoscape (Fig. 5). Similarly, a network including 1,232 miRNA-DEG pairs from downregulated DEMs and upregulated DEGs was also generated (Fig. 6).
Figure 5.
Regulatory network of upregulated miRNAs and downregulated DEGs in ependymoma. The red triangles represent the upregulated miRNAs. The green circles indicate downregulated DEGs. miRNA, microRNA; DEGs, differentially expressed genes.
Figure 6.
Regulatory network between downregulated miRNAs and upregulated DEGs in ependymoma. The green triangles represent downregulated miRNAs. The red circles indicate upregulated DEGs. miRNA, microRNA; DEGs, differentially expressed genes.
Following analysis, upregulated homo sapiens (hsa)-miR-34a-5p, hsa-miR-449a and hsa-miR-106a-5p, were revealed to be the 3 most significantly upregulated miRs, whilst hsa-miR-124-3p, hsa-miR-128-3p and hsa-miR-330-3p were identified as the top 3 downregulated miR. Finally, 4 miRNA-DEG pairs (hsa-miR-449a-SYT1, hsa-miR-34a-5p-SYT1, hsa-miR-330-3p-CD44 and hsa-miR-124-3p-PCNA) were identified (Table I).
Table I.
miRNA-target pairs of top 3 miRNAs and hub genes from the miRNA-mRNA regulatory network.
Significant differences were observed in the expression levels of 3 mRNAs (SYT1, CD44 and PCNA) and 4 miRNAs (hsa-miR-449a, hsa-miR-34a-5p, hsa-miR-330-3p and hsa-miR-124-3p) between EPN and normal tissues (P<0.05; (Fig. 7). The independent sample t-test and AUC analyses confirmed that each of mRNA and miRNA (including hsa-miR-449a, hsa-miR-34a-5p, hsa-miR-330-3p and hsa-miR-124-3p) possessed high specificity and sensitivity values, indicating their potentials as prognostic biomarkers (Fig. 8).
Figure 7.
Boxplots of association analyses. Boxplots of (A) SYT1, (B) CD44 and (C) PCNA in the GSE25604, GSE50161, GSE66354, GSE74195 and GSE86574 datasets. (D) Boxplots of hsa-miR-449a, hsa-miR-34a-5p, hsa-miR-330-3p and hsa-miR-124-3p from the GSE42657. *P<0.05 and **P<0.01. SYT1, synaptotagmin 1; CD44, CD44 molecule (Indian blood group); PCNA, proliferating cell nuclear antigen; hsa, homo sapiens; miR, microRNA.
Figure 8.
ROC analyses for efficacy evaluation. ROC analyses of (A) SYT1, (B) CD44 and (C) PCNA in the GSE25604, GSE50161, GSE66354, GSE74195 and GSE86574 datasets. (D) ROC analyses of hsa-miR-449a, hsa-miR-34a-5p, hsa-miR-330-3p and hsa-miR-124-3p from the GSE42657. AUC, area under the curve; hsa-miR, homo sapiens microRNA; ROC, receiver operating characteristic; SYT1, synaptotagmin 1; CD44, CD44 molecule (Indian blood group); PCNA, proliferating cell nuclear antigen; hsa, homo sapiens; miR, microRNA.
RT-qPCR analysis
Compared with the normal tissue samples, the expression levels of CD44, PCNA, hsa-miR-124-3p and hsa-miR-330-3p were significantly upregulated (P<0.01), and those of SYT1, hsa-miR-34a-5p and hsa-miR-449a were significantly downregulated (P<0.01) in EPN tissues. This was consistent with the results of the integrated bioinformatics analyses (Fig. 9).
Figure 9.
Boxplots of hub genes and miRs analyzed using reverse transcription-quantitative polymerase chain reaction. Significance between normal (n=10) and EPN tissues (n=10) for genes and miRs, including (A) CD44, (B) PCNA, (C) SYT1, (D) hsa-miR-124-3p, (E) hsa-miR-330-3p, (F) hsa-miR-34a-5p and (G) hsa-miR-449a, was determined using an independent sample t-test. **P<0.01. miRs, microRNAs; hsa, homo sapiens; EPN, ependymoma; CD44, CD44 molecule (Indian blood group); PCNA, proliferating cell nuclear antigen; SYT1, synaptotagmin 1.
Discussion
In previous years, microarray and next-generation sequencing technologies have developed rapidly, promoting more in-depth investigation into the molecular mechanisms of EPN. Evidence to indicate that miRNAs negatively regulate mRNA expression, and subsequently affect the development and progression of tumors, has been described (6,7). In the present study, 948 DEGs and 129 DEMs were selected from the GEO in an attempt to identify the molecular mechanisms and potential biomarkers of EPN.GO and KEGG analysis results provided the functional and pathway information for 948 DEGs and the top 5 modules identified. These genes were primarily associated with ‘cell division’, ‘mitotic nuclear division’ and ‘G2/M transition of mitotic cell cycle’, which implied a regulatory function in mitosis. Additionally, ‘retrograde endocannabinoid signaling’, ‘morphine addiction’ and ‘nicotine addiction’, identified using KEGG pathway analysis, indicated that EPN may be significantly associated with cannabinoids, morphine and nicotine. The ‘response to drug’ term, obtained from the GO analysis, was associated with 28 upregulated DEGs, which indicated a potential drug resistance mechanism and provided potential targets for patients with EPN exhibiting chemotherapeutic tolerance. The PI3K-Akt signaling pathway was has been suggested to be important in glioblastoma and medulloblastoma, which supports the results of the present study (32,33).Furthermore, following the integrated analysis, 6 hub genes (including CDK1, CD44, PCNA, MYC, SYT1 and KIF4A) were identified from the PPI network of DEGs. PCNA is an immunohistochemical marker of cellular proliferation in tumors, and is therefore able to predict the survival outcome of patients with EPN (34). The myc oncogenes comprise 3 principal genes: C-myc, N-myc and L-myc. Fluorescence in situ hybridization and immunohistochemistry data have indicated that C-myc, a targeted gene of the Notch pathway, was significantly correlated with the development of adult onset EPNs (35). Also, increased C-myc expression has been associated with poor prognosis in patients with low-grade tumors (36). Although not previously confirmed, an increasing number of studies have demonstrated the importance of the 4 hub genes (CDK1, CD44, SYT1 and KIF4A) in diseases, particularly CNS tumors (37–41). CDK1, a member of the cyclin-dependent kinase (CDK) family, serves an important role in G2-M phase transition, which is consistent with the GO analysis results of the present study (37). Following cDNA array and immunohistochemical analyses, a significant positive correlation was revealed between the expression level of CDK1 and glioma oncogenesis (38). CD44 was also suggested to be a confirmed biomarker for distinguishing the molecular subtype of glioblastoma multiforme (39). In conclusion, the present study identified 6 hub genes that serve essential roles in EPN, and which may function as notable diagnostic and therapeutic biomarkers.miRNAs are able to bind to the 3′-UTRs of specific genes to inhibit translation or promote the degradation of the corresponding mRNAs (7). Increasing evidence has indicated that miRNA dysregulation is responsible for the pathogenesis of EPN. In the present study, 6 miRNAs and 4 miRNA-DEG pairs were identified. Previous data has demonstrated that the expression level of miR-34a-5p was downregulated in colorectal cancer, which repressed apoptosis and cell cycle arrest at the G1 phase, and promoted p53 transcription to suppress tumor recurrence (42). CD117 was also identified as a direct target of miR-34-5p in the progression of osteosarcoma (43). Additionally, the overexpression of miR-449a in glioblastoma inhibited myc-associated zinc-finger protein activity through the PI3K/AKT pathway, which highlighted a novel miRNA biomarker (44). Zhi et al (45,46) revealed that miR-106a-5p served as a tumor suppressor in astrocytoma by inhibiting Fas-activated serine/threonine kinase, and that its expression was negatively correlated with clinical outcome. miR-124-3p was identified as a potential therapeutic and prognostic biomarker of EPN, which is consistent with the results of the present study, and further inhibited its target, TP53INP1, affecting clinical outcome (8). A study has indicated that miR-128-3p overexpression promoted neuronal survival in ischemia-induced brain injury (47). Additionally, miR-128-3p was demonstrated to be a suppressive biomarker in various malignancies, including lung cancer, acute lymphoblastic leukemia and hepatocellular carcinoma (48–50). miR-330-3p knockdown was also discovered to inhibit tumor growth, in contrast to the effect produced by its overexpression (51). Notably, Pantaleo et al (52) identified 3 miRNA-mRNA regulatory networks, including the miR-330-3p-CD44 pair, as biomarkers for the treatment of tyrosine protein kinase KIT/platelet derived growth factor receptor α wild type succinatedehydrogenase-deficient gastrointestinal stromal tumors (GISTs), which supported the results of the present study. However, the complex regulatory mechanisms of miRNAs in EPN remains to be fully elucidated.In conclusion, using integrated bioinformatics analysis, the present study revealed 948 DEGs and 129 DEMs associated with EPN. The results indicated 6 hub genes, 5 main modules, 6 crucial miRNAs and 4 miRNA-DEG pairs as novel potential biomarkers, which may facilitate further understanding of the molecular mechanisms of EPN. However, the functional value of these conclusions in EPN requires additional study.
Authors: Azra Krek; Dominic Grün; Matthew N Poy; Rachel Wolf; Lauren Rosenberg; Eric J Epstein; Philip MacMenamin; Isabelle da Piedade; Kristin C Gunsalus; Markus Stoffel; Nikolaus Rajewsky Journal: Nat Genet Date: 2005-04-03 Impact factor: 38.330
Authors: Hui Meng; Kai Wang; Xuedan Chen; Xingying Guan; Liwen Hu; Gang Xiong; Juan Li; Yun Bai Journal: Am J Cancer Res Date: 2015-02-15 Impact factor: 6.166
Authors: Yulia Margolin-Miller; Natalia Yanichkin; Keren Shichrur; Helen Toledano; Anat Ohali; Theophilos Tzaridis; Shalom Michowitz; Suzana Fichman-Horn; Meora Feinmesser; Stefan M Pfister; Hendrik Witt; Uri Tabori; Eric Bouffet; Vijay Ramaswamy; Cynthia Hawkins; Michael D Taylor; Isaac Yaniv; Smadar Avigad Journal: Genes Chromosomes Cancer Date: 2017-05-19 Impact factor: 5.006
Figure 1.
Figure 2.
Figure 3.
Figure 4.
Figure 5.
Figure 6.
Figure 7.
Figure 8.
Figure 9.