Meta-Analysis and Experimental Validation Identified FREM2 and SPRY1 as New Glioblastoma Marker Candidates

Glioblastoma (GB) is the most aggressive brain malignancy. Although some potential glioblastoma biomarkers have already been identified, there is a lack of cell membrane-bound biomarkers capable of distinguishing brain tissue from glioblastoma and/or glioblastoma stem cells (GSC), which are responsible for the rapid post-operative tumor reoccurrence. In order to find new GB/GSC marker candidates that would be cell surface proteins (CSP), we have performed meta-analysis of genome-scale mRNA expression data from three data repositories (GEO, ArrayExpress and GLIOMASdb). The search yielded ten appropriate datasets, and three (GSE4290/GDS1962, GSE23806/GDS3885, and GLIOMASdb) were used for selection of new GB/GSC marker candidates, while the other seven (GSE4412/GDS1975, GSE4412/GDS1976, E-GEOD-52009, E-GEOD-68848, E-GEOD-16011, E-GEOD-4536, and E-GEOD-74571) were used for bioinformatic validation. The selection identified four new CSP-encoding candidate genes—CD276, FREM2, SPRY1, and SLC47A1—and the bioinformatic validation confirmed these findings. A review of the literature revealed that CD276 is not a novel candidate, while SLC47A1 had lower validation test scores than the other new candidates and was therefore not considered for experimental validation. This validation revealed that the expression of FREM2—but not SPRY1—is higher in glioblastoma cell lines when compared to non-malignant astrocytes. In addition, FREM2 gene and protein expression levels are higher in GB stem-like cell lines than in conventional glioblastoma cell lines. FREM2 is thus proposed as a novel GB biomarker and a putative biomarker of glioblastoma stem cells. Both FREM2 and SPRY1 are expressed on the surface of the GB cells, while SPRY1 alone was found overexpressed in the cytosol of non-malignant astrocytes.


Introduction
Glioma is the cancer of the brain connective tissue. Histological characteristics of tumor cells are used for the classification of gliomas in several types, the most common of them being astrocytomas and oligodendrogliomas. Furthermore, gliomas are graded according to their malignancy in four

Results and Discussion
We performed meta-analysis of genome-scale mRNA expression data from three unrelated public repositories. Four genes were identified as potential new CSP-encoding glioblastoma markers. After bioinformatic validation, two of them were selected for experimental validation and one of them was confirmed as a potential novel cell-surface glioblastoma biomarker. Its expression was also experimentally confirmed in GSC.

Selection of Glioblastoma Datasets
The repositories Gene Expression Omnibus (GEO) [16], ArrayExpress [17], and GLIOMASdb [18] were browsed in order to obtain datasets appropriate for selection or validation of new GB/GSC marker candidates (see Table 1). The single GLIOMASdb dataset complied with the requirements for appropriate datasets (see Materials and Methods). Search in GEO (query "glioblastoma", filters "homo sapiens" and "DataSets") yielded four appropriate datasets: GSE4290/GDS1962, GSE23806/GDS3885, GSE4412/GDS1975, and GSE4412/GDS1976. In ArrayExpress, the query »glioblastoma, grade« (with filters »homo sapiens«, »RNA assay«, and »array assay«) yielded three appropriate datasets (E-GEOD-52009, E-GEOD-68848, and E-GEOD-16011), the query »glioblastoma, stem cells« (same filters) yielded two datasets (E-GEOD-4536 and E-GEOD-74571), while the last query »glioblastoma, neurospheres« (same filter) produced no appropriate results.  The best of the appropriate datasets-i.e., the datasets that contained curated data and/or permitted multiple comparisons relevant for selection of new GB/GSC marker candidates (see Table 1)-were chosen to be used for selection of new marker candidates. These datasets were the GEO datasets GSE4290/GDS1962 and GSE23806/GDS3885, as well as the single GLIOMASdb dataset (see Table 2). The rest of the appropriate datasets (Table 2) were used for bioinformatic validation of novel GB/GSC candidates that were identified in the selection stage. Category: Grade-datasets with samples from gliomas of different grade; Neuro/GSCL-datasets with samples of neurospheres and/or GSCL; NPC = number of possible comparisons; Rank-NPC-based rank within category; *** If rank is equal, curated datasets have precedence.

Selection of New Glioblastoma Marker Candidates
We conducted eight pairwise comparisons between groups of samples from the selection datasets (i.e., GSE4290/GDS1962, GSE23806/GDS3885, and GLIOMASdb). The aim of these selection tests was to identify new GB/GSC marker candidates. Similarly, 13 pairwise comparisons were conducted for bioinformatic validation of new GB/GSC marker candidates between groups of samples from the seven independent validation datasets (i.e., GSE4412/GDS1975, GSE4412/GDS1976 , E-GEOD-52009,  E-GEOD-68848, E-GEOD-16011, E-GEOD-4536, and E-GEOD-74571). Additional data about the tests-e.g., the size and composition of the compared groups-are presented in Table 3. Detailed  test results for each test are available in the Tables S1−S22 (see Supplementary Materials), while the  Figures S1−S20 contain sample codes and boxplot distributions of sample Table 4 lists the top ten genes according to the combined results of the eight tests of the selection stage ( Table 3). The entire list is available in the Table S23. The selection of genes was based on false discovery rate (FDR)-adjusted p values (cut-off value p < 0.001), and on log2 of the ratio of average expression levels in the two compared groups (log change or lfc; cut-off value lfc > 1). Three of the top ten genes, including one CSP gene (CD276), met the criteria of all eight tests. Furthermore, seven other genes passed 7/8 tests, and three of them are CSP genes: SPRY1, FREM2, and SLC47A1. All these genes scored both p < 0.05 and lfc > 0 in the single failed test (where they had either p > 0.001 or lfc < 1), indicating that their expression was at least slightly elevated (lfc > 0) in the first test group and their p value was within conventionally used significance limits (p < 0.05). SPRY1, FREM2, and SLC47A1 were thus included along CD276 in further analysis.

Assessment of the Effectiveness of Selection Tests
Six established candidate genes coding for cell surface proteins (CD133-also known as PROM1, CXCR4, MELK, PTCH1, CD9, and CD44) were used as a benchmark for assessing the effectiveness of selection tests that identified SPRY1, FREM2, SLC47A1, and CD276 as new glioblastoma marker candidates. Housekeeping genes cytochrome c1 (CYC1), hypoxanthine phosphoribosyltransferase 1 (HPRT1), H3 histone family 3A (H3F3A), ribosomal protein lateral stalk subunit P0 (RPLP0), and endoplasmatic reticulum membrane protein complex subunit 7 (EMC7) were used as negative control. A points-based system (described in detail in Materials and Methods) was applied in order to indicate different levels of significantly elevated gene expression. A scale from four to zero points was used-four points for the most significantly elevated expression, zero points for non-significant expression differences or genes with decreased expression. Average points per test (APPT) was used as a performance measure.
Results of our selection tests (Table 5) show that all six established candidates scored three points or more in a majority of the tests and that all the established candidates (except PTCH1) scored on average more than two points in each test. These results show that our selection tests were able to expose the elevated expression of the established candidates. On the other hand, four of five housekeeping genes had APPT less than one and only HPRT1-purportedly universally reliable as negative control in cancer-related transcriptome analysis [23]-had APPT equal to one.

Validation Tests with the Data from Other Datasets
The comparisons used in selection tests that identified four new glioblastoma marker candidates were recreated with data from seven other datasets (Table 2), resulting in 13 additional validation tests (Table 3) of the same type as selection tests. The points-based system-described in Materials and Methods-was used to evaluate different levels of elevated expression. The established candidates and housekeeping genes from the previous section were also used here with the same purpose. Results of these validation tests are presented in Table 6.   APPT = average points per test (rounded to three decimals); 4 points: FDR-corrected p < 0.001 and lfc > 1 (highly significant and highly overexpressed); 3 points: FDR-corrected p < 0.05 and lfc > 0, while criteria for 4 points are not met (significant and overexpressed); 1 point: p < 0.05 and lfc > 0 (overexpressed but borderline non-significant after application of FDR correction); 0 points: either p > 0.05 or lfc < 0 (non-significant or not overexpressed); /: gene not featured on the platform or in the results table.  FREM2 had the highest point average (3.250) among all 15 featured genes, but it also had the highest number of missing results (5/13). In the second place among all twelve genes, CD276 had only slightly lower point average (3.167) than FREM2 and only one missing value. In fact, all four new candidates outperformed each of the established candidates. The worst among the established candidates was PTCH1, whose APPT (0.615) was lower than APPT of each control gene except EMC7 (0.545). Two of the housekeeping genes-H3F3A and RPLP0-had APPT equal to or slightly above 1, but their APPT was still notably lower than that of CD133 (1.917), which was after the outlier PTCH1 the second worst-placed established candidate. All new candidates and all established candidates, except PTCH1, had point averages generally similar to those in the selection tests. This proves that the results of selection tests were not a consequence of the particular dataset choice.

Novelty Assessment of the Identified Marker Candidates
A brief review of relevant literature (Table 4) revealed that CD276 had already been linked to carcinogenesis of glioblastoma and other cancers [19,24]. However, for the other three candidates (FREM2, SPRY1, and SLC47A1) no direct association with glioblastoma has been reported so far. Since SLC47A1 had the worst performance among these three genes in the previous stage (Table 6), it was omitted from experimental validation, which therefore included FREM2 and SPRY1. Their already known association with various diseases is discussed below.

FREM2-Expression Patterns and Association with Disease
FREM2 is a protein localized mostly in various epithelia or in the extracellular matrix, and is widely expressed in embryonic tissues. Smaller fractions of FREM2 are also present in the cytosol and endoplasmic reticulum, which is localized close to the nuclear membrane and is the place of protein synthesis. FREM2 mRNA has been localized, among other tissues, in the embryonic neural stem cells (NSC) that are predecessors of some brain and spinal cord regions [25]. However, there are no reports of FREM2 detection in NSC in the adult brain. Biological processes in which FREM2 is involved include epithelial formation, cell adhesion and communication, and development of heart and inner ear. FREM2 missense mutations have been associated with congenital conditions such as Fraser syndrome [26] and unilateral renal agenesis [27]. According to the GeneCards database [25], low levels of FREM2 expression at the protein level have been detected in the Ga-MG cell line one of conventional GB cell lines [25].

SPRY1-Expression Patterns and Association with Disease
SPRY1 is a protein with cellular location dependent on cell activation via epidermal growth factor (EGF). In cells not stimulated by EGF, it is localized mainly in the cytoplasm, while in stimulated cells it migrates to the edge of the outer cell membrane where it becomes anchored as a peripheral membrane protein. SPRY1 has been reported to inhibit the differentiation of murine NSC into neurons, so it most likely acts in the same way with regard to human embryonic NSC, where its expression has already been confirmed [28]. There are no reports of SPRY1 detection in NSC in the adult brain [29].
Beside its pathological role in non-malignant legius syndrome-a condition similar to neurofibromatosis-and alongside its involvement in glioma (which will be discussed later) SPRY1 has pathological role in non-malignant Legius syndrome and has been implicated in some malignancies, such as adrenal cortical adenocarcinoma, prostate cancer, and rhabdomyosarcoma (RMS) [29,30]. In RMS, SPRY1 is present in cancer stem cells [30].

Reported Relations between Glioblastoma Multiforme and FREM2/SPRY1
FREM2 and SPRY1 were briefly mentioned in relation to GB in a few studies. However, none of these studies was focused on identification of new GB/GSC biomarkers. Ivliev et al. [31] identified EGFR gene amplification as a genetic signature in glioma and suggested its cause was the EGF signaling deregulation mediated by Sprouty family genes. Since SPRY1 is a member of this family, it may be implicated alongside other related genes in glioma carcinogenesis. Likewise, FREM2 gene amplification [32] has been discovered in gliosarcoma, a mixed tumor with a glioblastoma multiforme-resembling glial component and a sarcoma-like mesenchymal component. The mesenchymal part contains glial cells that have acquired mesenchymal properties by mutation. FREM2 amplification is part of genetic signature of these cells and may be causatively linked to glial-mesenchymal transformation [32].
In a study by Forte et al. [33], SPRY1 was among a number of genes with elevated expression in the undifferentiated glioblastoma-derived neuropsheres as compared to the differentiated glioblastoma-derived primary cultures. An unrelated study Aldaz et al. [34] demonstrated that SPRY1 was over-expressed at the protein level in undifferentiated cells as opposed to cells that had undergone four days of differentiation. Expression of SPRY1 protein was decreased in undifferentiated cells transfected with microRNA-21 (miR-21)-a small non-coding RNA that supposedly acts as a tumor suppressor-and the silencing of SPRY1 mRNA by miR-21 was subsequently confirmed as the mechanism behind this decrease [34].
Berezovsky et al. [35] knocked down the expression of Sox2 and reported (among many other genes) decreased expression of FREM2 in Sox2-negative glioblastoma-derived neurospheres, as well as decreased expression of SPRY1 in Sox2-negative serum-differentiated glioblastoma cells. Recently, Okawa et al. [36] analyzed the proteomes and secreted proteomes (secretomes) of NSC and cells obtained from four GSCLs. FREM2 was found to be among 447 proteins differentially expressed between the proteomes of NSC and GSCL, as well as among 138 proteins differentially expressed between the corresponding secretomes. However, this study did not reveal which of the differentially expressed genes had elevated expression and which had decreased expression in GSCL. On the other hand, SPRY1 expression levels demonstrated no significant difference between the proteomes, and SPRY1 was not even detected in the secretomes.
SPRY2, a gene from the Sprouty family-related but not identical to SPRY1-has already been proposed as a general biomarker of glioblastoma cells and a potential therapeutic target [37]. SPRY2 is a known regulator of oncogenic receptor tyrosine kinases (RTK), which are involved in signaling processes that regulate glioblastoma cell tumorigenicity and survival. Analysis of glioblastoma specimens [37] confirmed SPRY2 expression at the protein level in tumor tissue and further demonstrated that SPRY2 was over-expressed in glioblastoma tumors expressing EGFR variant III (EGFRvIII), a mutant EGFR variant that is constitutively active and whose presence in the tumor is deemed to be a negative prognostic indicator. As the Sprouty family genes are closely related, the function of SPRY1 in glioblastoma carcinogenesis might resemble that of SPRY2. SPRY2 was also reported to be one of extracellular matrix components that prevent EGFR downregulation [38], and constitutive EGFR activation was in turn implicated as one of necessary steps for the transformation of non-malignant NSC into malignant GSC [39]. Non-malignant glial cells such as astrocytes can also give rise to GSC after EGFR activation [39]. Thus, elevated expression of extracellular matrix components such as Sprouty proteins and FREM2-whose amplification in glial cells has already been linked to their transformation into mesenchymal cells in gliosarcoma [32]-may be involved in the NSC-GSC and astrocytes-GSC transformations via prevention of EGFR deactivation. Gene expression of FREM2 and SPRY1 was measured in the mature glioblastoma cell lines U251MG and U87MG, in the stem-like NCH cells-a co-culture of the glioblastoma stem-like cell lines NCH644 and NCH421K -as well as in neural stem cells (NSC) and astrocytes. mRNA was successfully isolated from all cell lines with mean concentration 1357 ng/µL and A 260 /A 280 ratio from 1.82 to 1.99. Primer efficiency for FREM2 was 1.71 and for SPRY1 1.87. Statistical analysis of the results was performed in GraphPad Prism 6. All data sets were tested for normality (Gauss distribution) and outliners were identified. In the cases where no expression was observed (e.g., FREM2 in U87MG and astrocytes) zeroes were added to the graph data sets. Outliers are presented on the graphs, but were removed using the 'Identify outliers' option of column analysis in GraphPad Prism 6 software and were thus excluded from the statistical analysis (one-way ANOVA, unpaired t-tests, Kruskal-Wallis, and Mann-Whitney tests).
Multiple group comparisons revealed a difference in FREM2 expression between NCH and U87MG, as well as between NCH and NSC ( Figure 1A, left). For SPRY1, significantly different expression was observed between U251MG and non-malignant neural stem cells ( Figure 1A, right), as well as between both conventional glioblastoma cell lines, U251MG vs. U87MG ( Figure 1A, right). Mann-Whitney test showed change in expression of SPRY1 in NCH/U251MG/U87MG compared to NSC, whereas no change was observed for FREM2 for the same analyzed groups ( Figure 1B). Statistical analysis with unpaired t-test with Welch's correction highlighted significant difference in FREM2 expression between NCH and the mature glioblastoma cell lines U251MG and U87MG ( Figure 1C, left). Comparison between NCH and NSC ( Figure 1D) revealed that both FREM2 and SPRY1 are differentially expressed in malignant (NCH) when compared to non-malignant (NSC) stem cells.
In a separate experiment, we analyzed all glioblastoma cell lines (NCH, U251MG, and U87MG) in combination with astrocytes as reference cells. Multiple group comparisons showed significant difference in expression of FREM2 in NCH compared to astrocytes, between NCH and U87MG and between the two mature glioblastoma cell lines U251MG and U87MG (Figure 2A, left). For SPRY1, multiple group comparisons showed strong difference in expression of all glioblastoma cell lines (NCH, U251MG, and U87MG) compared to astrocytes where higher expression was observed (Figure 2A, right). For FREM2, Mann-Whitney test showed difference in expression of the combined glioblastoma cell lines (NCH, U251MG, and U87MG) compared to astrocytes, and unpaired t-test with Welch's correction also showed change in expression of SPRY1 for the tested groups. When comparing NCH to both U251MG and U87MG mature cell lines, significant change in expression was observed for FREM2 only ( Figure 2C). At last, comparison of NCH and astrocytes showed significant change in expression for both analyzed genes FREM2 (overexpression in NCH) and SPRY1 (overexpression in astrocytes) ( Figure 2D).  These results indicate that FREM2 was significantly over-expressed in glioblastoma stem-like cell line (NCH) in comparison to mature glioblastoma cells of the U87MG cell line representing relatively less aggressive neuronal type GB [40], and in comparison to neural stem cells and astrocytes. Its expression was also higher in the NCH cell line than in mature glioblastoma cells of the U251MG cell line representing most aggressive mesenchymal type GB [40]. In addition, FREM2 was found significantly over-expressed in NCH in relation to the combination of mature GB-related cells of U87MG/U251MG, and significance of its over-expression was also found in U251MG compared to U87MG. Since the difference in FREM2 expression between NCH and U251 was not significant, although FREM2 expression was higher in the stem-like NCH cell line, this marker candidate may be considered a common indicator of aggressiveness of glioma cells.
SPRY1 was expressed in all glioblastoma-related cell lines. SPRY1 expression was significantly higher in the U251MG cell line towards both the U87MG and NSC cell lines (see Figure 1A, right). Its expression in glioblastoma-related malignant cells (NCH, U251MG, and U87 grouped together) was significantly higher than in the non-malignant NSC cell line ( Figure 1B, right). However, it was also remarkably expressed in neural stem cells and especially in astrocytes. For both genes, relatively high expression levels were indicated in the U251MG cell line, which recapitulates the most typical pathobiological features reported for human GBM [40]. For FREM2, these experimental findings also confirm the results obtained by meta-analysis, namely the results of Test 5 (GSCL vs. conventional glioblastoma cell lines), where FREM2 expression was significantly elevated in GSCL, with p < 0.001 and lfc > 1 (see Table 5).

Immunofluorescent Staining
Gene transcription levels are not necessarily matched by the expression of corresponding proteins, which are the actual factors of physiological functions in the cell. Thus, with immunofluorescent detection, we validated the expression of FREM2 and SPRY1 proteins in the glioblastoma cell lines U251MG and U87MG compared to primary astrocytes as a control or a model of healthy environmental brain tissue.
Experimental results are presented in Figure 3A,B. Cells were fixed in formaldehyde and immunofluorescently stained with commercially available anti-FREM2 ( Figure 3A) and anti-SPRY1 ( Figure 3B) antibodies. For better contrast, the left part of the image is shown in black and white, and the overlay in the original staining is shown on the right. Both target proteins are stained red, while the nuclei are blue. In the U251MG cell line, FREM2 was found highly expressed and clearly localized on the cell membrane, but also present close to the nuclear membrane (location of ER with ribosomes as the site of protein synthesis) and in cytosol, which may reflect significant overexpression of this protein in the aggressive U251MG cell line. Both cell membrane and cytosolic localization for FREM2 are in concordance with the data from the COMPARTMENTS database [25], while PROTEINATLAS (which also contains data from the U251MG cell line) [41] highlights its cytosolic and UNIPROT [42] its membrane localization. Differences in data from individual databases are probably due to different biological sources (cells, tissues). In our case, FREM2 was less membrane-expressed in the less aggressive U87MG cell line, and much less in non-malignant astrocytes, where no evidence of its membrane localization was found. SPRY1 was also found to be localized on the membrane of U251MG cells while no evidence of membrane localization was found in U87MG and non-malignant astrocytes.
These results highlighted FREM2 and potentially also SPRY1 as GBM marker candidates that may be useful for therapeutic targeting since they are highly expressed on the surface of the most typical GBM cells while not expressed on the surface of non-malignant glial cells which are found in tumor surroundings.

Western Blot Validation
In order to support and quantify the results obtained with immunofluorescent staining and considering qRT-PCR findings, and to more precisely define if either of the markers is more specific to glioblastoma stem cells, we performed Western blot analysis on all cell lines used in this study. The results-which for FREM2 broadly support the previous results from qRT-PCR analysis and immunofluorescent staining-are presented in Figure 4. FREM2 was significantly overexpressed in the aggressive mesenchymal U251MG cell line in comparison with non-malignant astrocytes, while no such observation could be made for the less aggressive neural U87MG cell line. FREM2 expression was also higher in NCH when compared to both U251MG and U87MG, even when considering confidence intervals, but these differences were not statistically significant. The results for FREM2 expression in astrocytes confirm the results previously obtained with RT-qPCR and immunofluorescence, which both showed very low expression levels that correlate with the literature data that indicate the absence of this marker in the adult brain [29]. Conversely, the results for SPRY1 confirmed its substantial expression at the protein level in astrocytes, which was also demonstrated at the gene expression level with qRT-PCR.
At the protein level, both FREM2 and SPRY1 were expressed in NSC, with no significant difference to NCH, although in the case of FREM2, the expression in NCH was remarkably higher than in NSC. Expression of both potential marker candidates in non-malignant neural stem cells is not unexpected, since it is also a common trait of established cell-surface GSC marker candidates, such as CD133. Western blot analysis confirmed qPCR results that revealed the presence of FREM2 and SPRY1 in neural stem cells and in stem-like glioblastoma cells. In addition, while expression of SPRY1 was found to be relatively high in all glioblastoma-related cell types, FREM2 expression was remarkably higher in GSC compared to U251 and especially U87, and higher than in NSC.
Western blot results indicate that FREM2 is promising GBM marker candidate for typical GB tumors, represented by the U251MG cells, but is less suitable as a marker for atypical, less aggressive GB tumors, represented by the U87MG cells. FREM2 could also be a potential GSC marker candidate even though expression differences in the comparisons NCH vs. U251MG and NCH vs. U87MG fell short of being statistically significant. Thus, our experimental validation confirmed FREM2 as a novel glioblastoma marker candidate, and additionally a potential GSC marker candidate. On the other hand, Western blot analysis did not support the findings of immunocytochemistry for SPRY1, which however still remains as a marker candidate to distinguish GB-related cells (lower expression) from non-malignant astrocytes (remarkably higher expression) and neural stem cells (lower expression) at the qRT-PCR level.  NCH421K). U251, U87-protein extracts from (respectively) glioblastoma U251MG and U87MG cell lines. NSC-protein extract from neural stem cell cells. ASTROCYTES-protein extract from primary human astrocytes. Quantification of Western blotting for FREM2 and SPRY1 expression: relative band intensities were calculated as the ratios between the antigen band intensities for FREM2 and SPRY1 and those of the whole protein loaded on the gel (both measured in arbitrary units [AU]). GAPDH was used as loading control for both western blots. Statistical analysis of Western blot. Mean protein expression value is presented on the graph with error bars representing standard deviation (SD). x-axis: samples; y-axis: relative protein expression. * p ≤0.05, ** p < 0.01. Kruskal-Wallis test with Dunn's correction: NCH vs. ASTROCYTES: FREM2-significant difference (p = 0.0013); U251 vs. ASTROCYTES: FREM2-significant difference (p = 0.0405); U251MG vs. U87MG: SPRY1-significant difference (p = 0.0027).

Selection of Glioblastoma Datasets
In order to find appropriate datasets, we searched the data repositories GEO (NCBI, Bethesda MD, USA) and ArrayExpress (EMBL-EBI, Hinxton, UK), as well as the GLIOMASdb database (Beijing Neurosurgical Institute, Beijing, China).
We devised seven pairwise comparisons relevant for selection of new GB/GSC marker candidates ( Table 1). Each dataset had to contain groups of samples that would permit at least one of these comparisons. Additionally, each compared group had to contain at least 10 samples (therefore, we considered only datasets with at least 20 samples, as this is the minimal number that would theoretically permit comparison of two groups with 10 samples each). This cut-off sample size is derived from the findings of Stretch et al. [43] that the results of differential gene expression analysis based on less than 20 samples (in two equal-sized groups) lack prediction accuracy.
GEO and ArrayExpress datasets that fulfilled the requirements from the previous paragraph (subsequently referred to as appropriate datasets) were divided into two categories: (i) datasets with samples from gliomas of different grade, and (ii) datasets with samples of neurospheres and/or GSCL. Within each category, datasets were ranked by the number of comparisons relevant for GB/GSC marker selection ( Table 1) that can be made by using data from a particular dataset. If two or more datasets enabled the same number of comparisons, priority was given to curated datasets (i.e., GEO datasets marked with "GDS"), for which NCBI assures that background processing and measurement are consistent across the dataset, that all value measurements are calculated in an equivalent manner, and that samples are biologically and statistically comparable [44]. The best datasets from each category were chosen to be used for the selection stage, while all other appropriate datasets were used for bioinformatical validation. The ranking of datasets obtained by using the above criteria is presented in Table 2. The single GLIOMASdb dataset was also chosen for use in the selection stage since it is curated, enables one of the relevant comparisons, and has a large size (325 samples).

Differential Gene Expression Analysis with GEO2R
In each test used for selection or validation (see Table 3)-except in the Test 8, which was based on GliomasDB, where results had already been provided by the dataset authors [18]-we used the program GEO2R (NCBI, Bethesda MD, USA) to find differentially expressed genes. The program GEO2R [45]-based on the limma R package [20]-was used to perform a series of two-tailed t-tests (see Table 3), which used variance stabilization [46] and FDR-corrected p values [47]. Additional data about the tests-such as detailed lists of samples included in each of the compared groups for each test, and graphical presentations of distribution of sample values (boxplots)-are available in the Supplements.
We used the cut-off values of p < 0.001 for FDR-adjusted p values and lfc > 1 (indicating the gene expression ratio of at least 2:1 between the two compared groups) for log change. A gene had to fulfill both the p value and lfc requirements in order to be listed as differentially expressed in an individual test. If the analyzed dataset included data from alternative probes that hybridize to the same gene transcript, a gene was deemed to fulfill the requirements if at least one of the hybridizing probes scored appropriate lfc and p values. The lists of differentially expressed genes in each individual test were compared by the Draw Venn Diagrams [48] application (University of Ghent, Belgium). Genes were checked for localization of proteins that they encode (see the next subsection) and only those genes that encode cell-surface proteins (CSP genes) were taken into further consideration. CSP genes differentially expressed in 8/8 tests immediately qualified for further consideration, while CSP genes that passed 7/8 tests qualified if they had both FDR-adjusted p < 0.05 and lfc > 0 in the single failed test.

Analysis of Localization Data of Corresponding Proteins with COMPARTMENTS
The database COMPARTMENTS [49] was used to obtain quantified localization data for proteins encoded by the genes that had been highlighted by the selection procedure. This database provides confidence scores that a protein can be found at a particular location, e.g., in the plasma membrane. A score of 5 denotes the highest probability, while a score of 0 indicates a complete lack of evidence for this location [49]. For each protein that corresponds to one of highlighted candidate genes, we looked for its COMPARTMENTS score indicating the probability of its localization at the membrane. If this score was high-i.e., 5 or 4-the protein was deemed to be a cell-surface protein (CSP).

Assessment of the Effectiveness of Selection Tests
We were interested in the performance of established GSC markers candidates in our selection tests, since established marker candidates are expected to demonstrate significantly elevated expression in any relevant test. Six established CSP-encoding GSC marker candidates (CD133/PROM1, CXCR4, MELK, PTCH1, CD9, and CD44) were used as positive control to check the relevance of selection tests. Housekeeping genes known for their stable expression either across all cancer types (HPRT1 [23] and EMC7 [50]) or particularly in glioma (CYC10 [51], H3F3A [52,53], and RPLP0 [52]) were used as negative control.
A class-based point system was used to reflect different levels of significance (measured by p values) and elevated expression (measured by lfc). Based on p values and lfc, genes were classified into four classes: (i) a gene had both FDR-adjusted p < 0.001 and lfc > 1; (ii) a gene did not meet the first criteria but still scored both FDR-adjusted p < 0.05 and lfc > 0; (iii) a gene did not meet the second criteria but still had lfc > 0 and p < 0.05; (iv) a gene failed all the three criteria-i.e., it had either p > 0.05, or lfc < 0, or both. The classes (i) and (ii) contain genes with statistically significant expression differences under conditions of multiple testing (FDR-correction of p values), while the classes (iii) and (iv) contain genes with non-significant differences after FDR-correction is applied.
The system was designed to give additional weight to classes (i) and (ii), which contain genes with statistically significant expression differences after FDR-correction. Genes from the class (i) received 4 points, those from the class (ii) received 3 points, those form the class (iii) received 1 point, while genes from the class (iv) received 0 points. Since some test results for a particular gene may be missing (in most cases because the gene is not featured on the microarray platform that was used to generate the data on which the test is based), we used average points per test (APPT) instead of sum of points as the performance measure for tested genes. When the APPT for an individual gene was calculated, tests with missing values for the gene were disregarded.

Validation Tests with the Data from Other Datasets
These tests were based on the seven relevant comparisons listed in the section Selection of glioblastoma datasets. Test numbers refer to numbers used in this listing: e.g., test v3-2 is the second (2) validation (v) test based on comparison 3. These tests used data from all appropriate datasets (see Selection of glioblastoma datasets) that were not chosen for the selection stage ( Table 2). Classification of test results (points-based system), as well as use of established candidates and negative controls were the same as in the performance section of established marker candidates in selection tests.

Novelty Test of Genes Identified as Potential Markers
To assess the novelty of genes identified as potential markers, for each gene revealed by the selection procedure, we used the query ("name of the candidate, glioblastoma") in the Google Scholar browser (Google, Mountain View CA, USA; accessed on 14 November 2017) and browsed the top 20 results. If this search revealed that a gene had already been linked to glioblastoma carcinogenesis, the gene was excluded from further consideration.
The outliers were removed using the 'Identify outliers' option of column analysis in GraphPad Prism 6 software. The ROUT method, which was recommended by the program and can identify one or more outliers, was used. The ROUT method is based on a false discovery rate (FDR). The Q value, which is the maximum allowed FDR, was set to 1%, which means that there was a 1% chance of falsely discovering an outlier. The graphs as presented in Figures 1 and 2 contain all data that were obtained, including the identified outliers. However, these values were excluded during the statistical analysis (one-way ANOVA, unpaired t-tests, Kruskal-Wallis, and Mann-Whitney tests).
3.6.3. Immunofluorescent Staining 6 × 10 4 U251MG cells (CLS) [58] were seeded on coverslips and incubated for 24 h at 37 • C and 5% CO 2 . 1 × 10 5 U87MG cells (ATCC) [59] and 1 × 10 5 astrocytes (Ixcells Biotech) [54] were seeded on poly-D-Lys (Sigma-Aldrich) pre-treated coverslips and incubated for 48 h at 37 • C and 5% CO 2 . Afterwards cells were washed three times and fixed in 4% formaldehyde for 15 min at room temperature. The fixative was washed with PBS three times for 10 min. Cells were permeabilized with 0.1% Triton X-100 for 15 min. Subsequently, cells were incubated for one hour in 1% BSA. After removal of blocking buffer, anti-FREM2 antibody (SAB3500517 produced in rabbit; Sigma-Aldrich) was added in dilution 1:30 and anti-SPRY1 (WH0010252M1 produced in mouse; Sigma-Aldrich) was added in dilution 1:25. The cells were incubated overnight at 4 • C. Afterwards cells were washed three times for 10 min with PBS and incubated with anti-rabbit secondary antibody (1:250 dilution; SAB4600399, Sigma-Aldrich) and anti-mouse antibody (1:250 dilution; SAB4600396, Sigma-Aldrich) for one hour at room temperature. Cells were washed three times for 10 min with PBS, stained with 300 nM DAPI, washed again three times for 10 min in PBS and then mounted on slides with ProLong Glass Antifade (ThermoFisher). Axio Imager M2 (Zeiss, Oberkochen, Germany) equipped with ZEN software was used to acquire images. Images were analyzed using ImageJ [62].

Western Blot Analysis
From all cultured cell lines proteins were extracted using commercially available kit ProteoExtract ® Transmembrane Protein Extraction Kit (Merck, Darmstadt, Germany and protein concentration was determined with BCA assay. The expression of FREM2 and SPRY1 at the protein level in the cell lines U251MG, U87MG, NCH, NSC and primary human astrocytes was examined using Western blot analysis. Approximately 20 µg of protein lysate from each cell line was denaturated in 2% SDS, 60 mM dithiothreitol, 60 mM Tris-hydrochloric acid (Tris-HCl, pH 6.8), and 0.1% bromophenol blue, and loaded onto a NuPAGE 4% to 12% Tris-glycine gel. After SDS-PAGE proteins were transferred onto polyvinyldidene fluoride (PVDF) membranes. The membranes were stained using Ponceau S to show total protein amount loaded on the gel and quantified with the Multi Gauge 3.2 software (Fuji, Tokyo, Japan). The PVDF membranes were blocked for 1 h at room temperature in PBS containing 5% non-fat dried milk, and incubated overnight at 4 • C in 1% PBS-milk containing the primary antibody (rabbit polyclonal anti-FREM2, 1:500 dilution Novus Biologicals, Littleton, CO, USA; or mouse monoclonal anti-SPRY1, 1:500 dilution; Sigma-Aldrich). The membranes were washed in 0.1% PBS-Tween, incubated with the secondary antibody (anti-rabbit/mouse IgG HRP-conjugated antibody; 1:5000 dilution; Sigma Aldrich) for 1 h at 4 • C, and washed in 0.1% PBS-Tween. Membranes were treated with SuperSignal West Pico Chemiluminescent Substrate (ThermoFisher). Bands were visualized using a LAS-4000 CCD camera (Fuji) and analyzed with the Multi Gauge 3.2 software (Fuji).
Experiments of Western blot were performed in triplicates. The data were analyzed for statistical significance using one-way ANOVA followed by Dunn's multiple comparison tests. P values < 0.05 were considered to indicate statistical significance. The means ± standard deviation were calculated for all variables.

Conclusions
Meta-analysis of genome-scale expression data from three unrelated public repositories identified four genes as potential new CSP gene markers of glioblastoma multiforme. After bioinformatic validation, two-FREM2 and SPRY1-qualified for experimental validation at the transcriptional and protein levels. At the protein level of expression" FREM2 was significantly over-expressed in the GBM-relevant cells in comparison to non-malignant brain cells, which was also confirmed at the gene level by qRT-PCR. In addition, FREM2 was over-expressed at the protein level in glioblastoma stem-like cells in comparison to mature GBM cells as well as neural stem cells. Moreover, qRT-PCR results proved significantly elevated expression of FREM2 in glioblastoma stem-like cells when compared to mature GBM cells. These results confirmed FREM2 as a novel GBM marker candidate and also as a potential GSC marker candidate. On the other hand, SPRY1 was not confirmed as a GSC marker candidate at the protein level. However, its strong expression in astrocytes in relation to the studied GB-related cell lines highlights its potential value to distinguish GB-related cells from the non-malignant tissue.
Meta-analysis may have some disadvantages, such as heterogeneity of input data, different ways of data normalization, and the use of groups with a low number of samples. These shortcomings were avoided by using curated data that were uniformly normalized and by not using very small groups in our statistical comparisons. The results of meta-analysis have been experimentally validated with three different in vitro methods, which also have their own limitations, such as the heterogeneity of cell lines used for in vitro experiments. In our study, there was a disproportion between the results of gene/protein expression in both conventional GBM cell lines. However, in the U251MG cell line, which recapitulates the most typical pathobiological features reported for human GBM, the experimental results for FREM2 were consistent with the findings of bioinformatic analysis. It is also difficult to compare gene/protein expression in cell lines with different growth speed, such as GSCL (e.g., NCH) and NSC. Unlike the GSCL, NSC cells grow very slowly and large number of passages is needed in order to acquire enough material for comparative expression studies. The NSC population likely became heterogeneous due to slow growth and consisted of unevenly differentiated cells, and this prevented determination of statistically significant differences in FREM2 expression at the protein level.
In conclusion, we identified two new cell-surface glioblastoma marker candidates, FREM2 and SPRY1. At the protein level, which matters the most for potential therapeutic use of markers, FREM2 was significantly overexpressed in the GBM-relevant cells in comparison to non-malignant brain cells and is thus proposed GBM marker candidate. In addition, FREM2 was overexpressed in glioblastoma stem-like cells in comparison to mature GBM cells as well as neural stem cells which highlights this protein as a putative biomarker of glioblastoma stem cells. As concerning SPRY1, despite Western blot results, it may also be interesting as a therapeutic target since immunofluorescent staining, as in the case of FREM2, indicated its presence on the membranes of GB cells belonging to the aggressive mesenchymal type, while it was not present on the surface of non-malignant astrocytes, where it was expressed exclusively in the cytosol and close to the nuclear membrane. Thus, SPRY1 could also be used as a target for therapeutic antibodies, which act by binding with cell-surface antigens and do not penetrate in cell interior. In order to assess the potential value of FREM2 and SPRY1 as GB/GSC-specific therapeutic targets, further research is needed at the level of immunohistochemistry and animal models. In this context, combination of both markers may be considered in research towards a multitargeting anti-glioblastoma therapeutic strategy.  Table S12: Results  of Test v2-3; Table S13: Results of Test v3-1; Table S14: Results of Test v3-2; Table S15: Results of Test v3-3; Table  S16: Results of Test v3-4; Table S17: Results of Test v3-5; Table S18: Results of Test v4-1; Table S19: Results of Test v5-1; Table S20: Results of Test v6-1; Table S21: Results of Test v7-1; Table S22: Platform decoding for tests v2-3 and v3-5; Table S23: Results of Venn diagrams analysis; Figure S1: Selection of samples and boxplot for Test 1; Figure S2: Selection of samples and boxplot for Test 2; Figure S3: Selection of samples and boxplot for Test 3; Figure S4: Selection of samples and boxplot for Test 4; Figure S5: Selection of samples and boxplot for Test 5; Figure S6: Selection of samples and boxplot for Test 6; Figure S7: Selection of samples and boxplot for Test 7; Figure S8: Selection of samples and boxplot for Test 8; Figure S8: Selection of samples and boxplot for Test v1-1; Figure S9: Selection of samples and boxplot for Test v2-1; Figure S10: Selection of samples and boxplot for Test v2-2; Figure S11: Selection of samples and boxplot for Test v2-3; Figure S12: Selection of samples and boxplot for Test v3-1; Figure S13: Selection of samples and boxplot for Test v3-2; Figure S14: Selection of samples and boxplot for Test v3-3; Figure S15: Selection of samples and boxplot for Test v3-4; Figure S16: Selection of samples and boxplot for Test v3-5; Figure S17: Selection of samples and boxplot for Test v4-1; Figure S18: Selection of samples and boxplot for Test v5-1; Figure S19: Selection of samples and boxplot for Test v6-1; Figure S20: Selection of samples and boxplot for Test v7-1. qRT-PCR quantitative real-time polymerase chain reaction RMS Rhabdomyosarcoma RPLP0 ribosomal protein lateral stalk subnit P0 RTK receptor tyrosine kinase