Glioma-Associated Microglia/Macrophages Display an Expression Profile Different from M1 and M2 Polarization and Highly Express Gpnmb and Spp1

Malignant glioma belong to the most aggressive neoplasms in humans with no successful treatment available. Patients suffering from glioblastoma multiforme (GBM), the highest-grade glioma, have an average survival time of only around one year after diagnosis. Both microglia and peripheral macrophages/monocytes accumulate within and around glioma, but fail to exert effective anti-tumor activity and even support tumor growth. Here we use microarray analysis to compare the expression profiles of glioma-associated microglia/macrophages and naive control cells. Samples were generated from CD11b+ MACS-isolated cells from naïve and GL261-implanted C57BL/6 mouse brains. Around 1000 genes were more than 2-fold up- or downregulated in glioma-associated microglia/macrophages when compared to control cells. A comparison with published data sets of M1, M2a,b,c-polarized macrophages revealed a gene expression pattern that has only partial overlap with any of the M1 or M2 gene expression patterns. Samples for the qRT-PCR validation of selected M1 and M2a,b,c-specific genes were generated from two different glioma mouse models and isolated by flow cytometry to distinguish between resident microglia and invading macrophages. We confirmed in both models the unique glioma-associated microglia/macrophage phenotype including a mixture of M1 and M2a,b,c-specific genes. To validate the expression of these genes in human we MACS-isolated CD11b+ microglia/macrophages from GBM, lower grade brain tumors and control specimens. Apart from the M1/M2 gene analysis, we demonstrate that the expression of Gpnmb and Spp1 is highly upregulated in both murine and human glioma-associated microglia/macrophages. High expression of these genes has been associated with poor prognosis in human GBM, as indicated by patient survival data linked to gene expression data. We also show that microglia/macrophages are the predominant source of these transcripts in murine and human GBM. Our findings provide new potential targets for future anti-glioma therapy.


Introduction
Malignant glioma (WHO grade III and IV), such as glioblastoma multiforme (GBM), are highly aggressive and account for almost 50% of all brain neoplasms. As of today no successful treatment exists, offering glioma patients an average survival time of about one year after diagnosis, despite surgical tumor resection, radio-, and chemotherapy [1]. Glioma are highly invasive and the cellular and genetic inter-and intra-tumor heterogeneity has lead to the failure of current anti-glioma therapy [2,3,4,5].
Among other immune cells, brain-resident microglia and peripheral macrophages/monocytes are attracted toward glioma in large numbers and they can amount up to 30% of the cells in the tumor tissue [6,7,8]. Experimental findings by us and others show that a tumoricidal activity of glioma-associated microglia/macrophages (GAMs) is counteracted by the glioma cells, which reprogram GAMs into tumor supportive cells. GAMs actively support tumor growth, e.g. by secreting matrix-degrading enzymes as well as immunosuppressive factors [9,10,11,12]. Thus, these stromal cells might serve as an attractive target for anti-glioma therapy and the identification of factors produced by GAMs might help to better understand disease progression [11,13,14,15,16].
Based on in vitro activation states macrophage polarization is divided into classical inflammatory (M1) activation (responsible for Th1 responses, type I inflammation, killing of intracellular pathogens) and alternative (M2) activation. Alternative activation can be further divided into M2a (Th2 responses, type II inflammation, killing of pathogens, allergy), M2b (Th2 activation, immunoregulation), and M2c (immunoregulation, matrix deposition, tissue remodeling) activation [17]. Although the definition of these pure activation states is based on defined in vitro conditions, and macrophage activation is likely to be much more diverse and complex in the in vivo setting, several studies have addressed the expression of polarization marker genes in GAMs, either in vitro or in vivo [9,13,18,19]. Similar to solid tumors in other organs, GAMs produce factors associated with an alternative macrophage activation, such as increased production of anti-inflammatory molecules (e.g. TGF-β1, ARG1, and IL-10), and molecules supporting tissue remodeling and angiogenesis (e.g. VEGF, MMP2, MMP9, and MT1-MMP). However, GAMs also produce pro-inflammatory molecules (e.g. TNF-α, IL1-β, and CXCL10) [9,10,11,18,20,21]. To date no comprehensive comparison of GAMs expression pattern with M1, M2a,b,c macrophages has been performed. In this study, we isolated CD11b + cells from GL261 mouse gliomas and naïve control mice and performed a genome-wide microarray expression analysis. We compared the GAMs expression profile to datasets of M1, M2a, M2b, and M2c in vitro-stimulated macrophages. Using flow cytometry to distinguish between resident microglia and invading macrophages/monocytes, we investigated the expression of selected genes in resident and invading GAMs in two different glioma mouse models by qRT-PCR. To show the relevance for potential therapeutic applications, we investigated some key genes in CD11b + GAMs isolated from human GBM samples compared to cells isolated from lower grade brain tumors, control microglia, and CD11bcells in GBM. Last, we used the The Cancer Genome Atlas (TCGA) database for glioma-patient survival data linked to gene expression data.

Ethics statement
This study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. This study was approved by the local ethics committees for animal experiments: the Institutional Animal Care and Use Committees of Cleveland Clinic, Lerner Research Instituteapproved protocol number 2013-1029 (LRI, last approved June 25 2013) and the Committee on the Ethics of Animal Experiments in Berlin (LaGeSo, Permit Numbers: GO 268-10, GO 343-10, GO 438-12). All surgery was performed under anesthesia, and all efforts were made to minimize suffering.
Freshly resected patient samples were provided by the Department for Neurosurgery Charité University Hospital and the Helios Clinics (both Berlin, Germany). Handling and analysis of these tissues was performed according to the rules and with the approval of the Ethical Committee and with patient's written consents (Ethics Committee: "Ethikkommission der Charité-Universitätsmedizin Berlin", application number: EA4/098/11).

Animals
All in vivo work with GL261 glioma cells was done in C57BL/6 wild-type mice (Charles River Laboratories, Wilmington, MA, USA). Age matched naïve C57BL/6 wild-type mice were used as controls.
To broaden the relevance of our findings, we employed another murine model where the tumor is initiated by the overexpression of PDGFb in Nestin-expressing cells in vivo: Ntv-a/ Ink4a-Arf -/mice develop pro-neural high-grade gliomas 6 to 8 weeks following intracranial injection of RCAS-PDGFb-producing DF-1 chicken fibroblast cells at 4.5 to 10 weeks of age [22]. The genetic backgrounds of tv-a mice are FVB/N, C57BL/6, BALB/C, and tv-a. RCAS-PDGFb tumors were dissociated as described below and intracranially re-transplanted into Cx3cr1 GFP/wt Ccr2 RFP/wt mice (kind gift of Richard Ransohoff, [23]) to distinguish between microglia and peripheral monocytes that invaded the brain. Naïve age-matched Cx3cr1 GFP/wt Ccr2 RFP/wt mice were used as controls. 5 to 7 weeks after intracranial injection of RCAS tumor cells into Cx3cr1 GFP/wt Ccr2 RFP/wt mice these mice develop tumors that are histologically identical to the original tumors in Ntv-a/Ink4a-Arf -/mice.

Human tissue
Tumor and control tissue was taken during surgery while patients were under a general anesthetic, and was placed immediately in culture medium for transportation. Cells were isolated from the tissue as described below no later than 2 hours after resection. Cell sorting was performed via MACS isolation, using anti-CD11b microbeads (described below).

Cell Isolation
GL261-implanted mice were sacrificed 20 days post-injection, RCAS-PDGFb-implanted mice were sacrificed 4-5 weeks post-injection. Tumor-bearing and control mice were euthanized by i.p. injection of 200 μl pentobarbital-sodium (Narcoren, Pharmazeutischen Handelsgesellschaft) and perfused using a 0.9% NaCl solution. The brain and spleen were extracted and stored in ice-cold HBSS (Gibco-Invitrogen). For naïve mouse brains, the olfactory bulbs and the cerebellum were cut by a scalpel and discarded. The rest of the tissue was used for dissociation. In tumor-bearing mouse brains, only the visible tumor area around the injection site was used.
Human and mouse tissue was dissociated with the Neural Tissue Dissociation Kit (Miltenyi Biotec, Bergisch-Gladbach, Germany) according to the manufacturer's instructions. To remove the myelin we followed a protocol published elsewhere [24]. In brief, the brain cell suspension was mixed with a total of 25 ml of a 22% Percoll (Th.Geyer, Renningen, Germany) solution and a layer of 5 ml cold PBS (Gibco-Invitrogen) was added on top. Centrifugation at 950 g with slow acceleration and without breaks created a gradient that separated the cell pellet on the bottom of the tube from the myelin which was carefully aspirated. For the isolation of GAMs from RCAS-PDGFb tumors-bearing and corresponding naïve Cx3cr1 GFP/wt Ccr2 RFP/wt mice a 30%/70% Percoll gradient was used. After 25 min of centrifugation at 800 g GAMs and naïve microglia were enriched at the 30%/70% interphase. Cells were collected, washed once with PBS, and subsequently centrifuged again at 300 g for 10 min. Subsequently, the cell pellet was resuspended in sorting buffer for subsequent magnetic-activated cell sorting (MACS) or fluorescence-activated cell sorting (FACS) isolation.
Spleens were processed through a 70 μm cell strainer with a syringe plunger and the mesh rinsed with 10 ml of PBS per spleen. The cells were centrifuged and the pellet subjected to erythrocyte lysis by adding 5 ml of 1x RBC lysis buffer (Cat# 420301, Biolegend, San Diego, CA, USA). The lysis was carried out by shaking the tube mildly for 5 min at RT and subsequently stopped with 20 ml of PBS. The pellet was washed once with PBS and resuspended in PBS, containing 0.5% FCS and 2mM EDTA (FACS buffer) for subsequent FACS isolation.
For RCAS-PDGFb tumors blood monocytes were used as peripheral controls (instead of monocytes isolated from spleens). To obtain blood monocytes, 200 μl of blood was mixed with 500 μl of PBS, containing 2.5 mM EDTA, centrifuged, the clear phase aspirated and the remaining phase mixed with 1 ml 1x RBC lysis buffer and incubated for 3 min on ice. The reaction was stopped by adding PBS, cells were centrifuged, washed once in 1x RBC lysis buffer and resuspended in FACS buffer.

MACS sorting
The CD11b + samples for the microarray were generated using magnetic activated cell sorting (MACS). Following percoll gradient centrifugation, tumor and control cell pellets were resuspended in PBS, containing 0.5% FCS and 2 mM EDTA and labeled with anti-CD11b microbeads (Miltenyi Biotec, 130-093-634). The MACS isolation was carried out according to the manufacturer's instructions and cells were subsequently used for RNA isolation.

RNA Protocols
Total RNA from MACS or FACS-sorted cells was isolated using the RNeasy Plus Mini Kit (Qiagen, Hilden, Germany). On-column DNase 1 (Qiagen) digestion was performed and total RNA was eluted in RNase-free water. RNA yield was measured using a Nanodrop 1000 (Nanodrop, Wilmington, DE, USA) spectrophotometer and quality was assessed using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Samples were stored at -80°C until further use. For qRT-PCR, first strand cDNA synthesis of RNA was done using the Superscript II (Invitrogen) reverse transcriptase according to the manufacturer's instructions. For mRNA transcription oligo-dT primers (Invitrogen) were used. cDNA was stored at -20°C until further processing.
Sample preparation for microarray hybridization was carried out as described in the Ambion WT Expression Kit Protocol (Life Technologies, Carlsbad, CA, USA) and the Affymetrix WT Terminal Labeling and Hybridization User Manual (Affymetrix, Inc., Santa Clara, CA, USA). In brief, 300 ng of total RNA were used to generate double-stranded cDNA. 12 μg of subsequently synthesized cRNA was purified and reverse transcribed into sense-strand (ss) cDNA, whereat unnatural dUTP residues were incorporated. Purified ss cDNA was fragmented using a combination of uracil DNA glycosylase (UDG) and apurinic/apyrimidinic endonuclease 1 (APE 1) followed by a terminal labeling with biotin. 3,8 μg fragmented and labeled ss cDNA were hybridized to Affymetrix Mouse Gene 1.0 ST arrays for 16 h at 45°C in a rotating chamber. Hybridized arrays were washed and stained in an Affymetrix Fluidics Station FS450, and the fluorescent signals were measured with an Affymetrix GeneChip Scanner 3000 7G.
Sample processing was performed at an Affymetrix Service Provider and Core Facility, "KFB-Center of Excellence for Fluorescent Bioanalytics" (Regensburg, Germany; www.kfbregensburg.de). Data sets are available at http://www.ebi.ac.uk/arrayexpress under the accession number E-MTAB-2660.

Bioinformatic Analysis
Two additional data sets consisting of three samples each that were generated from CD11b + MACS-isolated peritoneal myeloid cells were downloaded from http://www.ebi.ac.uk/ arrayexpress [25]-E-MEXP-3623; [26]-E-GEOD-25585) and used as peripheral controls. Two independent approaches were used for bioinformatic data analysis. For the first method the Affymetrix Expression Console Software Version 1.0 was used to create summarized expression values (CHP-files) from the expression array feature intensities (CEL-files). The Robust Multichip Analysis (RMA) algorithm was applied [27]. Integrative analysis of genomewide expression activities from DNA-microarray datasets was performed with the Gene Expression Dynamics Inspector (GEDI), a Matlab (Mathworks, Natick, MA) freeware program which uses self-organizing maps (SOMs) to translate high-dimensional data into a 2D mosaic [28]. Each tile of the mosaic represents an individual SOM cluster and is color-coded to represent high or low expression of the cluster's genes, thus identifying the underlying pattern.
For the second approach basal expression values (from CEL-files) were processed in R using Bioconductor package Affy [29]. The Affy-scale-value-expression-set-function was set to 500 and data was normalized with expresso-function, set to quantile normalization, and pmonly correcte.method and medianpolish summary method. Gene annotation was done with mogen-e10sttranscriptcluster [30], and only annotated probes were included in the analysis. Subsequently, collapseRows, was used to obtain a single representative probe for each gene, resulting in 13,943 genes taken into the analysis [31]. Signed network analysis was done with Weighted Gene Coexpression Network Analysis (WGCNA [32]). An adjacency network was made using a softpower (β) value of 20, based upon the scale free criteria, and low connected genes were filtered out of the analysis, resulting in 10,875 genes. Dendrogram formation and module determination was done by average linkage clustering and an arbitrary cut-off, determined by WGCNA. From each module, the Module Eigengene (ME) was calculated, which is the first principal component and functions as a representative of the expression profile of the module. Next, as a measure for intramodular connectivity, the ME was correlated to the expression profile of all intramodular genes, resulting in a kME-table, to determine which genes are most important, e.g. hub-genes. These module Membership values were multiple testing corrected using R package Stats. Genes significantly associated to these modules (FDR p value <0.01) were used for further analysis.
Heat maps were generated using heatmap.2 function of Bioconductor Package Gplots [33]. Functional annotation and transcription factor binding site enrichment analysis of the modules was done using Webgestalt [34,35]. The WGCNA UserlistEnrichment method was used to determine the significance of the overlap of several gene lists with the induced gene lists using a Hypergeometric test [31].
For the Gene Set Enrichment Analysis (GSEA) the 438 upregulated genes from our GAMs dataset were taken to form a gene set, and it was tested if this gene set is overrepresented in the upregulated or downregulated genes in the M1 or M2a,b,c macrophage data sets in comparison to M0. For each of M1, M2a, M2b, and M2c polarization, in total 4 runs of GSEA were executed. Taking M1 as an example, we fed the gene expression values of M1 and M0, each with 3 replicates, to the GSEA software; signal-to-noise ratios of each gene were calculated to rank the genes in a descending manner. An enrichment score (ES) for the gene set of 438 genes was then computed based on the gene ranking, which corresponds to the magnitude of the overrepresentation. Thereafter, the significance level of the ES was evaluated by a gene-set permutation test with 1000 permutations.

qRT-PCR
Quantitative RT-PCR reactions were performed using the SYBR Select Mastermix (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's instructions on a 7500 Fast Real-Time PCR System (Applied Biosystems).
Primers for target genes were designed to recognize all validated mRNA splice variants of the gene, relying on the RefSeq sequences of the UCSC genome browser (http://genome.ucsc. edu/). The primer sequences were generated using the Primer-Blast tool (http://www.ncbi.nlm. nih.gov/tools/primer-blast/). Amplification of unintended targets was excluded by BLAST search (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Calculation of unfavorable secondary structures and primer-dimer amplification was done with the IDT oligo analyzer (http://eu.idtdna. com/analyzer/applications/oligoanalyzer/). For absolute quantification of target genes, nested primers flanking the original primers were designed. Nested amplicons were purified from agarose gels, and serial dilutions were used to generate standard curves for each gene. Primer sequences can be found in S1 Table. Survival outcome analysis Cbioportal (http://www.cbioportal.org/public-portal/) was used to access patient survival information and tumor gene expression data of GBM patients (Study: Glioblastoma multiforme, TCGA Provisional, mRNA Expression z-Scores (microarray), accessed in August 2014) [36,37]. Subtype information was retrieved from [38]. The standard deviation of the gene expression values of all patients was calculated. Patients with a gene expression lower than the negative standard deviation were clustered into the low expression group, whereas patients with a gene expression higher than the positive standard deviation were clustered into the high expression group. For subtype analysis the patients were clustered according to the standard deviation of gene expression within that subtype. Survival plots were generated and Kaplan-Meier curves were produced and statistical significances between high and low expression groups were calculated with the Log Rank (Mantel-Cox) test using GraphPad Prism 5. To test for a proportional risk increase with each unit of gene expression, we fitted a Cox proportional hazards regression model to our data including a test for the proportional hazards assumption [39]. Analysis was done using the statistical script language R [40].

Gene expression profile of glioma-associated microglia/macrophages (GAMs)
To identify glioma-regulated transcripts in glioma-associated microglia/macrophages, we isolated microglia/macrophages from GL261 glioma-bearing brains using MACS sorting with anti-CD11b antibodies. As controls we isolated microglia from age-matched healthy mouse brain and performed an Affymetrix Genechip Mouse Gene 1.0ST microarray on these samples. Additionally, we also used two data sets that were generated from CD11b + MACS-isolated peritoneal myeloid cells as peripheral controls. These data sets were downloaded from http:// www.ebi.ac.uk/arrayexpress [25]-E-MEXP-3623; [26]-E-GEOD-25585). The two additional data sets consisted of three samples each, generated from either CD11b + MACS-sorted [26] or CD11b + MHC -II hi B220 -Gr1flow-sorted cells [25]. When comparing the GAMs data set to the three naïve control data sets, 783 genes were significantly upregulated at least 2-fold and 198 genes downregulated at least 2-fold (S2 Table).
We first applied the Gene Expression Dynamics Inspector (GEDI) on these datasets to visualize the global patterns of gene expression in GAMs versus naive microglia and macrophages. GEDI uses self-organizing maps to capture genome-wide transcriptome activity via 'gestalt' recognition [28]. GEDI facilitates the identification of genome-wide patterns with each mosaic tile in the map representing a gene cluster that is expressed at similar levels. The four GEDI maps, with blue colour indicating low and red colour high mRNA expression levels, show that the gene expression patterns of GAMs and naive microglia are more similar to each other than both macrophage datasets. Closer inspection of both microglia GEDI patterns then revealed a central cluster of highly expressed genes in GAMs that is different from naive microglia (white circle) (Fig. 1).

Weighted Gene Coexpression Network Analysis of the data sets
Using Weighted Gene Coexpression Network Analysis (WGCNA), we clustered genes into different modules, according to their co-expression pattern (Fig. 2). Within each module genes were ranked according to how closely they correlated to the Eigengene expression of this module (termed Module Membership). The Module Eigengene is the first principal component, representing the gene expression patterns of all genes clustered into this module [41]. Two modules contained genes that were upregulated in glioma-associated samples compared to the three control data set (labeled red and brown in Fig. 2). Table 1 shows the genes with the strongest upregulation in GAMs and the modules they were clustered into. S3 Table lists all 10,875 genes and the corresponding modules they were clustered into.
In addition, the gene clustering revealed genes that were specific for either microglia or macrophages. The orange module contains genes that were highly expressed in GAMs and already showed high expression in naïve microglia, but low expression in naïve macrophages (Fig. 2). Genes that showed high expression in GAMs and naïve macrophages, but low expression in naïve microglia were clustered into the green module (Fig. 2). Two studies used RNA sequencing and microarray analysis to determine microglia and macrophage-specific genes [42,43]. Genes that were found to be microglia-specific in these reports predominantly clustered into our orange module, whereas macrophage-specific genes are enriched in our green module, validating our findings.
We used Webgestalt to perform gene ontology (GO) enrichment analysis on all >2 fold upregulated genes that clustered into the glioma-regulated (red and brown) modules (438 genes in total) to identify overrepresented GO terms ( Table 2). Several GO terms were overrepresented in the GAMs data set that can be grouped into three main groups: regulation of immune response/activation, programmed cell death, and response to other organism/to virus. All genes included in the overrepresented GO terms are listed in S4 Table. Next, we used Webgestalt to identify transcription factor binding sites that are enriched in the same set of  (Table 3). This analysis revealed that enriched binding sites included known sites for IRF1, IRF2, IRF7, IRF8, IRF9/STAT1/STAT2, STAT5A, STAT5B, and NFAT. Both Stat1 and Stat2 clustered into the glioma-regulated brown module and were around 2.5 fold upregulated. In contrast, Stat5a and Stat5b were not significantly regulated at the transcriptional level. Irf7 and Irf9 were 6.9-fold and 1.4-fold upregulated and clustered into the brown and red module, respectively.

GAMs transcriptome only partially resembles an M1 or M2 polarization
To further analyze the phenotype of GAMs, we compared our glioma-associated data set with typical markers for an M1 and M2 macrophage phenotype (Table 4). Furthermore, we compared our data set with data sets of macrophages that were polarized in vitro into an M1 or M2a,b,c phenotype. We used data sets (Data set: E-GEOD-32690; [44]) that contained data of macrophages stimulated in vitro for 24 h with LPS and IFN-γ (M1 polarization), IL4 (M2a polarization), IFN-γ and complexed Ig (M2b polarization), and Dexamethasone (M2c polarization)in comparison to unstimulated M0 macrophages.
For this analysis we considered genes in our GAMs dataset that were glioma-regulated (genes that clustered into the red and brown module) and were >2-fold upregulated and WGCNA gene clustering reveals glioma-regulated gene modules. Each color represents a different module and each module contains genes with similar expression patterns over all four sample sets. The glioma-regulated modules (labeled as red and brown) were further analyzed, as they contained genes that were upregulated in the glioma-associated set when compared to all three control sets. A third module (labeled as green module) contained genes that were highly expressed in GAMs when compared to naïve microglia, but were already highly expressed in peripheral macrophages. The module that is depicted as orange contains genes that were highly expressed in GAMs, and also highly expressed in naïve microglia, however expressed at lower levels in peripheral macrophages. All other clusters contained genes that did not seem to be glioma-regulated. M2a (23 out of 438 genes, 5.3%, p value: 6.88e-09) was less strong (Fig. 3A). Table 5 lists all genes that are significantly enriched in our GAMs data set and are specific for either M1, M2a, M2b, or M2c polarization. Using Gene Set Enrichment Analysis (GSEA) we found that the 438 upregulated genes in GAMs are significantly enriched in the upregulated genes of all four macrophage polarization sets (Fig. 3B). Furthermore we tested the enrichment of the >2 fold upand downregulated M1, M2a,b,c genes in the entire GAMs data set that resulted from the WGCNA analysis (10,875 genes) (S1 Fig.). Six out of eight sets were significantly enriched in the up-and downregulated genes of the GAMs data set. However, 59.6% of the genes that were upregulated in GAMs (261 out of 438 genes) were not upregulated in any of the four macrophage phenotypes (Fig. 3C). This indicates that the GAMs phenotype shows only partial overlap with the classical M1 or M2 macrophage phenotype.

Validation of M1 and M2a,b,c marker gene expression in GAMs
In addition we selected five genes that were upregulated in GAMs, clustered into the gliomaregulated (red and brown) modules, and were specifically upregulated in either M1, M2a, M2b, or M2c-polarized macrophages and validated the expression of these genes by qRT-PCR in FACS-sorted samples. We selected Il1rn, Isg20 (both specific for M1-polarized macrophages), Clec7a (M2a), Tgfbi (M2b), and Cxcr4 (M2c). We FACS-isolated GAMs from GL261 tumors and RCAS-PDGFb tumors and measured the expression by qRT-PCR in comparison to microglia and peripheral macrophages isolated from control animals. For GAMs isolated from the GL261 mouse model, we distinguished between resident microglia (CD11b + /CD45 low ) and invading macrophages/monocytes (CD11b + / CD45 high /Ly6G -/Ly6C high ) by FACS sorting (Fig. 4A; [45]). In addition to microglia derived from naïve animals we also used spleen-derived macrophages/monocytes as additional controls.
As a further validation we used the RCAS-PDGFb tumor model, which produces pro-neural high-grade glioma, as a second glioma model system [22]. For this tumor model primary RCAS-PDGFb tumors were re-transplanted into Cx3cr1 GFP/wt Ccr2 RFP/wt mice, which allowed us to FACS-isolate RFP -/GFP + microglia and RFP + /GFP low macrophages/monocytes from tumors, as well as RFP -/GFP + microglia from control brains and RFP + /GFP low monocytes from peripheral blood without further antibody staining ( Fig. 4B; [23]).
In the GL261 model a significant higher expression of all selected genes could be confirmed in glioma-associated macrophages/monocytes and for Il1rn, Clec7a, and Cxcr4 in glioma-associated microglia, compared to the respective control cells (Fig. 5). In the RCAS model we saw a significant higher expression of all selected genes when comparing glioma-associated microglia to naïve microglia. When comparing glioma-associated macrophages/monocytes to naïve monocytes, only the expression of Il1rn, Tgfbi, and Cxcr4 was significantly higher. For Isg20 we detected a lower expression in glioma-associated macrophages/monocytes when compared to naïve monocytes. Furthermore, the expression of Clec7a was unchanged in glioma-associated macrophages/monocytes when compared to naïve monocytes. All investigated genes were expressed at higher levels in naïve monocytes when compared to naïve microglia. In part the expression was as high as in glioma-associated microglia. The difference of gene expression between the RCAS and the GL261 model might be partially caused by the different isolation procedures. We used GFP/RFP expression in Cx3cr1 GFP/wt Ccr2 RFP/wt mice for the RCAS model, as compared to antibody staining for CD11b/CD45/Ly6G/Ly6C in the GL261 model. Furthermore, the expression of the selected genes was lower in the RCAS model when compared to the GL261 modelthis was most notable for Isg20, Il1rn, and Clec7a. This might represent differences in the tumor biology of both models, but might also be owed to the different growth pattern and kinetics of the tumors. Mice injected with GL261 tumors were sacrificed 20 days post-injection, whereas mice implanted with RCAS tumors were sacrificed 4-5 weeks post-operation. Thus, GAMs in RCAS tumors were exposed to the tumor environment over longer times and the transcriptional program might have changed over the duration of the stimulus. Next, we collected human glioma samples and isolated CD11b + microglia/macrophages via MACS. We analyzed the gene expression of GAMs (CD11b + sorted cells) in GBM samples, the flow through of these tumor samples (the CD11b-negative fraction after MACS-isolation), control microglia samples (taken from brain resections of hippocampus, epilepsy, and trauma patients), and blood monocytes (Fig. 6).
We detected high expression of all genes in CD11b + GAMs isolated from GBM samples, however only the gene expression of IL1RN, ISG20 (both M1-specific), and TGFBI (M2b-specific) was significantly higher in CD11b + cells isolated from GBM samples, when compared to CD11b + cells isolated from control brain samples. The expression of CXCR4 was not significantly higher, whereas the expression of CLEC7A was significantly lower in CD11b + cells isolated from GBM samples. This might be in part due to the fact that the control brain specimen  were not taken from healthy patients, but patients suffering from epilepsy (3 samples) and trauma injury (1 sample). In these conditions the microglia might have already been polarized toward an M2-like phenotype.  Flow cytometry isolation strategies for generation of mouse qPCR samples. A) Plots depicting the strategy for isolation of microglia and macrophages from naïve brains and GL261 glioma-bearing brains using antibody staining for CD11b, CD45, Ly6C, and Ly6G. CD45 staining was used to distinguish between CD11b + /CD45 low resident microglia (gate P9) and CD11b + /CD45 high /Ly6G -/Ly6C high invading macrophages/monocytes (gates P8, P10, and P11), which were mostly absent in naïve brain samples. B) GAMs from RCAS-PDGFb tumors were isolated relying on an antibody-independent approach. Primary tumors from Ntv-a/Ink4a-Arf -/mice were reimplanted into Cx3cr1 GFP/wt Ccr2 RFP/wt mice and GAMS were sorted according to GFP and RFP expression. In naïve Cx3cr1 GFP/wt Ccr2 RFP/wt mice only GFP + cells were present in the brain. In tumor-bearing mice, resident microglia were GFP + , whereas invading macrophages/monocytes were RFP + /GFP low or GFP + /RFP + . RFP + /GFP low naïve blood monocytes were used as peripheral controls.
GAMs express the pro-tumorigenic genes Gpnmb and Spp1 Two features of high grade glioma are the aggressive invasion into the brain parenchyma and the immune-suppressive environment which prevents tumor rejection. Previously, we and other groups reported the expression of several pro-tumorigenic genes expressed in GAMs. We could confirm the expression of most of these genes in our dataset and a selection of these genes can be found in Table 4.
To identify novel glioma-regulated genes expressed in microglia/macrophages that might play a role in tumor-progression, we screened our dataset for genes that have been reported to play a pro-tumorigenic role in peripheral tumors, but have not been reported in GAMs. The genes Gpnmb and Spp1 were two of the highest upregulated genes in our GAMs data set (Table 1), clustered into the glioma-regulated (red) module, and have been implicated in immune-suppression (Gpnmb) or tumor cell invasion (Spp1) in peripheral tumors. To investigate the expression of these genes in GAMs via qRT-PCR in FACS-sorted samples from GL261 and RCAS-PDGFb tumors in comparison to naïve control cells we used the same samples as described above for the validation of the M1 and M2a,b,c-specific genes.
We were able to validate the expression of both genes using qRT-PCR. Furthermore, the regulation of the genes was similar in both brain tumor models -GL261, as well as RCAS-PDGFb tumors (Fig. 7). This shows that the up-regulation of these genes was not specific for Fig 5. qPCR validation of selected M1 and M2a,b,c-specific genes in murine GAMs. We selected 5 genes that were upregulated in GAMs and specific for either M1 (Il1rn and Isg20), M2a (Clec7a), M2b (Tgfbi), or M2c polarization (Cxcr4) and investigated the expression of these using qRT-PCR. For this we isolated GAMs from GL261 and RCAS-PDGFb tumors using flow cytometry, in order to distinguish between resident microglia and invading macrophages/ monocytes and used microglia, and spleen-derived macrophages/monocytes from naïve mice as controls. CTR MG: naïve microglia, CTR Mph: naïve monocytes, Tum MG: glioma-associated microglia, Tum Mph: glioma-associated macrophages/monocytes, GL261: cultured GL261 cells, Tum Mix: cultured RCAS-PDGFb tumor cells. Bar graphs illustrate the absolute number of transcripts normalized to 100,000 transcripts of Actb (n = 4). Analysis was done by students t test. Error bars indicate the Standard Error of Mean (SEM). *, p<0.05; **, p<0.01; ***, p<0.001 doi:10.1371/journal.pone.0116644.g005 one tumor model, but was observed in two different independent models. In addition, the expression of both genes was higher in GAMs when compared to the tumor cells (except for Gpnmb in the GL261 model), indicating that GAMs may be the primary source. Both genes were differently expressed in glioma-associated microglia and macrophages, as invading macrophages were the major source for Gpnmb in the RCAS-PDGFb model, whereas resident microglia were the main source for Spp1 in both models. The expression level of Gpnmb in both, glioma-associated resident microglia and invading macrophages/monocytes, was similar in both glioma models. In contrast, the expression of Spp1 was stronger in resident microglia and invading macrophages/monocytes isolated from RCAS-PDGFb tumors, when compared to the GL261 model.
To investigate the difference of gene expression in GAMs in the GL261 and the RCAS-PDGFb model in a broader set of genes we investigated the expression of six additional highly upregulated genes that clustered into the red or brown module in both models (S2 Fig.). From these six genes two were higher expressed in GAMs derived from the RCAS model (CD300lf and Cd200r4), two genes were expressed at similar levels in GAMs in both tumor models (Trem1 and Sh2d1b1), and two genes were higher expressed in GAMs derived from the GL261 model (Uck2 and Creb5).  (IL1RN and ISG20), M2a (CLEC7A), M2b (TGFBI), and M2c-specific (CXCR4) genes in CD11b + and CD11bcells isolated from human GBM (CD11b + n = 13, CD11bn = 5), control brain (CD11b + n = 5, CD11bn = 2) and blood monocyte samples (n = 2). Bar graphs illustrate the absolute number of transcripts normalized to 100,000 transcripts of ACTB. Analysis was done by students t test. Error bars indicate the SEM. *, p<0.05; **, p<0.01; ***, p<0.001 doi:10.1371/journal.pone.0116644.g006 Expression of GPNMB and SPP1 is upregulated in human GBMassociated microglia/macrophages As described above, we collected human samples and isolated CD11b + microglia/macrophages via MACS and used qRT-PCR to determine the expression of GPNMB and SPP1 in these  GPNMB and SPP1 expression is upregulated in human GAMs. We determined the expression of the genes GPNMB and SPP1 in CD11b + and CD11bcells isolated from human GBM (CD11b + n = 15, CD11bn = 9), meningioma (CD11b + n = 5, CD11bn = 3), grade III anaplastic astrocytoma (CD11b + n = 2, CD11bn = 2), control brain (CD11b + n = 5, CD11bn = 2) and blood monocyte samples (n = 2). The expression of GPNMB and SPP1 was significantly higher in CD11b + cells isolated from GBMs compared to CD11b + cells isolated from control brain, benign meningioma samples, blood monocytes, and CD11bcells in GBM. Bar graphs illustrate the absolute number of transcripts normalized to 100,000 transcripts of ACTB. Analysis was done by students t test. Error bars indicate the SEM. *, p<0.05; **, p<0.01; ***, p<0.001 samples. In addition to CD11b + and CD11bcells isolated from GBM samples and naïve brain tissue we also investigated meningioma and anaplastic astrocytoma (grade III glioma) samples.
We found that the expression of GPNMB and SPP1 was significantly higher in human GAMs isolated from GBM samples, when compared to non-tumor-associated control microglia and blood monocytes (Fig. 8). Furthermore, the expression of these genes was significantly higher in GAMs than in the tumor flow through, indicating that GAMs are the predominant source for these transcripts (S3 Fig.). The source and rate of expression of GPNMB seems to be dependent on the grade of malignancy of the tumor. In GBM GAMs were the main source for the expression, showing a 5 times higher expression than the GBM flow through, which comprises tumor cells and other stromal cells. However, when looking at meningioma, a benign tumor of the meninges, the expression of GPNMB was higher in the CD11bcell fraction, compared to the CD11b + cell fraction. The overall expression of GPNMB was higher in GBM than in meningioma and grade III astrocytoma samples. The expression level of SPP1 was also dependent on tumor grade. Similar to GPNMB, the expression of SPP1 was higher in GAMs isolated from GBM specimen, compared to the GBM flow through. Whereas, in lower grade tumors the expression of SPP1 was generally lower, which might indicate a possible role for SPP1 in higher grade tumors.

High expression of GPNMB and SPP1 in human GBM tissues is associated with poorer survival outcome
We used the cBioPortal database (http://www.cbioportal.org/public-portal/), to access TCGA data which links gene expression data to patient data, to investigate the effect of GPNMB and SPP1 expression on patient prognosis [36,37]. We grouped patients into low and high expression (gene expression lower than the negative standard deviation or higher than the positive standard deviation, respectively) and determined the overall survival for these patients. High expression of each gene has a negative effect on patient prognosis. Median survival was 19.77 months (low GPNMB expression) vs. 12.92 months (high GPNMB expression) and 15.31 months (low SPP1 expression) vs. 8.82 months (high SPP1 expression) (Fig. 9). This dataset also included G-CIMP + tumors that are mostly proneural subtype tumors and generally have a better overall survival prognosis when compared to G-CIMPtumors. We excluded these G-CIMP + tumors from our analysis and reanalyzed the data. Most G-CIMP + tumors exhibited a low expression of both, GPNMB and SPP1. Accordingly, the overall survival prediction for tumors with low expression of either GPNMB or SPP1 was less favorable after removing G-CIMP + tumors from the analysis. Median overall survival for GBM with low SPP1 or GPNMB expression was 14.16 months (low SPP1 expression including G-CIMP + tumors) vs. 12.56 months (low SPP1 expression excluding G-CIMP + tumors) and 17.7 months (low GPNMB expression including G-CIMP + tumors) vs. 14.82 months (low GPNMB expression excluding G-CIMP + tumors). We also plotted the overall survival of patients with intermediate expression of each gene (S4 Fig.) In addition, we grouped patients into the four molecular subtypes (proneural incl. G-CIMP + tumors and w/o G-CIMP + tumors, neural, mesenchymal, and classical) to investigate the effect of high GPNMB and SPP1 expression in these subtypes. Low GPNMB expression is associated with the most positive effect on patient prognosis in the proneural subtype if G-CIMP + patients are included (p = 0.09, median survival of 20.66 months (low expression) High GPNMB and SPP1 expression is associated with worsened survival prognosis in human GBM patients. Data taken from the TCGA database, showing survival probability of glioma patients grouped according to high and low expression of our target genes GPNMB (A) and SPP1 (B). High expression of both genes has a negative effect on patient prognosis. Patients were in addition grouped into the four molecular subtypes (proneural, neural, mesenchymal, and classical). Low GPNMB expression seems to have the most severe effect on patient prognosis in the proneural subtype when including G-CIMP + tumors, but no significant effect in the other subtypes. Low SPP1 expression seems to have the highest effect on patient prognosis in the neural and classical subtypes. vs. 9.21 months (high expression)), but not in the other subtypes. Low SPP1 expression is associated with the most beneficial effect on patient prognosis in the neural (p = 0.044, median survival of 25.90 months vs. 5.65 months) and the classical subtype (p = 0.09, median survival of 15.31 months vs. 9.11 months). Furthermore, both genes are differently regulated within the four subtypes. Mesenchymal tumors have a higher probability of high GPNMB and SPP1 expression, whereas classical tumors have a higher probability of low expression of both genes. G-CIMP + tumors, which generally have a better survival prognosis, display the lowest expression of both genes, indicating a possible role of these genes in advanced malignancy.
To test for a proportional risk increase for survival with the gene expression level, we fitted a Cox proportional hazards regression model to our data including a test for the proportional hazards assumption. The expression level of both genes had a significant effect on survival of glioblastoma patients (all glioblastoma) when G-CIMP + cases are included, as well as in the proneural subtype when G-CIMP + cases are included. All data sets, except neural subtype tumors when tested for SPP1 expression, passed the test for the proportional hazards assumption (S5 Table).

Discussion
Glioma-associated microglia/macrophages actively support glioma growth by the release of factors that stimulate angiogenesis, invasion, or suppression of immunity [9,10,11,15]. In the present study we performed a genome-wide expression analysis of these cells isolated from an experimental glioma mouse model. We have identified genes which are specifically regulated in glioma-associated microglia/macrophages when compared to naïve microglia or peripheral macrophages. We selected some of these genes and validated their expression in two independent glioma mouse models as well as in human glioma samples. In addition we could show that the expression of these genes differs in microglia and macrophages/monocytes. Our results are in accordance with previous reports that investigated the expression of selected target genes in GAMs via qRT-PCR [9,10,11], and markers known for tumor-associated macrophages in peripheral tumors. Several genes that are implicated in angiogenesis (Vegfa, Hgf), suppression of immunity (Arg1, Tgfb3), and tumor invasion (Mmp2, Mmp14, Ctgf) were also highly expressed in our GAMs data set (Table 4).
Previous studies reported that GAMs express both markers of the M1 and M2 macrophage phenotype [9,18]. However, no extensive comparison of these phenotypes has been performed yet. By comparing our expression data with those generated from M1, M2a, M2b, and M2c-polarized macrophages, we show that the GAM phenotype shows only partial overlap with the M1, M2a, M2b, and M2c phenotype. To date the term "M2-like polarization" has often been used for describing the polarization of GAMs; our data indicate that GAMs do not fit into a classical M1 or M2 phenotype, but represent a unique phenotype. As an alternate explanation it might reflect heterogeneity among GAMs. One might speculate that dependent on the location in the tumor tissue, some GAMs are polarized toward an M1-like phenotype, whereas other GAMs possess a more M2-like phenotype and another population of GAMs is not polarized toward M1 or M2-like states. Therefore, single-cell sequencing of GAMs would help to better understand the activation status of GAMs.
GAMs have been shown to promote tumor growth, rather than inhibiting it, e.g. by secreting factors that support glioma invasion or immunosuppressive factors [9,11,13,15]. Therefore, these cells represent an attractive target for anti-glioma therapy, as modulation of their activation state might be useful to inhibit glioma progression [14,16]. Here we present GPNMB and SPP1 as possible new targets and could show that these genes are highly expressed in GAMs in different glioma mouse models, in human GBM, and that high expression of these genes is correlated with shorter glioma-patient survival. The expression of neither of these genes was previously linked to GAMs.
The expression of GPNMB has been reported in tumor cells for different cancers, including glioma [46,47]. Furthermore, GPNMB expression was detected in microglia of non-neoplastic rat brains and increased with inflammation [48]. The data from our GL261 glioma model suggests that GL261 cells express GPNMB at a very high level. However, in the RCAS-PDGFb glioma mouse model, as well as in human GBM samples GAMs were the predominant source for GPNMB expression in all tested paired samples (S3 Fig.). GPNMB (also called Osteoactivin) is a transmembrane protein, but is also localized in the phagosome and can also be secreted, and might have different functions in GAMs and in the tumor. Ripoll et al., have shown that GPNMB acts as a negative regulator of pro-inflammatory macrophage activation in RAW264.7 cells [49]. Thus, high GPNMB expression in GAMs could participate in the modulation of the pro-tumorigenic phenotype of GAMs. Furthermore, GPNMB has been shown to inhibit T cell activation via direct cell-cell interaction of antigen-presenting cells and T cells, and could thus contribute to the immunosuppressive milieu in gliomas [50,51,52]. Finally, anti-Gpnmb antibodies conjugated with a cytotoxic agent are under investigation for the treatment of malignant glioma, breast cancer, and cutaneous melanoma [53,54,55].
Furthermore, we found that GAMs highly express SPP1. SPP1 (also called Osteopontin) is a secreted protein that has been postulated to increase tumor cell invasion in vivo and migration in vitro, and was found to be highly expressed in different types of cancers, such as lung cancer, ovarian cancer, and also glioma [56,57,58,59,60,61]. Furthermore, SPP1 has recently been identified as a ligand for CD44. SPP1-CD44 interaction was shown to increase stemness of CD44expressing glioma-initiating cells [62]. Here we show that GAMs and not other cells of the tumor microenvironment are the predominant source for SPP1 expression in glioma.
Two studies have previously performed screens of freshly-isolated GAMs. We used these data sets to investigate the expression of Gpnmb, Spp1, Il1rn, Isg20, Cxcr4, Tgfbi, and Clec7a in these studies (S6 Table). Huang et al., performed a screen of mouse GFP + chimeric GL261-associated and naïve (also isolated from the brain) bone-marrow-derived myeloid cells. However, using this approach they did not target the brain-resident microglia population that is also present in the tumor (Data set: E-GEOD-38283) [63]. All of the investigated genes, except Cxcr4, were also upregulated in our screenhowever to a lesser degree. This was partially due to an already high expression in the control cells. Furthermore, Murat et al., performed a microarray experiment on a paired sample of human GAMs and whole tumor lysate from the same patient (Data set: GSE16119) [64]. All of the 7 genes were higher expressed in GAMs when compared to the whole tumor lysate. However, it should also be noted that the whole tumor lysate was not depleted of GAMs.
Taken together, our findings show that GAMs are polarized toward a phenotype that has only partial overlap with the M1 or M2a, M2b, and M2c phenotypes. Furthermore, we identified GAMs as the predominant source for the pro-tumorigenic proteins GPNMB and SPP1 in murine and human malignant gliomahighlighting the importance of macrophages and microglia as therapeutic targets in anti-tumor treatment regimens.   Table. Comparison with other screens on GAMs. The expression values of Gpnmb, Spp1, Il1rn, Isg20, Clec7a, Tgfbi, and Cxcr4 in our microarray screen and in two studies that performed microarrays on GAMs are listed [63,64]. (DOCX)