Two hundred and forty-four patients (94 females; mean age 59.5 years (SD = 14.29)) were included from the TCGA database (only patients with, at least, reversed phase protein array [RPPA] data). To specifically select those patients with primary glioblastoma, the presence of isocitrate dehydrogenase 1 (IDH1) and/or isocitrate dehydrogenase 2 (IDH2) mutations were discarded. The mutational annotation for both genes was analyzed and those patients with any kind of mutation in those genes were excluded (1 patient with IDH1 mutation and no patients with IDH2 mutation). Furthermore, those patients whose probes came from a previously treated glioblastoma or this information was not available were also excluded (n = 10). Therefore, a cohort of 233 patients (92 females; mean age 59.7 (SD = 14.18)) formed the definitive database for further analysis. Clinical and molecular features of the patients included in the study are shown in Table 1.

Data from TCGA was downloaded from Firebrowse (http://firebrowse.org/) (TCGA data version 2016_01-28). As mentioned above, only patients with RPPA data available were included in the study. Apart from RPPA data, clinical, mutational, copy number variations, RNA expression and methylation data (for determination the O6-methylguanine-DNA-methyltransferase (MGMT)) from the selected patients were also downloaded and included in the database. It should be noted that this information was not available for all patients. Details of this data generation are described elsewhere [24, 25].

Methylation probes were available for 198 patients from the cohort of TCGA patients that were included in the present study. In 102 of these patients, the methylation status was tested using the HumanMethylation27 (HM27) platform and in 112 patients with the HumanMethylation450 (HM450) platform. Two patients (0.9%) presented data from both platforms. To combine the data from both platforms, the absence of significant differences was first confirmed between those samples on the two Infinium platforms by calculating the p-value using a student-t test [25]. Afterwards, data was merged by averaging the beta-values of the CpG probes of interest.

The methylation status of MGMT was determined as explained in other works [26]. In brief, the beta values were transformed in M-values using this formula:

# M = log2 (Beta/(1-Beta))

Then we calculated the logit (y) using the model proposed by Bady et al. [26], where only the M-value of two MGMT CpG islands are considered (cg12434587 and cg12981137):

# logit (y) = 4.215 + 0.5271 $\times$ cg12434587 + 0.9265 $\times$ cg12981137

As proposed by Bady et al. [26] a probability cutoff of 0.358 was used which empirically maximized the sum of sensitivity and specificity.

AR expression data were extracted as explained above and were analyzed as described elsewhere [27]. These data showed a normal distribution (Kolmogorov-Smirnov; p = 0.093), and, thus, was included in the univariate Cox regression analysis (see below). Furthermore, AR expression was also dichotomized (using the median) and this variable was used to compare groups of low ($\le$p50) or high ($>$p50) AR expression for Kaplan-Meier curves and Log-Rank test analysis.

Apart from AR expression, the rest of the RPPA data was also included in the database to identify those proteins which showed a positive or negative correlation with AR activity (see below).

The RNAseq data was available for 127 patients (54.5% of the selected cohort of patients). These data were used for two approaches. Firstly, the expression profiles of a selected set of genes were also extracted to perform a molecular classification. As has been widely described in previous works, there are differences in gene expression in glioblastoma that allow their classification in proneural, neural, classical and mesenchymal transcriptomic subtypes [24, 28]. 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 (MORPHEUS, Broad Institute, http://software.broadinstitute.org/morpheus/) was performed and each patient was assigned to a molecular subgroup by cutting the resulted dendrogram (Supplementary Fig. 1).

Secondly, in order to infer the activity of the AR, the expression profile of previously validated androgen response genes (ARGs) was included in the database [29]. Specifically, a list of 13 genes that are induced by the activation of AR in both HPr-1AR (normal prostate cell line) and LNCaP (prostate cancer cell line) cells [29] were included (Supplementary Table 1). The activity of the AR was defined by the quantification of the composite expression of this 13-gene signature in each sample. As in other works [30], a Z-score was computed for the expression of each gene in each sample by subtracting the pooled mean from the RNAseq expression values and dividing the result by the pooled standard deviation. The AR putative activity (called here AR-score) for each sample was then computed as the sum of the Z-scores of the ARGs signature. The AR-score showed a normal distribution (Kolmogorov-Smirnov; p = 0.866) and was included in the univariate Cox regression analysis (see below). Furthermore, the AR-score expression was also dichotomized (using the median) and this variable was used to compare groups of low ($\le$p50) or high ($>$p50) AR-score for Kaplan-Meier curves and Log-Rank test analysis.

Copy number variation (CNV) and mutation annotation files were downloaded as described above and analyzed as described elsewhere [27]. After excluding those patients whose RPPA data was not available, the ten-most-common cancer related genes showing copy number alteration (CNA) and/or mutation were studied. Those genes (Supplementary Tables 2,3) were identified from the whole TCGA glioblastoma cohort whose analysis can be found in cBioPortal (https://www.cbioportal.org/). Comparative analysis of the distribution of these genetic events between the two groups of AR expression and AR-score was performed. Furthermore, CNVs and mutations of the AR were also analyzed.

Statistics has been performed as described in previous reports [27]. In brief, nonparametric statistical tests were used (Mann-Whitney U for continuous variables and the Chi-Square/Fisher exact test for discrete variables) to identify differences between low and high AR expression groups, as well as between low and high AR-score groups. Statistical significance 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 significant when FDR 0.1.

Univariate and multivariate Cox regression analysis were performed. Clinical and molecular variables were included in the univariate analysis and those with a p-value 0.05 were included in a multivariate model. Kaplan-Meier curves and the Log-Rank test were also used to study the differences in overall survival between different AR expression and AR-score groups. The statistical significance for survival analysis was considered when the p-value 0.05.

Finally, a correlation analysis between RPPA data and the AR-score of each sample was performed. The 5 proteins with the highest (positive or negative) correlation coefficient were studied. The genes that encoded these proteins were surveyed by shared Gene Ontology (GO) Biological Process data using the g: Profiler (http://biit.cs.ut.ee/gprofiler/gost).

Reversed Phase Protein Analysis (RPPA) data showed a mean AR abundance of 0.027 (SD = 0.38) in the probes of glioblastoma. As expected, a significant positive but moderate correlation between AR at the protein and mRNA levels was found (Spearman’s Rho; Correlation Coefficient (CC) = 0.557; p 0.0001). No differences in mean AR expression was identified between women and men (0.047 vs. 0.015; p = 0.670). Using the median of AR protein abundance (p50 = –0.017), a comparative analysis between patients with low AR expression ($\le$p50) and high AR expression ($>$p50) was performed. Patients from the high AR expression group presented a higher rate of MGMT methylation status (33.0 vs. 47.5%; p = 0.043) (Table 2). No other significant difference was identified between these two groups (Table 2). Regarding the median overall survival (OS) in each group, patients with high AR expression showed a longer OS than low AR expression patients (15.43 vs. 13.47 months), but this difference did not reach statistical significance (Log-Rank; p = 0.124) (Table 2, Fig. 1). No difference in OS was identified between low and high AR expression groups in women and men (see Supplementary Fig. 2A).

Fig. 1.

Survival analysis in Androgen Receptor (AR) and AR-score groups. (A) Kaplan-Meier curve for overall survival in patients with high and low AR expression groups (using the 50 percentile as cutoff). (B) Kaplan-Meier curve for overall survival in patients with high and low AR expression groups (using the 50 percentile as cutoff).

As regards the broad molecular information available in the TCGA patients, additional comparisons between low and high AR expression groups were performed. These comparisons were focused on DNA copy-number variation (CNV) and DNA mutations. Determination of the frequency of chromosomal gains or losses (Supplementary Table 4), analysis of focal amplifications and deletions in the top-10 glioblastoma cancer-related genes with focal CNVs (Supplementary Table 5), and the mutational signature in the top-10 glioblastoma cancer related genes (Supplementary Table 6) in both AR expression groups were all performed. No significant difference in all these molecular data between groups were identified, even when considering a non-corrected p-value (Supplementary Tables 4–6).

The mRNA expression of 13 validated androgen response genes (ARGs) were used to create an index which putatively reflects the activity of the AR (called AR-score) (see Methods, Supplementary Table 1). The mean AR-score in the studied cohort of patients (127 patients) was -0.024 (SD = 5.74). No correlation between the AR-score and AR at protein level expression was identified (Spearman’s Rho; CC = 0.060; p = 0.500), or between AR-score and AR expression at RNA level (Spearman’s Rho; CC = –0.007 p = 0.934). Males presented a higher mean AR-score (0.14, SD = 5.89) than females (–0.32, SD = 5.53), but this difference did not reach statistical significance (p = 0.850).

Using the median of AR-score (p50 = 0.05), a comparative analysis between patients with a low AR-score ($\le$p50) and high AR-score ($>$p50) was performed. No significant differences were identified between AR-score groups in the clinical and molecular variables included in the study (Table 3). Regarding chromosomal aberrations, a higher proportion of patients with a loss of 6p was identified in the high AR-score group (8 vs. 1 patient (13.6% vs. 1.7%)). This genomic loss was significantly different (p = 0.022), but did not reach significance with corrected p-values (FDR $>$0.01) (Supplementary Table 7). On the other hand, patients in the low AR-score group presented a higher proportion of CNVs involving the CDKN2A (60.9 vs. 37.5%; uncorrected p = 0.013) and CDKN2B (60.9 vs. 37.5%; uncorrected p = 0.013) genes. However, these differences were not significant when corrected p-values were considered (FDR $>$0.01) (Supplementary Table 8). Finally, no differences in the incidence of top-10 mutation glioblastoma cancer related genes between the low and high AR-score groups were identified (Supplementary Table 9).

Finally, in the survival analysis, patients with a high AR-Score showed a significantly worse OS than low AR-Score patients (13.1 vs. 15.6 months; Log-Rank; p = 0.025) (Table 3, Fig. 1).

Univariate Cox regression analysis was performed to analyze the value of AR protein expression and AR-score as prognostic factors. The AR protein expression did not show any association with OS (HR = 0.868; 95% Confident Interval (C.I.) (0.600–1.255); p = 0.452) (Table 4). The AR-score was associated with a worse prognosis in OS (HR = 1.070, 95% C.I. (1.033–1.108); p = 0.000177). The univariate Cox regression analysis also included other variables that have previously been associated with the prognosis of GB (Table 4). The association of AR-score with a worse OS was still significant in the multivariate analysis (HR = 1.054, 95% C.I. (1.008–1.103); p = 0.020), where variables that showed statistical significance in the univariate Cox regression analysis (p 0.05) were included (Table 4). Furthermore, the activity of the AR in glioblastoma was associated with a worse prognosis in both genders, but it seems to have more impact in females than males (Supplementary Fig. 2B).

Correlation analysis between AR-Score and RPPA data was performed. The aim of this approach was to identify positive and negative significant relationships between AR activity and other cancer related proteins included in the RPPA data from the TCGA. Among the five proteins with the highest correlation coefficients, four had a positive relationship with the AR-score (Fig. 2): PAI-1 (CC: 0.51), followed by Caveoline-1 (CC = 0.43), phosphorylated NDRG1 (NDRG1 pT346) (CC = 0.42) and Caspase-8 (CC = 0.38). The other selected protein was phosphorylated (CDK1 pY15) and its expression had a negative correlation with the AR-score (CC = –0.39) (Fig. 2).

Fig. 2.

Correlation analysis between AR-Score and the top-5 correlation coefficient proteins in the RPPA from TCGA data. (A) Plots showing the relationship between the AR-score and each selected protein expression. (B) Table showing the information from the correlation analysis plotted in A.

In order to infer biological processes regulated by these proteins, the genes that encoded these 5 selected proteins were surveyed by Gene Ontology (GO) Biological Processes data. This analysis showed that these proteins are mainly associated with the regulation of apoptotic process (adjusted p = 0.0007).

As stated in the introduction section, many studies have reported a high expression of AR at mRNA [23] and protein levels [20, 22] in glioblastomas. Furthermore, the expression of the AR in gliomas has been related to the tumor grade, thus its expression increases as the tumor grade increases [22] and the activation of the AR seems to have a promoting effect in GB [20, 21, 23]. In this respect, Yu X et al. [20] reported that promoting the TGF$\beta$signaling pathway in GB by adding TGF$\beta$1 to a GB cell line with high AR expression, decreases the rate of cell growth, thereby increasing apoptosis. However, when they subsequently increased DHT levels, the effects of TGF$\beta$ in GB were inhibited, which would suggest that AR activation inhibits the effect of TGF$\beta$ in GB, thus the AR indirectly promotes cell growth and reduces apoptosis [20]. In 2018, Nomi Zalcman et al. [21] observed that adding AR inhibitors, like enzalutamide, to glioma cell cultures produced a dose-dependent cell death, as well as when transfecting GB cells with low interfering RNA targeting to the full-length AR-RNA. Finally, Rodríguez-Lozano et al. [23], found that the addition of testosterone to GB cell cultures increased GB cell proliferation, tumor cell migration and invasion through AR activation.

The present work tries to determine AR activity by calculating the AR-Score, an index calculated by using the expression of 13 androgen responsive genes (ARGs), which has been previously validated for prostate cancer. According to the results of previous studies, higher activation of the AR (here a higher AR-Score) was associated with a worse prognosis (Tables 3,4; Fig. 1). The difference in prognosis does not seem to be associated with other molecular or clinical variables and the worse prognosis of higher AR-Score was maintained in multivariate analysis. Furthermore, a recent study has described an implication of the androgenic pathway activation in the immune response against GB, facilitating the infiltration of regulatory T-cells which are immunosuppressive T-cells [31]. Therefore, one can assume the importance of considering the AR activation when studying its pathogenic role in GB, not only the AR expression at the mRNA and/or protein level.

A correlation analysis was performed between the AR-Score and RPPA data with the aim of identifying the relationship between the AR activity and proteins that have a recognized role in oncogenesis. Among the 5 proteins with the highest correlation with the AR-Score, four had a positive correlation (PAI-1, Caveoline-1, NDRG1 and caspase-8) and one of them had a negative correlation (phosphorylated CDK1).

PAI-1 (plasminogen activator inhibitor 1) is a glycoprotein which is synthesized by endothelial and tumor cells, among others. It belongs to the Serpin superfamily and acts as a protease inhibitor [32]. High levels of PAI-1 have been associated with tumor growth, necrosis development, micro-thrombosis and, consequently, with a lower OS in malignant brain tumors [33]. The expression of PAI-1 is regulated by the TGF$\beta$ signaling pathway [34], which, as previously stated, is modulated by the activation of the AR [20]. Furthermore, there are experimental results which have shown that PAI-1 can positively regulate the PI3K/AKT pathway, another signaling pathway implicated in GB pathogenesis [33]. Thus, the activation of the AR may promote the expression of PAI-1, which has an oncogenic role in GB.

The expression of Caveolin-1 (Cav-1) also showed a positive correlation with the AR activity. Cav-1 is a protein responsible for the formation of caveolae, which are complex plasma membrane structures involved in the signal transduction process, cell-cycle regulation and cell migration [35]. The relationship between Cav-1 expression and AR activity in GB has not been established until now, but studies have analyzed this relationship in prostate cancer. In this regard, Cav-1 seems to act as a coactivator of the AR, increasing the androgen/AR signaling [36, 37]. Although the association of Cav-1 expression and AR activity has not been previously studied until now, the role of Cav-1 in gliomas has already been analyzed. Firstly, higher Cav-1 expression in GB versus normal brain tissue has been reported [38, 39]. Secondly, the expression of Cav-1 in gliomas seems to increase proportionally with the tumor grade [35]. Finally, Cav-1 participates in multiple processes that promote GB development and facilitate tumor cells invasion [39].

As explained for PAI-1 and Cav-1, the positive relationship between the expression of NDRG1 and Caspase-8 with AR activity described in the present work are also supported by the literature. On the one hand, NDRG1 (N-Myc Downstream Regulated 1), is associated with an unfavorable prognosis in gliomas and prostate cancer [40]. NDRG1 expression in GB cells was mainly described in those cells exposed to hypoxia [40, 41]. Furthermore, NDRG1 overexpression in GB may inhibit the TFG-$\beta$ pathway [40] in a similar way as previously described for PAI-1. On the other hand, Caspase-8 (an initiator protease of the signaling cascade that leads to apoptosis) is also overexpressed in GB and its expression is associated with a worse prognosis [42, 43, 44]. Caspase-8, in the context of GB, promotes angiogenesis, tumorigenesis and the expression of various cytokines such as VEGF (vascular endothelial growth factor), IL-6 and IL-8 [43, 44]. IL-6 and IL-8 are particularly interesting because they promote the activation of the AR in an independent-ligand manner [45].

Finally, a negative correlation between the CDK1 (activated by phosphorylation, CDK1 pY15) and the AR activity was identified. CDK1 plays a fundamental role in the cell cycle by modulating mitosis onset [46, 47]. Phosphorylated CDK1 leads to a phosphorylation of the AR in prostate cancer (phosphorylation in Ser-308) and this modification of the AR has been associated with a decrease in its transcriptional activity [48, 49] and with a longer OS [50]. Although these previous reports in prostate cancer may explain the findings of the present work (i.e., a negative relationship between AR-Score and phosphorylated CDK1), no evidence about the phosphorylation status of the AR and its relationship with AR activity and/or prognosis in GB have been reported.

Overall, the relationship between the selected proteins and the AR activity is supported by previous reports on GB and other tumors (mainly in prostate cancer). However, more studies are needed to specifically analyze the interaction of these proteins and AR in the context of GB.

Some limitations should be considered in the present work. Firstly, AR activity has been determined by using the expression of an ARGs set that was previously validated in prostate cancer cells and normal prostate cells, but not in glioma cells. In this respect, new studies, using Chip-Seq methods, focusing on specific ARGs in gliomas are required. In any case, the ARGs used in this work should also be validated in primary GB cultures exposed to androgens.

Secondly, using the AR-Score as a measure of AR activity, one cannot infer whether the activation of the AR is ligand dependent or independent. Bearing this in mind, a recent work has reported an activation of the AR by ligand-independent signaling through EGFR [51]. Thus, future works should also be focused on this issue, both in vitro and in vivo, measuring blood androgen levels and evaluating the effect on tumor growth or recurrence.