Differential and Common Signatures of miRNA Expression and Methylation in Childhood Central Nervous System Malignancies: An Experimental and Computational Approach

Simple Summary Epigenetic mechanisms, that are modifications of the genome without the presence of mutations, are known to play a crucial role in central nervous system (CNS) tumors during childhood. Two well-known epigenetic regulatory mechanisms include methylation and miRNA regulatory mechanisms. Therefore, in the present study we have investigated the presence of methylated genes in childhood CNS tumors, along with miRNA expression. We have searched for correlations between gene methylation and miRNA expression. In addition, we have investigated mRNA expression in order to search for possible miRNA targets. Such approaches could prove useful for the improvement of CNS tumor prognosis, as well as for the discovery of new therapeutic targets. Abstract Epigenetic modifications are considered of utmost significance for tumor ontogenesis and progression. Especially, it has been found that miRNA expression, as well as DNA methylation plays a significant role in central nervous system tumors during childhood. A total of 49 resected brain tumors from children were used for further analysis. DNA methylation was identified with methylation-specific MLPA and, in particular, for the tumor suppressor genes CASP8, RASSF1, MGMT, MSH6, GATA5, ATM1, TP53, and CADM1. miRNAs were identified with microarray screening, as well as selected samples, were tested for their mRNA expression levels. CASP8, RASSF1 were the most frequently methylated genes in all tumor samples. Simultaneous methylation of genes manifested significant results with respect to tumor staging, tumor type, and the differentiation of tumor and control samples. There was no significant dependence observed with the methylation of one gene promoter, rather with the simultaneous presence of all detected methylated genes’ promoters. miRNA expression was found to be correlated to gene methylation. Epigenetic regulation appears to be of major importance in tumor progression and pathophysiology, making it an imperative field of study.


Introduction
Cancer is one of the leading fatal diseases in the western world. Tumors of the central nervous system (CNS), are considered to be a complex, heterogeneous disease that is often

CNS Tumor Biomarkers
An important diagnostic tool for the determination and staging of CNS tumors is still the investigation of protein biomarkers, which are detected by immunohistochemistry (IHC) [3][4][5]. Numerous protein markers are used for the determination of tumor staging, where some include glial fibrillary acidic protein (GFAP), P38, neuromicrofilaments (NF), tumor protein p53 (p53), Bcl2 apoptosis regulator (BCL2), β-tubulin III (TUBB3), RNA binding Fox-1 homolog 3 (Neu-N), β-catenin (CTNNB1), pre-mRNA-splicing factor ini1 (INI1), marker of proliferation Ki-67 (MKI67), keratin, mucin 1 cell surface associated (MUC1 or EMA), epidermal growth factor receptor (EGFR), P27, S100 calcium binding protein (S100), actin, desmin (DES), myelin (MBP), and others. Immunohistochemistry consists of an invaluable tool for the differential diagnostic process and the accurate determination of tumors of the CNS [6]. Out of the numerous immunohistochemical biomarkers available, several have been highlighted as "independent prognostic markers" [7,8]. For example, the ki-67 is an antigen whose expression is related to the tumors' aggressiveness in the malignancies of the CNS [7,8]. The biological properties of ki-67 made it a significant biomarker in terms of tumor prognosis. Currently, the most prevalent therapeutic approach concerns the craniospinal irradiation, which is administered for both the local and metastatic disease in children older than three years of age [9]. Yet, disease prognosis is still poor and CNS tumors of the childhood are, unfortunately, still manifesting high mortality rates and therapy resistance. Thus it is possible that improved treatments will come from the in-depth understanding of the tumor's molecular machinery.

Gene Expression in CNS Tumors
There is overwhelming evidence on the role of gene expression and regulation in CNS tumors, both those of adults, as well of the children. miRNAs and mRNAs, have been studied thoroughly in the literature, where each has an abundant portion of literature dedicated to CNS tumors. A search in the principal literature databases returns more the 6000 articles on the topic of both miRNAs and mRNAs in CNS tumors.

Epigenetic Mechanisms in CNS Tumors
A very significant factor that has to be considered is the basic understanding of tumor biological mechanisms. In that sense, a very significant part of tumor mechanics is the epigenetic modification of the genome. Epigenetic mechanisms include the regulation of gene expression through methylation of genes or post-transcriptional modifications. Epigenetic mechanisms of biological systems has been neglected throughout time. Yet, it has been shown that they consist of a very significant regulatory mechanism.

DNA Methylation and Cancer
The role of DNA methylation in cancer has been the topic of intensive study during the recent years. Since the discovery of DNA methylation mechanisms, a large part of the literature has reported that DNA methyltransferase aberrant activity is present in tumors [10][11][12][13][14][15]. It has been found that malignant cells often manifest increased total DNA methyltrasferase activity, significantly extensive loss of methylation, from otherwise physiologically methylated promoters, but also hypermethylation of normally unmethlylated DNA sites [16][17][18][19][20]. Several studies have reported the significance of gene promoter methylation. For example, the genes SPARC, UCHL1, NPTX2, PENK, and PDAC were investigated in pancreatic cancer, where they were found to be hypermethylated in patients with pancreatic cancer [21]. Similarly, studies have reported the presence of hypermethylated tumor suppressor gene's promoters in lung cancer [22,23], gastric cancer [24], and breast cancer [25].

DNA Methylation and CNS Tumors
Recent findings have highlighted the role of methylation on CNS tumors. In particular, several reports have shown that mutations of the isocitrate dehydrogenase 1 (IDH1) and IDH2, and H3 histone family member 3A, are strongly associated with DNA and histone methylation, with a frequent methylation aberration being the O6-methylguanine-DNA methyltransferase promoter in human diffuse gliomas [26]. Similarly, the promoter of MGMT has been shown to be hypermethylated in non-malignant tumors of the CNS [27]. In addition, a recent report showed that in glioblastomas the promoter of GBX2, PDGFRA, and GLI2 were hypermethylated and also linked to poor prognosis and resistance to chemotherapy [28]. In a recent study, the methylation of RASSF1 and CASP8 were reported as significant markers for the separation between childhood ependymoma and choroid plexus papilloma [29].

Patient Administration and Stratification Based on New Biomarkers
One very important and interesting aspect that arises from studies as the one presented here, is the potential use of molecular factors that could be used for patient prognosis, diagnosis, and most of all therapy. A quick search in the literature up-to-date, shows that there is a great interest in the scientific community for the role of miRNAs in CNS tumors. For example, in a very recent study it has been reported that the signatures of a miRNA set, was able to distinguish between primary CNS and non-malignant tumors. The highlight of this study, was that the distinction could be facilitated by miRNAs detected in the cerebrospinal fluid (CSF) [30]. The great advantage of miRNAs (but limited to) is that besides their role in tumor tissue biology, they are also found in circulation, being able to be detected more easily. Another interesting issue, outlined in the literature (on which we also agree) is that clinical advantages (prognosis, diagnosis, therapy) are not derived from a single miRNA, but rather from a cohort or repertoire of miRNAs. This brings about another very important aspect in tumor biology, which is its multi-factorial and therefore complex nature [31]. On the other hand, the most important aspect that has being discussed recently, is the possibility of therapeutic interventions using epigenetic mechanisms, such as miRNAs. Currently, there are no known miRNA therapies for human brain tumors. Yet, several studies have highlighted the use of these molecules in animal models. In particular, miR-370-3p [32], miR-142-3p [33], miR-181 [34], miR-124, miR-128, and miR-137 [35] are molecules that, when used in animal models, manifested potential anti-tumor and chemotherapy-enhancement properties.
Another interesting aspect is the type of approach in the study of CNS biology and the role of epigenetic regulatory mechanisms that would be whether investigate the tumor on the tissue level or plasma level. Both approaches are significant and provide significant information on the tumor's biology. The investigation of miRNA levels in plasma/blood (thus circulating) could provide information for the prognostic and diagnostic course of the disease, while the investigation of epigenetic mechanisms on the tissue level provides information on the tumor's biology along with the diagnostic and therapy-related knowledge. A quick search on the bibliographic databases shows that most studies are concerned with biology of the disease on the tumor's site, while a smaller part is concerned with the circulating miRNA molecules. Yet, both approaches still remain quite significant for our understanding of the disease. From that perspective, we have chosen to investigate the biology of the tumors from the tissue approach, as we attempted to identify those miRNAs and epigenetic mechanisms that are common to the investigated tumors irrespective of their stratification and diversity.

Design and Aim of the Present Study
The current study, was performed on four levels; the first level included the diagnosis and staging of childhood CNS tumors through microscopy and immunohistochemical methods, the second level included the determination of miRNA expression levels using high throughput methodologies, the third level included the determination of mRNA levels using high throughput methodologies and the fourth level included the determination of the methylation status on specific genes. Our strategy consisted of two approaches; the first was to examine the relation of clinical variables and differential expressed genes and the second was to examine those miRNAs and mRNAs that were globally up-or down-regulated in all tumor samples. The present work is summarized in the flow chart of Figure S1.
Therefore, based on the aforementioned strategy, we aimed at determining the differential and common signatures of miRNA, mRNA expression and methylation in different childhood CNS tumors, in order to discover common regulatory mechanisms.

RNA Extraction
Samples (n = 61) were processed for both total RNA, as well as miRNA extraction. In brief, total RNA and miRNAs were extracted using the Trizol standard protocol (Invitrogen, Carlsbad, CA, USA) and the mirVANA miRNA isolation kit (Ambion, Austin, TX, USA). The RNA quantity and quality were evaluated using a spectrophotometer (NanoDrop ® ND-1000 UV-vis, Nanogen Inc., San Diego, CA, USA), as previously described [3][4][5].

miRNA Profiling
MicroRNA profiling has been described previously in detail [3][4][5], which we reproduce in quotes as follows: "In brief, total RNA and miRNAs were extracted using the Trizol standard protocol (Invitrogen, Carlsbad, CA, USA) and the mirVANA miRNA isolation kit (Ambion, Austin, TX, USA). Labelling and hybridization were performed using the LabelIT miRNA labelling kit (Mirus Bio LLC, Madison, WI, USA) according to manufacturer's instructions. Samples were hybridized to Applied MicroArrays (miRlink Bioarray 300054-3PK) platform. This array contained 1211 human miRNAs. Hybridization was performed at 37 • C with rotation at 145 rpm for 16 h. Images were scanned using Agilent Microarray Scanner (G2565CA) controlled by Agilent Scan Control 7.0 software. The total gene signals were extracted using the Imagene 6.0 software (Biodiscovery Inc., El Segundo, CA, USA) that contains summarized signal intensities for each miRNA by combining intensities of replicate probes and background subtraction" [3][4][5]. In total, 49 CNS tumor samples (the complete cohort described in Section 2.1.) and 13 control samples were investigated for their miRNA expressional profile. All microarray data are MIAME compliant.

cRNA Profiling
Oligos microarray chips (~57k genes) were obtained from GE HealthCare (Chicago, IL, USA) and Applied Microarrays (Tempe, AZ, USA) (formerly Amersham Biosciences, Buckinghamshire, UK) (CodeLink 57k Human Whole Genome) [38][39][40]. Hybridization was performed with the CodeLink RNA amplification and labeling kit as described by the manufacturer, utilizing the Cy5 fluorescent dye. Slides were scanned with a microarray scanner (ScanArray 4000XL, Northville, MI, USA). Images were generated with ScanArray microarray acquisition software (GSI Lumonics, Northville, MI, USA). cRNAs from three experimental setups were used in single experiments with internal spikes as controls. The scanned images were further processed with the CodeLink Expression Analysis Software v5.0 from Amersham Biosciences (presently GE Health Care Inc., Chicago, IL, USA). The experimental setup was analyzed based on the reference-design, as described previously [41][42][43]. Gene expression values of tumor samples were compared against the mean value of the control samples. In total, 4 CNS tumor samples (one PA, one MB, one EP, and one ATRT) and 3 control samples were investigated for their mRNA expressional profile. All microarray data are MIAME compliant. This type of experimentation was performed in order to be used as a reference population for the applied ontological annotations in the present study. In other words, each ontological annotation analysis, requires a reference population of genes to be used, which can be used the complete genome or a certain gene cohort. In our case, we have used the DE genes detected from mRNA expression analysis as the reference gene cohort and this was used for any further bioinformatics analyses.

Methylation
From the total patient cohort, children diagnosed with pilocytic astrocytomas (PA) (n = 13), ependymomas (EP) (n = 2), medulloblastoma (MB) (n = 8), and atypical teratoid rhabdoid tumor (ATRT) (n = 1) were evaluated for the possible methylation on specific genes. Twelve non-tumor brain samples were used (n = 12) as controls, as described in the previous section. Methylation-specific multiplex ligation-dependent probe amplification (MLPA) (MS-MLPA) is a semi-quantitative method for methylation profiling. In this study we used two MS-MLPA probe mixes (ME001-D1 tumor suppressor mix 1 and ME002-C1 tumor suppressor mix 2, MRC Holland (Amsterdam, NL) (www.mlpa.com, accessed 15 September 2021)) that contain methylation specific probes for 37 different tumor suppressor (TS) genes. Using the MS-MLPA technique we can distinguish the TS genes that keep the methylated status and therefore are silenced [44].

Methylation-Specific PCR
For the validation of methylated genes, a methylation-specific PCR (MS-PCR) protocol was designed [45]. In particular, samples underwent complete bisulfate conversion using the EZ DNA Methylation-LightningTM kit (Zymo Research, Irvine, CA, USA) before being used in the MS-PCR according to the manufacturer's instructions. The MS-PCR primers were designed using the MethPrimer software [46] and the PCR was performed at 55 • C. In particular, primers were designed for RASSF1 and CASP8 genes, in order to test more samples. Primer sets are summarized in Table 2.

Statistical Analysis
Continuous variables are expressed as mean ± standard deviation, unless indicated differently. Variable differences and association with clinical variables were conducted with the Kruskal-Wallis test. Patient's characteristics are presented with absolute and relative frequencies (%). Chi-square test of independence was used to evaluate the association between having multiple methylated genes or aberrations and patients' characteristics. The characteristics that were found statistically significant were entered in a logistic regression model in order to evaluate the probability of having multiple positive reactions. The level of statistical significance was set to a = 5%. For comparisons between groups, Student's t-test and one-way analysis of variance (ANOVA) were performed for the continuous variables and Chi-square tests were used for the categorical variables. Post hoc comparisons (adjusted with Bonferroni criterion) were also performed when significant differences (p < 0.05) of the variables in ANOVA tests were identified. The characteristics that were found statistically significant were entered in a logistic regression model in order to evaluate the probability of having multiple positive reactions. The modeling of a quantitative variable based on one or more qualitative and quantitative parameters, was performed through linear regression. Multiple logistic regression was performed in order to evaluate the probability of having multiple positive reactions. The relative risk (RR), odds ratio (OR), absolute risk (AR) were calculated. Methylation analysis and copy number variation (CNV) was performed with the Coffalyzer software (MRC Holland (www.mlpa.com, accessed 15 September 2021)). Indicative analysis results are presented in Figure S2. Statistical analysis has been performed with the Matlab ® computational environment (The Mathworks, Inc., The Natick, MA, USA).

Microarray Data Analysis Microarray Data Pre-Processing
The extracted microarray data were entered into a Microsoft Excel ® file.

Microarray Data Post-Processing
Microarray data were initially background-corrected using the multiplicative background correction (MBC) approach [47]. In particular, MBC subtracts the logarithmic estimates of the background intensity from the logarithmic foreground intensity. After background correction, negative values were removed and replaced with "NaN" values. Our intention was to find miRNA and mRNA expression, even those of low expression values. It is possible that not only those values that are of great difference are of importance, but also those that have very low expression values could be of biological importance. Microarray data post-processing has been performed with the Matlab ® computational environment (The Mathworks, Inc., The Natick, MA USA)

Microarray Data Normalization
Microarray data normalization was then performed using three algorithms; (a) Loess [48], (b) rank invariant, and (c) quantile algorithm. To account for differences across series we used the log2-transformed ratio: where, R i,j is the global mean-transformed ratio, x i,j is the expression value of gene i and sample j, X i,j is the restructured value of the ith gene and jth sample. The three algorithms were compared for their efficiencies. In general, the quantile algorithm performed better as comparted to the other two. Microarray data normalization has been performed with the Matlab ® computational environment (The Mathworks, Inc., The Natick, MA, USA).

Elimination of Duplicate Transcripts
To reduce the complexity of the dataset, we followed the replicate averaging approach proposed by Uzman et al. [49]. Calculations of duplicate genes was performed with the Matlab ® computational environment (The Mathworks, Inc., The Natick, MA, USA).

Detection of Differentially Expressed Genes (DEGs)
We used the Student's t-test [50] to identify the differentially expressed miRNA and mRNA genes (DE miRNAs, DE mRNAs) across all tumor samples and compared them to all control samples. The false discovery rate (FDR) was calculated, as previously described [51][52][53]. The DE miRNA and mRNA genes per experiment were identified at a confidence level of 95%. DE miRNAs and mRNAs were treated in two different ways. Data were further processed and analyzed as "ratios", i.e., as gene expression values calculated as the log2-transformed ratio of each tumor sample over the mean of all control samples, using the following formula: where E is the expression value, F Tumor,i,j is the expression value of tumor sample i and miRNA/mRNA j, F Controls,j is the mean expression value of all controls and miRNA/mRNA j. In addition, data were also analyzed as "naturals", meaning DEGs that are non-log2transformed and thus including the control and tumor samples separately. DEGs have been estimated with the Matlab ® computational environment (The Mathworks, Inc., The Natick, MA, USA).

Unsupervised Classification Methods
DEGs were further analyzed for common expression patterns using classification methods. To gain further insight into the gene expression data, we used unsupervised hierarchical clustering (HCL) and k-means classification [54,55]. HCL with dendrogram, was used and correlations were calculated with Euclidean distance. K-means classification [54,55] was recently reported as one of the best performing clustering approaches for microarray class discovery studies [56]. We applied the squared Euclidean as a distance measure, since it is generally considered to be a more appropriate measure for use with k-means and found to outperform for ratio-based measurements [57]. We used 100 iterations and the optimal cluster number for the k-means algorithm was estimated using the Calinski-Harabasz criterion. Complete k-means clusters, centroids, and sorted centroids [58] were utilized. The DE miRNAs/mRNAs, were also classified based on their diagnosis categorical variable. In particular, the mean values of all samples with respect to the diagnosis was estimated and the resulting descriptive statistical measure was utilized for further k-means classification. Gene expression was also analyzed with respect to the chromosomal distribution of the DE miRNAs/mRNAs. We explored the mean expression per chromosome and heat-maps of chromosomal-related expression. Unsupervised classification methods have been performed with the Matlab ® computational environment (The Mathworks, Inc., The Natick, MA, USA).

Unsupervised Classification Methods: Descriptive K-Means
In order to further investigate the expression patterns of DE miRNAs, we have used a variation of the k-means algorithm, where instead of using the complete sample cohort, i.e., each sample separately, we used the categorical properties of the sample cohort. In the classical k-means algorithm a matrix of m × n is used, where m is the number of genes/miRNAs and n is the number of samples. In our variation, a new matrix is formed with dimensions m × q, where q is the number of unique values of a categorical variable. For example, in the present cohort 49 tumor samples were utilized and could be divided into five distinct diagnostic categories (e.g., ependymoma, astrocytoma, medulloblastoma etc.) and, hence, q = 5. K-means clustering was implemented with the same parameterization as the aforementioned typical k-means approach.
Common Expression Patterns in DE miRNAs/mRNAs DE miRNAs/mRNAs were examined for possible common expression patterns, i.e., miRNAs/mRNAs that were either down-or up-regulated in all CNS tumor samples, irrespective of the tumor diagnosis. The clusters revealed by unsupervised classification were examined separately. Each miRNA was counted for its occurrences for up-or downregulation in all samples and the result was divided by the total number of samples, giving the percentage of up-or down-regulated samples of the respective miRNA/mRNA. We have looked for miRNAs/mRNAs that were either up or down-regulated in all samples (100%), in 90-99% of all samples, 80-89% of all samples, and 75-80% of all samples.

Gene Ontology (GO) Enrichment Analysis
We performed GO enrichment analysis using the gprofiler [59], and WebGestalt webtools [60]. Relations of the differentially expressed genes and the transcription factor binding motifs were further investigated using the Pubgene Ontology Database (www. pubgene.org, accessed 12 May 2020). Gene definitions and functions were based on the National Institute of Health databases (http://www.ncbi.nlm.nih.gov/sites/entrez/, accessed 27 September 2020).

Receiver Operating Characteristic (ROC) Analysis
ROC curves and naïve-Bayes classification were used to investigate the diagnostic ability of the co-deregulated miRNAs/mRNAs between CNS tumors and control samples. In the case of naïve-Bayes classification, the algorithm used the Bayes theorem, and (naïvely) assumes that the predictors are conditionally independent, given a class. Naïve Bayes classifiers assign observations to the most probable class (in other words, the maximum a posteriori decision rule). ROC analysis has been performed with the Matlab ® computational environment (The Mathworks, Inc., The Natick, MA, USA).

Ethics Statements
All experiments were conducted in compliance with the international biomedical studies stipulations, with reference to the Declaration of Helsinki of the World Medical Association. No personal data of patients were kept, while it was impossible to trace back any personal data from the data collected for the present study. Informed written consent was obtained from the parents of all children included in the study. The present study was conducted with the approval of "Aghia Sophia" Children's Hospital Ethics Committee (Protocol No. 35/19.16/09/13).

Results
In the present study, we have analyzed a cohort of pediatric patients suffering from different types of CNS tumors, including ependymoma, medulloblastoma, ATRT, and astrocytoma.

Patients and Tumor Samples
Our patient cohort included 39 males and 25 females with ages mean 6.59 ± 4.59 years, median 6.20 (min 0.03, max 16.06 years). The mean age of PA patients was 5.86 ± 3.42 years with a male-female ratio of 1:1, the median age for the EP patients was 5.13 ± 5.40 years with a male-female ratio of 6:1, the mean age for the MB patients was 6.92 ± 4.70 years with a male-female ratio of 1:1, the ATRT mean age was 2.46 ± 3.54, while the nonmalignant cohort was 9 years consisted of males only (Table 1). Among the patients' cohort, 24 patients (85.7%) remained in complete remission (after therapy); 17 PA, 5 EP patients, and 2 CD patients, and 4 patients (14.3%) succumbed from the disease; 2 PA patients and 2 EP patients. The laboratory and clinical parameters of all patients were examined in relation to the miRNA/mRNA profiles obtained. For part of the patient cohort, the information concerning the time of death was not known and, thus, it was not possible to calculate the overall survival for all patients.

miRNA and mRNA Expression
miRNA expression levels of all samples (n = 61) were estimated and were sorted in ascending order in order to identify expression patterns of DE miRNAs. In total, 75 miRNAs were detected as DEGs, with a FDR < 0.05 at p < 0.05 (Table S1). In addition, 869 mRNAs were detected as DEGs with a FDR < 0.01 at p < 0.05. DE mRNAs were used further on, in the present analysis, in order to identify miRNA targets with relation to the mRNA expression profile. Mean miRNA DEGs naturals in tumor and control samples separately are presented in Figure 1A, while miRNA ratios of tumors over controls are presented in Figure 1B Figure 2F). Interestingly, all six miRNAs were found to be up-regulated in female patients as compared to male patients. In the case of gender, we have estimated the ratios of gene expression and not the "natural" values due to the fact that gender "natural" values included both control and tumor values.

miRNA Expression Profiling with Respect to Gender
In the investigation of diagnosis, we have used both the "natural" values, as well as ratios of gene expression. Significant differences were observed between cortical dysplasia and ependymoma (p = 0.00066), as well as controls and ependymomas (p = 0.0055) with respect to miR-1234 ( Figure A). In addition, significant differences were observed between controls and ependymoma (p = 0.00003), as well as medulloblastoma and ependymomas (p = 0.0018) with respect to miR-183 ( Figure B). Similarly, significant differences were observed between controls and cortical dysplasias (p = 0.0044) with respect to miR-25* ( Figure C). Significant differences were observed between cortical dysplasia and pilocytic astrocytomas (p = 0.001), as well as cortical dysplasias and medulloblastomas (p = 0.002) with respect to miR-3675-5p ( Figure D), whereas miRNA expression manifested also a descending pattern from cortical dysplasias to ATRT ( Figure D). Further on, significant differences were observed between controls and cortical dysplasias (p = 0.0004) with respect to miR-612 ( Figure E) with a similar descending pattern as in the case of miR-3675-5p. By including control samples and calculating the expression ratios, significant differences were observed between medulloblastomas and ependymomas (p = 0.007) with respect to miR-122* (Figure F), between cortical dysplasia and ATRT (p = 0.004) with respect to miR-130b ( Figure G), ependymomas and ATRTs (p = 0.0009) with respect to miR-2909 ( Figure H). Interestingly, significant differences were observed between cortical dysplasias and ependymomas (p = 0.001) with respect to miR-302b ( Figure I), with an ascending pattern as we move from cortical dysplasias to ATRTs, indicating a possible relation to tumor aggressiveness. Subsequently, significant differences were observed between pilocytic astrocytomas and ependymomas (p = 0.004) with respect to miR-4251 ( Figure J), between medulloblastomas and ependymomas (p = 0.009) with respect to miR-576-5p ( Figure K), between medulloblastomas and ependymomas (p = 0.004) with respect to miR-600 ( Figure L), and, finally, between cortical dysplasias and ATRT (p = 0.0003) with respect to miR-96* ( Figure M).

miRNA Expression Profiling with Respect to Tumor Grade
Six miRNAs were found to be DE with respect to tumor grade. In particular, miR-1182 was found to be DE between controls and grade I tumors (p = 0.009) ( Figure 4A). Significant differences were also observed between controls and grade II tumors (p = 0.002), between controls and grade III tumors (p = 0.002), between grade I and grade III tumors (p = 0.004), between grade II and grade IV tumors (p = 0.005), as well as between grade III and grade IV tumors (p = 0.005) with respect to miR-183 ( Figure 4B), as well as significant differences were observed between controls and grade III tumors (p = 0.006) with respect to miR-194* ( Figure 4C). miR-3202 manifested a descending pattern from controls to grade IV tumors with significant differences between grade II and grade III tumors (p = 0.008) ( Figure 4D). Further on, significant differences were observed between grade I and grade II tumors (p = 0.008), as well as between grade II and grade IV tumors (p = 0.009) with respect to miR-4251 ( Figure 4E) and between controls and grade II tumors (p = 0.006) with respect to miR-592 ( Figure 4F). Finally, when including control samples, calculating the expression ratios significant differences were observed between grade I and grade III tumors (p = 0.0003), between grade II and grade III tumors (p = 0.007), as well as between grade III and grade IV tumors (p = 0.005) with respect to miR-4251 ( Figure 4G).

Descriptive and Association Statistics of Tumors Relations to Methylation
Tumor samples manifested a methylation pattern. In particular, methylated apoptotic genes included CASP8, RASSF1, MGMT, MSH6, GATA5, ATM1, TP53, CADM1, and RB1. Methylation was detected in one target gene, yet in some cases more than one genes were found to be simultaneously methylated. More specifically, it appeared that apoptotic gene methylation was significantly dependent on the samples' status that is if the sample was a neoplasm or a control (χ 2 = 16.19, p = 0.040). Similar dependence was observed with respect to malignant, benign tumors and controls (χ 2 = 30.90, p = 0.014), as well as with respect to MB, PA, ATRT, EP, and control samples (χ 2 = 65.86, p = 0.0004). Further on, similar dependence was manifested for tumor grade (χ 2 = 55.55, p = 0.0061). Interestingly, it appeared that there was no significant dependence observed with a particular methylated gene, but rather with the simultaneous presence of all detected methylated genes. The results of this analysis are summarized in Table 3.

Relations to Methylation
As in the case of gene methylation we have examined dependence of clinical parameters with respect to CNV. There were no significance dependence observed with any of the observed CNVs, which included duplication of CD27, deletion of CD6, deletion of CD6/CDKN2A/GATA5, and CD44 duplication. The results are summarized in Table 4.

Statistics of Total Number of Methylated Genes
Deceased patients (meaning that the patient did not survive at the time of the study) had significantly more methylated genes as compared to living patients, as well as deceased patients has more methylated genes as compared to control samples ( Figure 6A). The aforementioned observations were confirmed by the finding that there was a significant difference in the total number of methylated genes with respect to the tumor type that is if the tumor was malignant or benign and control samples. In particular, it appeared that malignant tumor types had significantly more methylated genes as compared to benign tumor types ( Figure 6B). No significance was observed between malignant tumors and control samples.

Methylation-Specific PCR
Since the most significant signal was obtained for the genes CASP8 and RASSF1 in order to test more samples, a methylation specific PCR (MS-PCR) protocol was designed [45]. The MS-PCR primers were designed using the MethPrimer software [46], and the PCR was performed at 55 • C. Representative results are shown in Figure 7. . DE miRNAs with respect to diagnosis. Significant differences were estimated using ANOVA with Bonferroni correction. Significant differences were observed between cortical dysplasia (n = 2) and ependymoma (n = 7), as well as controls (n = 13) and ependymomas (n = 7) with respect to miR-1234 (A). Significant differences were observed between controls and ependymoma (n = 7), as well as medulloblastoma (n = 16) and ependymomas (n = 7) with respect to miR-183 (B). Significant differences were observed between controls (n = 13) and cortical dysplasias (n = 2) with respect to miR-25* (C). Significant differences were observed between cortical dysplasia (n = 2) and pilocytic astrocytomas (n = 20), as well as cortical dysplasias (n = 2) and medulloblastomas (n = 16) with respect to miR-3675-5p (D). Significant differences were observed between controls and cortical dysplasias (n = 2) with respect to miR-612 (E). When including control samples, calculating the expression ratios significant differences were observed between medulloblastomas and ependymomas (n = 7) with respect to miR-122* (F). Significant differences were observed between cortical dysplasia (n = 2) and ATRT (n = 4) with respect to miR-130b (G). Significant differences were observed between ependymomas (n = 7) and ATRTs (n = 4) with respect to miR-2909 (H). Significant differences were observed between cortical dysplasias (n = 2) and ependymomas (n = 7) with respect to miR-302b (I). Significant differences were observed between pilocytic astrocytomas (n = 20) and ependymomas (n = 7) with respect to miR-4251 (J). Significant differences were observed between medulloblastomas (n = 16) and ependymomas with respect to miR-576-5p (K). Significant differences were observed between medulloblastomas (n = 16) and ependymomas (n = 7) with respect to miR-600 (L) and finally significant differences were observed between cortical dysplasias and ATRT with respect to miR-96* (  . DE miRNAs with respect to tumor grade. Significant differences were estimated using ANOVA with Bonferroni correction. Significant differences were observed between controls and grade I (n = 16) tumors (p = 0.009) with respect to miR-1182 (A). Significant differences were observed between controls and grade II tumors (n = 9) (p = 0.002), between controls and grade III tumors (n = 3) (p = 0.002), between grade I (n = 16) and grade III (n = 3) tumors (p = 0.004), between grade II (n = 9) and grade IV (n = 18) tumors (p = 0.005), as well as between grade III (n = 3) and grade IV (n = 18) tumors (p = 0.005) with respect to miR-183 (B). Significant differences were observed between controls and grade III (n = 3) tumors (p = 0.006) with respect to miR-194* (C). Significant differences were observed between grade II (n = 9) and grade III (n = 3) tumors (p = 0.008) with respect to miR-3202 (D). Significant differences were observed between grade I (n = 16) and grade II tumors (n = 9) (p = 0.008), as well as between grade II (n = 9) and grade IV (n = 18) tumors (p = 0.009) with respect to miR-4251 (E). Significant differences were observed between controls and grade II (n = 9) tumors (p = 0.006) with respect to miR-592 (F). Finally, when including control samples, calculating the expression ratios significant differences were observed between grade I (n = 16) and grade III tumors (n = 3) (p = 0.0003), between grade II (n = 9) and grade III (n = 3) tumors (p = 0.007), as well as between grade III (n = 3) and grade IV (n = 18) tumors (p = 0.005) with respect to miR-4251 (G) (Legend: R: Relative expression, where R = E tumors E controls and E is the expression levels of controls and tumors, respectively, C: Controls).

miRNA Differential Expression with Respect to Methylated Genes
Differential expression was also examined with respect to gene methylation. In particular, we have tested the presence of significant differences with respect to the methylation status of CASP8 and RASFF1 genes. In particular, miR-128, miR-183, miR-3202, miR-302e, miR-4307, miR-4330, and miR-491-5p manifested significant differences with respect to CASP8 methylation ( Figure 8A). Similarly, miR-128, miR-3202, miR-4251, miR-4307, and miR-576-5p manifested significant differences with respect to RASFF1 methylation ( Figure 8B).  Table 3. Relations of methylated genes with respect to sampling. The data represents the number (n) of samples found to be methylated in each subcategory (genes separated with '/' imply the simultaneous methylation of all genes). In particular, we have investigated the relation of methylation with respect to diagnosis and tumor grading. The table depicts the number of samples found to have their promoter methylated with respect to their diagnosis, i.e., between (a) neoplasms and controls; (b) malignant, benign, and control tissues; (c) medulloblastoma, astrocytoma, ATRT, ependymoma, and controls; and (d) tumor staging.

Hierarchical Clustering (HCL)
Unsupervised two-way hierarchical clustering (HCL) with Euclidian distance did not discriminate accurately between all tumor types and the normal control group, as well as between PAs and controls or EPs and controls (Figure 9).

Chromosomal Distribution of DE miRNAs
The first step in gene expression investigation, was the estimation of chromosomal distributions. In particular, we have investigated miRNA expression using the "natural" values, as well as mRNA ratios. It appeared that the most active chromosome was chromosome 22 ( Figure 10A) and, in particular, location 22q11.21 ( Figure 10B). It appeared also that the highest expression was mostly attributed to control samples ( Figure 10C,D). When accounting for miRNA ratios, chromosome 22 ( Figure 10E) and location 22q11.21 ( Figure 10F) were still the most active chromosome and location, respectively. Interestingly, miRNA expression appeared to be relatively equally distributed across all chromosomes ( Figure 10G,H)

K-Means Clustering
K-means clustering manifested four clusters ( Figure 11A), which is presented with the respective centroids ( Figure 11B) and sorted centroids ( Figure 11C). K-means clusters did not manifest any distinct patterns of expression, as well as sorted centroids did not manifest a certain pattern of expression. Figure 11. K-means clustering of DE miRNAs. K-means clustering manifested four clusters (A), which is presented with the respective centroids (B), and sorted centroids (C). K-means clusters did not manifest any distinct patterns of expression.

Functional Annotations of DE miRNAs Gene Ontology of DE miRNAs
Functional annotation included gene ontology (GO) annotation of DE miRNAs, which manifested major functions, such as angiogenesis, vasculature development, and developmental processes, such as tube morphogenesis and development ( Figure 12). No significant pathways were found for the annotated miRNAs and their gene targets.

Descriptive K-Means
Descriptive k-means were described in the Section 2. In search of a specific motif or pattern, we have clustered miRNAs with respect to tumor grade ( Figure 13A,B) and we have observed that DE miRNAs manifested in all cases a specific pattern, i.e., grade IV tumors manifested lower levels of expression, followed by and ascending order from grade I to III and cortical dysplasia manifesting the highest expression values ( Figure 13C). This observation involves probably a trigger mechanism, where DE miRNAs express to similar levels with control samples in near-control tissues (such as cortical dysplasias) and if down-regulated they are connected to aggressive tumor types (IV). The "trigger" mechanism arises probably from the fact that tumor grades I to III manifest an ascending pattern, indicating that those miRNAs operate with certain thresholds.  A,B). The first representation concerns the k-means classification of the miRNA expression data (A), whereas it is followed by the k-means centroids (B). The centroids in (B) represent the mean values of the k-means clusters depicted in (A). Those centroids (depicted in (B)) were sorted in an ascending order, in an attempt to find patterns related to tumor staging (C). Indeed, an interesting pattern was revealed with respect to tumor grade, where k-mean sorted centroids (C), showed a descending pattern from grade I to III, as well as an ascending pattern from grade III to I (C) (Legend: I: grade I (n = 16); II: grade II (n = 9); III: grade III (n = 3); IV: grade IV (n = 18); CD: Cortical dysplasia (n = 2)).
Further on, we have also performed k-means clustering with respect to patient survival ( Figure 14A,B). Interestingly, all DE miRNAs manifested an ascending pattern when moving from patients who achieved clinical remission, to relapse and deceased patients ( Figure 14C). This result indicated that the DE miRNA's expression levels followed a motif depending on the patient's survival status.

ROC Analysis
ROC classification manifested several factors that were able to discriminate between clinical parameters. This type of analysis is able to manifest significant correlations between two variables with respect to a classifier. Only significant classified variables are presented. This type of analysis, basically can indicate that a variable can be distinguished between the parameters of a classifier. As anticipated, Ki-67 expression manifested significant classification potential by discriminating between alive and deceased patients (survival) (Figure 15A), as well as between relapsed patients and patients in clinical remission (outcome) ( Figure 15B). Interestingly, miR-582 was also able to classify between outcome ( Figure 15C) and survival ( Figure 15D). Similarly, miRNA expression was able to discriminate the methylation status of specific genes and in particular miR-1246 ( Figure 15E) and miR-489 ( Figure 15F) were able to classify this status for CASP8, as well as miR-1246 ( Figure 15G) and miR-3614 ( Figure 15H) were able to classify for the methylation status of RASFF1.

Naïve-Bayes Analysis
Pairwise clustering of DE miRNAs was examined with a naïve-Bayes classification.
Pairs of DE miRNAs were tested for their classification potential and we found that miR-130b and miR-3672 ( Figure 16A), miR-106b and miR-130b ( Figure 16B), and miR-147 and miR-3672 ( Figure 16C) were able to discriminate between CD, ATRT, and EP. Similarly, miR-302b and miR-320e ( Figure 16D), miR-147 and miR-183 ( Figure 16E), and miR-516b and miR-95 ( Figure 16F) were able to discriminate between CD, tumor grade II, and tumor grade III. Although, we have tested all possible combinations of DE miRNAs for all clinical parameters, only few classifying miRNAs were obtained. Interestingly, the identified miRNAs, could separate ATRT, CD, and EP, while MB and PA manifested similar clusters, i.e., close expression values.

Common miRNA Signatures
DE miRNAs manifested a common pattern of expression. In particular, 31 miRNAs were found to be globally down-regulated in all Timor samples ( Figure 17A Figure 17C). MIR34A up-regulation was followed by MIR320E with over-expression in 68.75% of all tumor samples ( Figure 17D). MIR34A is estimated to have approximately 900 mRNA targets. Out of these, 17 genes were differentially expressed in our patient cohort. In particular, these genes were HYAL3, KIAA1210, MYO1C, FAM162B, OPN4, TP53INP2, ZNF281, TTC19, CRHR1, RAD9B, PAX8, JMJD1C, PPP1R11, PDGFRA, SCNN1G, SHKBP1, and ELMOD1. Similarly, MIR320E has 833 predicted mRNA targets, whereas six are differentially expressed in our patient cohort. More specifically, those genes included POLR1C, SIAH3, PLEKHA4, PTP4A1, TLL1, and RCN2. We have further searched for functional annotations of those genes, but no significant results were obtained, probably due to the fact that the gene sample size was small. Figure 15. ROC analysis. ROC analysis manifested significant classifiers with respect to clinical parameters. In particular, as expected Ki-67 expression, appeared to classify significantly between alive and deceased patients (described here as survival) (A), while it was also able to significantly classify between relapse and clinical remission (described here as outcome) (B). Further on, miR-582 appeared to discriminate between clinical remission and relapse (C), as well as between alive and deceased patients (D). In addition, CASP8 methylation status (i.e., if the gene is methylated or not), was significantly classified by miR-489, while in the case of RASFF1 methylation miR-1246 (G) and miR-3614 (H) were able to classify between methylated and un-methylated RASFF1.   A,B). On the other hand, one miRNA (miR-34a) was found to be globally up-regulated in~70% of all tumors (C), while miR-320e was the next miRNA with~68% up-regulation in all tumors (D).

Discussion
In the present study, we have performed high throughput experimentation, which included miRNA microarrays, mRNA microarrays and methylation studies with MS-MLPA. Our basic aim was to detect both differential, as well as common miRNA expression patterns and connect them to the expressed mRNA patterns and in addition, find their relations to specific gene methylation. Our analysis, revealed a total of 75 differentially expressed miRNAs between all CNS tumors and the control cohort. Amongst them, 31 were globally down-regulated and two were globally up-regulated in >68% of all tumors. More specifically, miR-34a and miR-320e were found up-regulated in~68% of all tumors when compared to the control group. The present observation was in agreement with our previous report on embryonal tumors, where we also found that miR-34a is probably globally up-regulated in tumors. MiR-34a is considered a tumor suppressor gene that it is known to regulate SIRT1 expression (silent information regulator 1). SIRT1, is known to be an oncogene, since it is a key-regulator of tumor suppressor proteins, such as the transcription factor p53 [66][67][68]. Consequently, one might have expected that the elevated expression levels of miR-34a observed in our study in both tumor types when compared with the control group would result in suppression of the SIRT1 expression, leading to apoptosis of cancer cells. A potential response for this came from Yamacuchi et al. (2008), who previously proposed that overexpression of miR-34a does not entirely suppress SIRT1 translation, possibly because it does not exactly match its SIRT1 binding site [67]. MiR-34a overexpression has also been observed in childhood ependymomas [69], pilocytic astrocytomas [70] and in low-and high-grade astrocytomas of childhood [71], suggesting potential global oncogenic roles in pediatric brain malignancies. Similarly, miR-34a regulates PDGFRA [72][73][74] and PAX8 [75], two genes which are known to be involved in childhood CNS tumors, as especially PDGFRA has been studied for its role in gliomas [72][73][74]. In general, there are still very few reports on the connection between miR-34a and those two genes. This indicates that there is a large field of research still open in order to comprehend CNS tumor biology. Yet, another interesting gene-target of miR-34a is the TP53 gene. Although there are numerous reports on the role and correlation between those two molecules, there is only one report concerning its role in the CNS. In particular, it is reported that TP53 is a direct target of miR-34a participating in microglia behavior and suppression of neuro-inflammation [76]. Hence, miR-320e and miR-34a might afford potential biomarkers related to inferior prognosis or even suggest possible global therapeutic targets.
Noteworthy, miR-34a overexpression observed in the present study is in agreement with previous studies, such as the study of Costa et al. (2011) [69], who reported that miR-34a was found to be highly overexpressed in ependymomas, however, in this study it was proposed that miR-34a is linked to a more favorable prognosis. In addition, we have found that miR-320e was significantly up-regulated in our patient population, highlighting a potential for this gene. According to our findings it can be not only as a possible diagnostic biomarker, but also a prognostic factor of a more favorable clinical outcome.
Through our analysis, using a ROC classifier, we have found a set of miRNAs that could pose potential biomarkers for tumor outcome, i.e., clinical remission or relapse. Interestingly, miR-582 was found to be differentially expressed between tumor samples and controls and it appeared to be able to separate samples with respect to survival and therapeutic outcome. The identification of miR-582 is described for the first time for childhood CNS tumors, since there are no previous reports concerning that miRNA. Moreover, we observed that miR-1246, miR-489, and miR-3614 could discriminate between the methylation status of CASP8 and RASFF1.
Interestingly, out of the globally deregulated miRNAs, miR-3202, miR-4251, and miR-4270, were found to be significantly different with respect to gender. To the best of our knowledge, there are no reports concerning the role of miRNAs in gender-specific CNS tumor ontogenesis and progression. Yet, few reports, which drew our attention, indicated that there is a gender-specific tumorigenesis for gliomas. In particular, it has been reported that females had a predominance in developing CNS tumors in case of previous cancer familial history, indicating a possible hereditary gender-specific risk [142]. Another interesting perspective of gender-specific visualization, has to do with the possible personalized treatments for patients. In particular, it has been reported that the adenine-to-inosine "inosinome" [143], is a potent gender-specific glioblastoma stratifier [144]. Although, there are no previous reports on the role of miRNAs with respect to gender, this approach could prove useful in towards a personalized treatment for CNS tumors.
Ki-67 protein is present during all active phases of the cell cycle making it an excellent marker to predict cell proliferation [145]. In the current setting we investigated the potential interactions between Ki-67 'positive' and miRNA expression patterns, in order to unravel and characterize the role of miRNAs underlying pediatric embryonal brain tumors. All potential miRNA oncogenes and tumor-suppressive genes that emerged from previous correlations between miRNA expression profiles and disease progression or patient clinical outcome were yet again further demonstrated as elevated and decreased expression levels were observed in the Ki-67 'positive' group of patients versus control tissues, respectively. In addition, identical findings were manifested regarding their putative prognostic role, either favorable or inferior. Yet, again, our prediction regarding their prognostic significance in pediatric embryonal CNS neoplasms was confirmed.
In the present work, we have also attempted to identify methylation and CNV patterns in childhood CNS tumors. Since the first reports on the role of DNA methylation and its significance in epigenetic regulation, numerous reports have highlighted its role in cancer. In particular, several reports have outlined its significance in pediatric CNS tumors [146]. In the present work, we have found that methylation played a role with respect to its number of methylated genes and not with respect to the individual methylated genes. This finding, to the best of our knowledge, is reported for the first time and it implies that epigenetic regulation is significant in tumor progression and it is probably the result of a multifactorial process.
Previous reports have indicated that supratentorial primitive neuroectodermal tumors (PNET), as well as ATRT manifest an aberrant RASSF1A methylation but not CASP8 [147]. In the present work, in a single ATRT sample, we have found methylation in the MGMT gene. This finding was different from a previous retrospective observational study, where it was found that most ATRT cases did not manifested a MGMT methylation [148].  [135][136][137][138][139][140][141] In our work, we have found that five MBs manifested a CASP8 methylation while two did not manifest any methylation. Our finding was in agreement with previous studies, where CASP8 was also found to be methylated indicating a significant role of CASP8 in MB pathophysiology [149,150]. Epigenetic regulation of CASP8 signifies the role of pro-apoptotic molecules in tumor progression. At the same time six samples were found to be methylated in the RASSF1 gene, which is in agreement with previous studies which indicated the significance of RASSF1 in MB pathology [151][152][153]. It appears that RASSF1 methylation leads to RASSF1 deactivation indicating a significant role of RASSF1 to MB pathology [151].
In the case of PA tumors, most samples appeared to be unmethylated and only three samples were found to be methylated in the CASP8, RASSF1, and MGMT genes. To the best of our knowledge, there are no previous reports on gene methylation of pediatric PA. This result was confirmed with the finding of our study that methylation status was significantly related to tumor grade.
Similarly, the two EP tumors we investigated were found to possess a methylation on MSH6 and RASSF1 genes. There are no previous reports concerning the methylation status of pediatric EP tumors.
Interestingly, we have found that several miRNAs including miR-128, miR-183, miR-3202, miR-302e, miR-4307, miR-4330, and miR-491-3p manifested significant differences between the methylation status of CASP8 and RASSF1. Although these miRNAs were not found to consist of a target for CASP8 and RASSF1, it was interesting to observe that those differences were manifested irrespectively of tumor type. To the best of our knowledge there are no previous reports on the correlation of those miRNAs to the methylation status of CASP8 and RASSF1 in CNS tumors.
Our understanding of the exact miRNA mechanisms in childhood CNS tumors' machinery is still limited. Thus, the comprehension of these mechanisms, regarding CNS tumor pathogenesis, could reveal candidate therapeutic targets. Further on, miRNAs could be used as possible biomarkers, for the prognosis, diagnosis, and treatment of childhood CNS tumors, but their role still remains largely unexplored.
MicroRNAs offer insights to many processes of the human body and are fascinating molecules to investigate. However, the study of miRNAs often involves hybridizationbased microarray technologies, a high-throughput technology, that may generate a large opportunity for errors when used to test miRNAs. Due to these limitations, all experiments must be regulated and controlled, to reduce the chances of error regarding the data produced. In addition, other technologies must be implemented as a way to confirm the results of the microarrays. A commonly used technique is RT-PCR, which amplifies specific genes as a way to validate the microarrays.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/cancers13215491/s1, Figure S1: Flow chart of the present work; Figure S2: Indicative analysis results as they are obtained from the special software; Table S1: The complete list of all differentially expressed miRNAs.