Low Expression of FOXP2 is Associated with Better Prognosis in Glioblastoma.

FOXP2 expression has been associated with the prognosis of some tumors, but the role of FOXP2 in glioblastoma has not been studied in-depth until now. The aim of the present work is to study the role of FOXP2 as a prognostic biomarker in glioblastoma. This is a retrospective observational case series study in which the expression of FOXP2 has been analyzed both at the protein level (immunohistochemistry) and at the mRNA level (RNAseq, in a cohort of glioblastoma patients from The Cancer Genome Atlas [TCGA] database). Other molecular and clinical data have also been included in the study, with special focus on miRNA expression data. Survival analysis using Log-Rank test and COX-regression have been used. Non-parametric statistical tests were also used to study differences between low and high FOXP2 expression groups. Patients with a high expression of FOXP2 protein showed a worse prognosis than those patients with low expression in both, progression free survival (PFS) (HR=1.711; p=0.034) and overall survival (OS) (HR=1.809;p=0.014). These associations were still statistically signicant in multivariate analysis. No prognostic association was found with FOXP2 RNA expression. Interestingly, two miRNAs that target FOXP2 (hsa-miR-181a-2-3p and hsa-miR-20a-3p) showed an interaction effect on OS with FOXP2 expression. A low level of these miRNAs expression was associated with a signicantly worse prognosis in patients with high FOXP2 RNA expression. Higher expression of FOXP2 at the protein level is associated with a worse prognosis. This protein expression may be regulated by the expression of specic miRNAs that target FOXP2 mRNA: hsa-miR-181a-2-3p and hsa-miR-20a-3p.


Introduction
Glioblastoma is the most common primary tumor in the adult brain. It is a tumor with an ominous prognosis, with a mean overall survival of approximately one year after diagnosis [1,2]. Despite the progress in the treatment of most neoplasms, advances in the management of glioblastoma are still scarce. Because of this, new strategies and approaches should be performed to gain a better understanding of glioblastoma biology. Bearing this in mind, the study of the molecular means used by the normal brain to achieve plasticity (e.g. FOXP2) could help to identify new molecular pathways implicated in glioblastoma development and recurrence. This knowledge would help to improve the management of this tumor.
FOXP2 is a transcription factor that belongs to the FOX family (forkhead/winged-helix). It is endowed with activating and, most often, repressing transcriptional activities [3]. It is highly expressed in the adult brain and during fetal central nervous system development [4][5][6]. Furthermore, this gene seems to be essential for the development of language [7,8]. In fact, mutations in this gene have been associated with congenital aphasia syndromes [9].
The expression of FOXP2 in some tumors has previously been reported. On the one hand, high expression of FOXP2 has been associated with a bad prognosis in prostate cancer [10], neuroblastoma [11] and multiple mieloma [12]. On the other hand, other studies have reported an association between FOXP2 expression and a better prognosis in breast cancer [13], liver cell carcinoma [14], osteosarcoma [15] and gastric cancer [16]. Thus, whether FOXP2 is a tumor suppressor or tumor-promoting gene is matter of controversy. In any case, the role of FOXP2 in the development and progression of glioblastoma has not been completely studied to date. A recent work reports that low expression of FOXP2 was associated with better survival outcomes and down regulation of this transcription factor resulted in inhibition of glioma cell proliferation [17]. Other studies also report a tumor-promoting effect of this transcription factor in glioma cell cultures [18,19]. However, more work is needed to establish the potential role of FOXP2 as a prognostic biomarker in glioblastoma.
The aim of the present study is to assess the role of FOXP2 expression as a prognostic biomarker in primary glioblastoma. Two cohorts of glioblastoma patients were analyzed. Firstly, the expression of FOXP2 was determined at the protein level using immunochemistry in a group of patients treated in our Neurosurgery Department. Secondly, the expression of FOXP2 at mRNA level was analyzed in a cohort of patients from the The Cancer Genome Atlas (TCGA) project.

Patients
Formalin-Fixed-Para n-Embedded (PFFE) tumor specimens were obtained from 62 patients ( TCGA data extraction Data from TCGA was downloaded from Firebrowse (http:// rebrowse.org/) (TCGA data version 2016_01-28). As mentioned above, only patients with RNAseq V2 data available were included in the study. Apart from RNA data, clinical, mutational, copy number variations, miRNA expression and methylation data (for determination the O6-methylguanine-DNA-methyltransferase [MGMT]) from the selected patients were also downloaded and included in a new database. It should be noted that miRNA expression and methylation data were not available for all patients. Details of this data generation are described elsewhere [20,21].

IDH1 and IDH2 mutations identi cation
In order to speci cally select those patients with primary glioblastoma, the presence of isocitrate dehydrogenase 1 (IDH1) and/or isocitrate dehydrogenase 2 (IDH2) mutations were studied. In the case of the TCGA group, the mutational annotation for both genes were analyzed and those patients with any kind of mutation in these genes (7 patients with IDH1 mutation and 1 patient with IDH2 mutation) were excluded. In the case of the local cohort of patients, the most common mutation in IDH1 (R132H) and IDH2 (R172K) were tested by PCR ampli cation and pyrosequencing.

MGMT methylation analysis in TCGA patients
Methylation probes were available for 117 patients in the cohort of TCGA patients included in the present study. In 67 of them the methylation status was tested using the HumanMethylation27 (HM27) platform and the rest of the patients were tested with the HumanMethylation (M450) platform. Some patients presented data from both platforms. Thus, the absence of signi cant differences between those samples on the two In nium platforms was con rmed by calculating the p-value using a student-t test [21]. Data was then merged by averaging the beta-values of the CpG probes of interest. The methylation status of MGMT was determined as explained in other reports [22]. In brief, the beta values were transformed in mvalues using the following formula: The logit(y) was then calculated using the model proposed by Bady et al (2012), where only the M-value of two MGMT CpG islands is considered (cg12434587 and cg12981137): #logit(y) = 4,3215 + 0,5271*cg12434587 + 0,9265*cg12981137 As proposed by Bady et al. (2012) a probability cutoff of 0.358 was used which empirically maximized the sum of sensitivity and speci city [22].

MGMT methylation analysis in local patients
The methylation status of the four most informative CpG islands in the MGMT promoter was studied in DNA puri ed from glioblastoma areas dissected from PPFE tissue sections using a commercial kit for pyrosequencing (Qiagen).

Immunohistochemistry (IHC) analysis
The expression of FOXP2 was tested in the PPFE specimens (originating from surgical resections) of the local cohort of patients by using a speci c antibody (Sigma-Aldrich, ref. HPA000382). Positivity was evaluated by two independent observers using the ImageJ® software in hotspot areas. A semiautomatic protocol, validated for ki-67, was used. Spearman's correlation coe cient was calculated between the positivity percentages of the two observers. FOXP2% positivity expression was included in regression analysis. Furthermore, this value was dichotomized (using different cutoffs: p25, p50 and p75) and this variable was used for the Kaplan-Meier curves and the Log-Rank test analysis.
RNAseq data from TCGA patients FOXP2 RNA expression data were extracted as explained above and they did not have a normal distribution, and thus were transformed by using base 2 logarithm (Log2 (x + 1)), where x is the expression of FOXP2. This new variable (Log2FOXP2) was included in the univariate Cox regression analysis. Furthermore, FOXP2 RNA expression was also dichotomized (using different cutoffs: p25, p50 and p75) and this variable was used for the Kaplan-Meier curves and the Log-Rank test analysis.
Apart from FOXP2 RNA expression data, other gene expression pro les were also extracted to perform a molecular classi cation. As is widely described in previous reports, there are differences in gene expression in glioblastoma that allow their classi cation in proneural, neural, classical and mesenchymal transcriptomic subtypes [20,23]. Using the list of input genes that are highly expressed in each subtype (http://tcga-data.nci.nih.gov/docs/publications/gbm_exp/), an unsupervised hierarchical cluster analysis was performed (MORPHEUS, Broad Institute, http://software.broadinstitute.org/morpheus/) and each patient was assigned to a molecular subgroup by cutting the resulted dendogram (supplementary Fig. 1).

MicroRNA expression data from TCGA patients
Bearing in mind that protein translation of FOXP2 mRNA is intensively regulated by miRNA [24], a number of these were selected to analyze their prognostic value. Firstly, miRNA public databases (mirDB [http://www.mirdb.org/] and mirANDA [http://www.microrna.org/microrna/home.do]) were used to identify the miRNA that target FOXP2. Secondly, from this list of miRNA (supplementary table 2), those with a proven role in glioblastoma pathogenesis acting as oncogenic or tumor suppressor miRNA were selected [25]. Finally, the expression of the miRNA from this list that was available in the TCGA data was analyzed (supplementary table 3).
As most of the miRNA expression did not show a normal distribution, a dichotomization of each one was performed by using the median as cutoff. These variables were included in a strati ed survival analysis with FOXP2 expression groups (supplementary table 4).

DNA Copy-Number Variation and DNA mutations
Copy number variation and mutation annotation les were downloaded as indicated above. After excluding those patients whose RNAseq data was not available, the ten-most-common cancer related genes showing copy number alteration (CNA) and/or mutation were studied. These genes (supplementary tables 5 and 6) were identi ed from the whole TCGA glioblastoma cohort whose analysis is available at cBioPortal (https://www.cbioportal.org/). A comparative analysis of the distribution of these genetic events between the two groups of FOXP2 expression was performed.

Statistical analysis
In order to study the differences between low and high expression of FOXP2, the median was used to generate two groups of expression. This approach was performed for both cohorts of patients. Non parametric statistical tests were used (Mann-Whitney U for continuous variables and Chi-Square/Fisher exact test for discrete variables). Statistical signi cance was considered when the p-value < 0.05. However, bearing in mind the high number of comparisons during molecular analysis, a corrected p value was used for these variables using the False Discovery Rate (FDR) method. Differences were considered statistically signi cant when FDR < 0.1.
Univariate and multivariate Cox regression analysis were performed in both cohorts of patients. Clinical, radiological and molecular variables were included in the univariate analysis and those with a p-value < 0.1 were included in a multivariate model. Kaplan-Meier curves and the Log-Rank test were also used to study the differences in overall survival (TCGA and local group) and progression free survival (only the local group) between different FOXP2 expression groups. The statistical signi cance for survival analysis was considered when the p-value < 0.05.
Finally, miRNA expression data was also considered for post-hoc analysis. A strati ed survival analysis bearing in mind the expression level of the selected FOXP2-targeting miRNAs and FOXP2 expression group was performed.

Results
Expression of FOXP2 in the local cohort of glioblastoma patients.
The IHC analysis, using a speci c FOXP2 antibody, showed that the protein was mainly located in the nucleus and the cytoplasm. The mean percentage of positivity in hot-spot areas was 28.33% (SD=32.29), with an adequate correlation between the results of the two independent observers (Spearman`s rho= 0.897; p<0.001).
Prognostic evaluation of FOXP2 expression in the local cohort of glioblastoma patients.
The univariate Cox regression analysis showed that higher FOXP2 expression was associated with a worse prognosis in both progression free survival (PFS) (HR=1.711, 95% C.I. [1.040 -2.814]; p=0.034) and overall survival (OS) (HR=1.809, 95% C.I. [1.129 -2.900]; p=.014). The univariate Cox regression analysis also included other variables that had previously been associated with the prognosis of glioblastoma (tables 1 and 2). The association of FOXP2 IHC positivity with worse PFS and OS was still signi cant when this variable was included in a multivariate Cox regression analysis with those variables that showed statistical signi cance in the univariate Cox regression analysis (p<0.1) (tables 1 and 2).
Comparison between high and low FOXP2 expression groupsin the local cohort of glioblastoma patients.
Using the median (p50=14.0) of FOXP2 % positivity expression, a comparative analysis between patients with low FOXP2 expression (<=p50) and high FOXP2 expression (>p50) was performed. This cutoff was selected because it was considered the best way to show prognostic differences ( gure 1). No clinical, radiological or molecular differences were identi ed between both groups ( Comparison between high and low FOXP2 expression groupsin the TCGA cohort of glioblastoma patients. A comparison between patients with high and low FOXP2 RNA expression was performed using the median expression of FOXP2 as cutoff (p50=7.1668). Differences were identi ed in the molecular classi cation distribution, where the "high RNA FOXP2" expression group presented a higher percentage of patients classi ed in the neural and proneural groups (table 5). These groups presented higher FOXP2 RNA levels than the classical and mesenchymal subtypes, but this difference was only statistically signi cant for the comparisons neural vs. mesenchymal and proneural vs. mesenchymal (ANOVA; p=0.024 and p=0.030, respectively) (supplementary gure 2).
Regarding the broad molecular information available in the TCGA patients, additional comparisons between low and high FOXP2 expression groups were performed. These comparisons were focused on DNA copy-number variation (CNV) and DNA mutations. Firstly, the frequency of chromosomal gains or losses was analyzed in both FOXP2 RNA expression groups. It should be noted that twelve patients (19.7%) from the high FOXP2 RNA expression group showed losses in 15q. This genomic loss was signi cantly different to those in the low FOXP2 RNA expression group (uncorrected p value=0.038), but did not reach signi cance with corrected FDR values (FDR>0.1) (supplementary table 7). No other statistical differences in the frequency of chromosomal alterations between both groups were identi ed.
Secondly, DNA CNV analysis was performed. Analysis of focal ampli cations and deletions in the top-10 glioblastoma cancer-related genes with focal CNVs showed a larger number of ampli cations of CDK4 (24.7% vs. 9.9%; p=0.027; FDR>0.1) in the high FOXP2 RNA expression group (supplementary table 8). CDK4 is located in chromosome 12 (cytoband: 12q14.1) and it encodes a protein member of the Ser/The protein kinase family. This kinase is responsible for the phosphorylation of retinoblastoma gene product (Rb). On the other hand, more patients in the low FOXP2 RNA expression group showed a notable deletion of CDKN2B (59.2% vs. 42.5%; p=0.048; FDR>0.1) (supplementary table 8). The protein encoded by this gene is also involved in cell growth and in the control of the cell cycle.
Finally, the mutational burden of the selected TCGA patients was analyzed. The most common driver mutations among the selected genes in the TCGA cohort of patients affected PTEN, TP53 and EGFR (supplementary table 5). Comparisons between the incidence of mutations in the two FOXP2 RNA groups are reported in supplementary table 9. No differences in the incidence of PTEN and TP53 mutations in low and high FOXP2 RNA expression groups were identi ed, but more patients with a missense mutation in EGFR gene were identi ed in the low FOXP2 RNA expression group (38.0% vs. 16.4%; p=0.042). Only one patient presented a driver mutation in ATRX. PCLO gene showed a higher incidence of passenger mutations in the low FOXP2 RNA expression group (15.5% vs. 2.8%; p=0.016; FDR=0.16). This gene encodes a protein of the presynaptic cytoskeletal matrix and is involved in establishing active synaptic zones and in synaptic vesicle tra cking. Its role in glioma pathogenesis has not been fully studied.
Regarding the results described above and bearing in mind the intense regulation of FOXP2 expression by different miRNAs, post-hoc analysis using the expression of some miRNAs were performed. Two miRNAs that targeted FOXP2 were selected (see Methods and supplementary table 4 . This difference was statistically signi cant (p=0.020) ( gure 3b). On the other hand, patients with a high expression level of has-miR-20a-3p did not show differences in OS between low and high expression groups of FOXP2 (333.0 vs. 454.0 days; Log-Rank test; p=0.511) ( gure 3b).

Discussion
The aim of the present work was to study the role of FOXP2 as a prognostic biomarker in glioblastoma.
Higher expression of FOXP2 at the protein level was signi cantly associated with worse outcomes (PFS and OS). However, this result was not replicated when considering the level of FOXP2 mRNA expression. No signi cant association between the expression of FOXP2 at mRNA level and prognosis was identi ed for glioblastoma. The explanation for these contradictory results may lie in the in uence of certain miRNAs that target FOXP2. All of these ndings are discussed below.

FOXP2 expression in glioblastoma
The expression of FOXP2 in glioma have been previously reported [17][18][19]. This expression is higher than the FOXP2 expression in normal brain and the level of expression seems to increase with the tumor grade [19]. In the present study, the majority of the patients showed some expression of FOXP2, either at the RNA or protein level. The expression of FOXP2 in glioma cells may promote the development and recurrence of the tumor. For example, He Q et al (2017) found an upregulation of FOXP2 in U87 gliomaexposed endothelial cells [18]. This upregulation was associated with an enhanced viability, migration and tube formation of these cells. Thus, FOXP2 seems to promote the angiogenesis in glioma [18]. Another in vitro study showed that FOXP2 overexpression in U87 and U251 glioma cells increases their migration, proliferation and invasion, with reduced apoptosis rates [19]. Finally, in a recent study, Zhang H et al (2019) reported that low expression of FOXP2 slowed down the cell cycle and inhibited the proliferation of U251 glioma cells [17]. Interestingly, these authors also con rmed these results in vivo in a nude mouse model [17]. Furthermore, it should be noted that FOXP2 has also been associated with the promotion of migration and invasion of tumor cellsin breast cancer [26,27], hepatocellular carcinoma [14], prostate cancer [10] and colorectal carcinoma [28], among others. Therefore, results from in vitro studies suggest that FOXP2 may promote the malignant biological behavior of glioma cells.
In the present study, higher expression of FOXP2 was associated with two speci c transcriptomic glioblastoma subtypes: proneural and neural. Although some molecular differences were identi ed (e.g. a different incidence of mutations in speci c genes or CNVs), they did not reach statistical signi cance when a corrected p-value (FDR) was considered. Furthermore, no other clinical or radiological data were distributed differently between the expression groups of FOXP2 of the two cohorts. Therefore, FOXP2 is expressed in glioblastoma and its expression is not clearly associated with speci c molecular, clinical and radiological features (apart from the mRNA expression pro le).

FOXP2 as a prognostic factor
Bearing in mind the results of in vitro studies, one could expect that those glioblastoma patients with a higher expression of FOXP2 would have a worse prognosis. Indeed, the present work showed an association between a worse OS and PFS in glioblastoma patients with higher expression at protein level (Table 1 and 2, Fig. 1). Zhang et al (2019), using a different method (western-blot) and a different antibody (Abcam, ab16046), also reported a worse prognosis for patients with higher FOXP2 expression in terms of mortality and tumor recurrences [17].
However, no signi cant association between FOXP2 expression at mRNA level and prognosis was identi ed (Table 4, Fig. 2). It should be noted that in many biological systems, gene expression at mRNA and protein levels are not perfectly correlated. In fact, signi cant differences have been detected by others [29][30][31]. In many cases, this discrepancy is mediated by the effect of microRNAs (miRNAs). miRNAs are small, highly conserved, non-coding RNA molecules found in most eukaryotes, including in humans [32]. miRNAs have a central role in the regulation of gene expression and, in the context of neoplastic cells, they are involved in malignant transformation [33,34]. A single miRNA may target several dozen or even hundreds of messenger RNAs (mRNA) and can cause an altered phenotypic state at the cellular level [35].
They normally affect gene expression (gene silencing) via translational inhibition and/or mRNA degradation [36]. Furthermore, data suggest that miRNAs play an essential role in hallmark biological features of glioblastoma [25], and they can act as potential oncogenes or tumor suppressors in gliomas [25,37].
Post-transcriptional regulation of FOXP2 mRNA is mainly mediated by miRNAs [24]. More than 280 miRNA are predicted to target human FOXP2 [24]. In the present study, miRNA public databases (mirDB and mirANDA) were used to identify miRNA that target FOXP2 and, from this list, a selection of those miRNAs that have a proven role in the pathogenesis of glioblastoma was compiled (supplementary table   4). An interaction effect of hsa-miR-181a-2-3p and hsa-miR-20a-3p expression with FOXP2 RNA expression was identi ed. Higher FOXP2 RNA expression was associated with worse OS in each of the groups with low hsa-miR-181a-2-3p and low hsa-miR-20a-3p expression. On the contrary, no difference in OS between FOXP2 RNA expression groups was identi ed in the high hsa-miR-181a-2-3p and hsa-miR-20a-3p expression groups. In other words, these miRNAs may negatively regulate FOXP2 expression and only when they are in a low level of expression FOXP2 mRNA could they be translated to protein and FOXP2 expression would have an oncogenic effect. In this regard, one can consider that the lack of signi cant prognostic value of FOXP2 mRNA expression levels (as opposed to FOXP2 protein expression level) is associated with the effect of these miRNAs. Both, hsa-miR-181a-2-3p and hsa-miR-20a-3p would act as tumor suppressors. hsa-miR-181a has already been considered as tumor suppressor in the context of glioblastoma [25]. Low expression of hsa-miR-181a has been associated with chemotherapy resistance in glioblastoma [38] and is considered a down-regulated miRNA in glioblastoma [39]. On the other hand, hsa-miR-20a-3p has not been speci cally studied in gliomas (unlike the 5p strand, which has been reported to show an oncogenic role [40,41] but does not target FOXP2). hsa-miR-20a-3p has been reported as a downregulated miRNA in gastric cancer [42] it seems to negatively regulate the TGF-b1 pathway or vascular endothelial growth factor (VEGF) expression in different diseases [43,44]. These pathways are important in glioblastoma development [45]. Therefore, the previously reported activity of these miRNAs would be compatible with the ndings of the present work.

Limitations
Some limitations should be mentioned concerning the present study. Firstly, the lack of FOXP2 protein expression information in the TCGA database hindered the direct validation of the results of the local cohort of glioblastoma patients. Although, the combination of information from RNAseq and miRNA expression were used to solve this lack of protein data, con rmation of the prognostic role of FOXP2 (at the protein level) in larger cohorts of patients is recommended.
Secondly, no information about the functional role of FOXP2 in glioblastoma has been provided here. Previous studies have reported on this issue, but only with commercial glioma cell cultures [17][18][19].
Bearing in mind that FOXP2 is a transcription factor, ChiP-Seq approaches (among others) in primary cultures should be considered in future studies to analyze the speci c target genes of FOXP2 in the context of glioblastoma.

Conclusion
Higher expression of FOXP2 at protein level is associated with a worse prognosis in glioblastoma patients. In the present study, results from protein expression were reproduced when considering the expression of FOXP2 at RNA level together with the expression of hsa-miR-181a and hsa-miR-20a-3p, two of the multiple miRNAs targeting FOXP2. Slight differences between patients with high or low FOXP2 expression in molecular features were identi ed. The expression of FOXP2 and its functional role in the pathogenesis of glioblastoma should be evaluated in future studies.

Declarations
Ethics approval and consent to participate The study was approved by the Hospital Universitario de Canarias Ethics Committee according to the Declaration of Helsinki.

Consent for publication
Not applicable.

Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that the research was conducted in the absence of any commercial or nancial relationships that could be construed as a potential con ict of interest.    Kaplan-Meier curves of low and high FOXP2 RNA expression for overall survival.