Single-cell characterization of macrophages in glioblastoma reveals MARCO as a mesenchymal pro-tumor marker
Genome Medicine volume 13, Article number: 88 (2021)
Macrophages are the most common infiltrating immune cells in gliomas and play a wide variety of pro-tumor and anti-tumor roles. However, the different subpopulations of macrophages and their effects on the tumor microenvironment remain poorly understood.
We combined new and previously published single-cell RNA-seq data from 98,015 single cells from a total of 66 gliomas to profile 19,331 individual macrophages.
Unsupervised clustering revealed a pro-tumor subpopulation of bone marrow-derived macrophages characterized by the scavenger receptor MARCO, which is almost exclusively found in IDH1-wild-type glioblastomas. Previous studies have implicated MARCO as an unfavorable marker in melanoma and non-small cell lung cancer; here, we find that bulk MARCO expression is associated with worse prognosis and mesenchymal subtype. Furthermore, MARCO expression is significantly altered over the course of treatment with anti-PD1 checkpoint inhibitors in a response-dependent manner, which we validate with immunofluorescence imaging.
These findings illustrate a novel macrophage subpopulation that drives tumor progression in glioblastomas and suggest potential therapeutic targets to prevent their recruitment.
Glioblastoma (GBM) is a devastating primary brain malignancy. Recurrence of GBM is inevitable despite the standard treatment of surgery, chemotherapy, and radiation—and the median survival is limited to around 16 months . Barriers to treatment include the complex interactions of macrophages in the tumor microenvironment, which play a variety of pro-tumor roles in gliomas [2,3,4]. While immunotherapies have been successful across a variety of other cancers , they are hindered in GBM by an immunosuppressive microenvironment including tumor-associated macrophages (TAMs) . Indeed, a recent study of checkpoint inhibitor therapy in GBM found an association between infiltration of HLA class II-deficient macrophages with tumor profiles unfavorable to immunotherapy . In addition, these TAMs are also hypothesized to contribute to the mesenchymal expression subtype, traditionally associated with poorer outcomes and treatment resistance .
Targeting these TAMs (such as with CSF1R inhibitors) is an intriguing therapeutic option but requires a better understanding of markers specific to and necessary for TAM functioning [9,10,11]. Our classical knowledge of macrophage polarization (M1 vs. M2) provides a simplified illustration of these different states . However, this is clearly not the whole picture, given that M1 and M2 genes are often co-expressed in individual TAMs . Furthermore, in GBMs, there is a diverse spectrum of monocytic lineage cells (which we broadly term as “macrophages”), comprising bone marrow-derived macrophages (BMDMs) recruited from the blood as well as tissue-resident microglia. Thus, the specific markers and pathways involved in pro-tumor macrophages in GBM still remain elusive.
Due to the diversity of cell types in the tumor microenvironment, bulk expression profiles are not ideal for investigating cellular subpopulations. Instead, single-cell RNA sequencing (scRNA-seq) has proven instrumental in understanding the heterogeneity within GBM [7, 13]. In recent years, there have been several published scRNA-seq datasets in gliomas [12, 14,15,16,17,18,19], which have helped elucidate the general differences between BMDM and microglial populations in the tumor microenvironment. However, each study has been limited in terms of patient numbers. Here, we combine 9 new scRNA-seq samples from glioblastomas with 57 previously published cases to explore the profiles of macrophages in gliomas at an unprecedented scale. In doing so, we identify a novel pro-tumor macrophage marker in GBM (MARCO) validated by immunofluorescence imaging. In addition to studying the clinical effects of MARCO expression on survival and during immunotherapy, we also characterize its associations with mesenchymal, hypoxic, and anti-inflammatory signatures.
Study design and data acquisition
Published single-cell RNA-seq count matrices were obtained from several repositories [20,21,22,23,24,25] and joined with 9 previously unpublished cases  to form a cohort of 50 GBM and 16 LGG samples (Additional file 1: Table S1). Bulk expression and survival data were obtained from 528 GBM cases from TCGA-GBM  and 75 GBM cases from Wang et al. . Expression subtyping and IDH1 mutation status for TCGA were acquired from Ceccarelli et al. . Gene sets were obtained from MSigDB v6.2 , with the exception of the BMDM versus microglia gene sets, which were obtained from Yuan et al. .
The GBM specimens were collected from surgical resections of de-identified patients at Columbia University Irving Medical Center who provided informed consent to participate in these studies through a protocol approved by the Columbia Institutional Review Board (IRB-AAAJ6163). Two different methods were used to dissociate the tissue specimens into single-cell suspensions. PDC001, PJ052, PJ053, PW032-706All, and PW032-710All were dissociated using the method previously published by Yuan et al. . PW039-705, PW035-710All, PW016-703_All, and PW017-703_All were dissociated into single cells using the Adult Brain Dissociation Kit (Miltenyi Biotec) on a gentleMACS Octo Dissociator with Heaters (Miltenyi Biotec) according to the manufacturer’s instructions. Single-cell suspensions were applied to an automated microwell-based platform for scRNA-seq library construction as previously described by Yuan et al. . scRNA-seq libraries for PDC001, PW039-705, PW035-710All, PJ052, PJ053, PW016-703_All, and PW017-703_All were sequenced on an Illumina NextSeq 500 with an 8-base index read, a 21-base read 1 containing cell-identifying barcodes (CBs) and unique molecular identifiers (UMIs), and a 63-base read 2 containing the transcript sequence. The raw sequencing data were processed as described by Yuan et al.  to generate the digital gene expression matrices. scRNA-seq libraries for PW032-706All and PW032-710All were pooled and sequenced on an Illumina NovaSeq 6000 with an 8-base index read, a 26-base read 1 containing CBs and UMIs, and a 91-base read 2 containing the transcript sequence. The raw sequencing data were first corrected for index swapping to avoid cross-talk between sample index sequences  and then aligned to generate the digital gene expression matrices as described by Szabo et al. .
All datasets were first filtered to remove mitochondrial and ribosomal proteins. Datasets were then merged together (separately for GBM and LGG) keeping the intersection of genes, with genes with zero total counts being discarded. Raw counts were then normalized to log2(1 + TPK), as described in Yuan et al. . The intersection of expressed genes with LM22 from CIBERSORT  was used as a filtered list to reduce the batch effect. For visualization, first, principal components analysis (PCA) was applied to reduce the total dimensionality to 5% of the number of genes then uniform manifold approximation and projection (UMAP)  with default parameters to non-linearly reduce that into a two-dimensional embedding. Three cell types were readily identified using the standard markers of CD14, CD3, and SOX2 for macrophages, T lymphocytes, and tumor cells , respectively. The macrophages were isolated, and the same dimensionality reduction procedure described above was used to generate an embedding for the macrophage population. The embedding of the macrophage population was separated through k-means clustering. Silhouette scores were calculated from k=2 to k=6, with k=2 providing the highest silhouette score. This filtered list was used to assess the top enriched genes characterizing each of these two macrophage clusters.
To validate the presence of MARCO in tissue from patients with GBM, we co-stained MARCO with CD163 using dual stain immunofluorescence (IF). Tissue specimens were evaluated from one patient with each tumor subtype, including IDH1-wild-type GBM, IDH1-mutant GBM, and grade III anaplastic astrocytoma (LGG). In addition, two patients treated with anti-PD1 inhibitors were compared pre- and post-anti-PD1, including a responder and a non-responder. IF staining was performed on 5-μm sections using anti-MARCO (1:100, PA5-64134, Invitrogen (Carlsbad, CA)) followed by biotinylated goat anti-rabbit IgG (1:200, BA1000, Vector Laboratories) and streptavidin-conjugated 594 (1:1000, S11127, Invitrogen), as well as anti-CD163 (1:100, 10D6, Biocare Medical (Pacheco, CA)) followed by biotinylated horse anti-mouse IgG (1:200, BA2000, Vector Laboratories) and streptavidin-conjugated 488 (1:300, Invitrogen). Stained sections were examined and photographed using a Nikon Eclipse TE2000-E. Composite images from emission wavelengths of 395, 488, and 568 were captured from each field of interest. All images in Fig. 4c were thresholded at a fixed value for the MARCO channel using ImageJ. For single-stain controls, all channels were thresholded at fixed values.
TCGA-GBM U133 microarray data and pre-treatment expression data from Wang et al.  were used to determine the expression of MARCO. Each cohort was first log- and Z-score normalized independently before being combined. These normalized values were merged with survival data, including overall and disease-free survival. Survival differences were assessed in two ways: (1) dichotomizing MARCO expression across the median into MARCO-high and MARCO-low then comparing the survival curves with the log-rank test and (2) generating a univariable Cox model directly based on the expression of MARCO and taking the Wald p-value of the MARCO covariate. In Fig. 2a and b, we only used cases that were known to be IDH1-wild-type.
Expression subtype score
The single-cell subtype score method from Patel et al.  was used to calculate the signature score of each Verhaak expression type . This same method was used to access the BMDM versus microglial signatures from Yuan et al.  as well as for assessing the combined score of the tumor-macrophage cross-talk genes.
Gene set enrichment analysis
GSEA was performed with the official client v4.0.1 for Linux [30, 36]. The normalized gene expression of all macrophages was used as the input data, with MARCO expression as the phenotype label, and Pearson correlation as the metric for ranking genes. Default weighting (p = 1) was used. Phenotype permutation with 1000 iterations was used to calculate p-values, and the final gene sets were ranked by enrichment score.
The mean normalized expression of MARCO within the macrophage population was calculated for each GBM sample, and compared to a selected list of recruitment factors  within the tumor population of those same samples: macrophage colony-stimulating factor (CSF1), granulocyte/macrophage colony-stimulating factor (CSF2), hepatocyte growth factor (HGF), monocyte chemotactic protein 1 (MCP-1), macrophage inhibitory factor (MIF), stromal-derived factor 1 (SDF-1), transforming growth factor β (TGF-β, including TGFB1, 2, and 3), interleukin 10 (IL-10), osteopontin (SPP1), and lactadherin (MFGE8). An expression subtype score was calculated using this gene set and averaged across the tumor cells from each case. This score was then compared to MARCO expression in the corresponding macrophage population via Spearman correlation. The same procedure was used to assess the correlation of MARCO expression in macrophages with expression subtype scores in paired tumor cells.
To quantify the utility of MARCO as a marker compared to conventional myeloid, BMDM, and microglia markers, we assessed the normalized dispersion coefficient of each gene. We used the highly_variable_genes function in SCANPY  with default parameters and batch correction. We determined p-values with the exact permutation test.
All statistical analyses were conducted in Python 3.6. In all boxplots, the center lines represent the median, lower and upper box limits are respectively the first and third quartiles, and whiskers represent the maximal values up to 1.5 times the interquartile range. All values extending beyond this range are considered fliers/outliers. Shaded regions in survival curves represent 95% confidence intervals. Violin plots use the Gaussian kernel to estimate densities, and the center lines represent the median. The two-sided Mann-Whitney U test was used throughout to non-parametrically compare two populations. p-values were adjusted for multiple comparisons using the Benjamini-Hochberg (FDR) procedure, and statistical significance was assessed at an adjusted p-value threshold of 0.05.
Single-cell identification of a MARCO+ subpopulation of macrophages in GBM
To understand the heterogeneity of cellular profiles in GBM, we collected single-cell RNA-seq data from nine GBM cases and combined it with 41 GBM cases from previously published studies for a total of 79,968 single-cell transcriptomes (Additional file 1: Table S1). To reduce the batch effect, we filtered our gene list to the 499 genes that overlapped with LM22, the reference matrix used by CIBERSORT  to differentiate immune cells (Additional file 2: Table S2). After filtering, log-normalization, batch effect reduction, and dimensionality reduction, visualization of these cells produced a distinct macrophage population of 17,132 cells characterized by CD14 (Fig. 1a, Additional file 3: Fig. S1). Unsupervised k-means clustering on this macrophage population revealed two subpopulations (validated by silhouette score, Additional file 3: Fig. S2). One of these subpopulations was enriched in inflammation-related genes, with CCL4 and IL1B as the top two genes on the filtered gene list (p = 0.004 and p = 0.008, respectively, exact permutation test, n = 499 genes; Fig. 1b). The other subpopulation, opposite to the inflammatory side, had MARCO (macrophage receptor with collagenous structure) as the top gene (p = 0.004, exact permutation test, n = 499 genes; Fig. 1b). MARCO remained the top enriched gene for this side even when left out during the dimensionality reduction, demonstrating its representativeness (p = 0.004, exact permutation test, n = 499 genes). Looking at the expression of MARCO in all cell populations, we found that it is specifically found in this subpopulation of macrophages (Additional file 3: Fig. S3A). MARCO is a scavenger receptor normally found on alveolar macrophages  with a variety of immunomodulatory roles [40,41,42]. Given the negative role of MARCO in other cancers [43,44,45], we chose to focus on these MARCO+ macrophages and their impact on the tumor microenvironment.
Absence of MARCO expression in LGG and IDH1-mutant GBM
We next investigated whether this macrophage subpopulation could also be observed in lower-grade gliomas (LGGs), which include grade II and III astrocytomas as well as oligodendrogliomas. We collected single-cell expression profiles of 18,047 cells from 16 LGGs and processed them in the same manner as for GBM above (Fig. 1c). Although 2199 macrophages were identified, only 1% (23 cells) had non-zero expression of MARCO (Fig. 1d) compared to the 12% in GBM (2092 cells out of 17,132, p < 0.001; chi-squared test). This is despite the GBM samples having lower library complexity (p < 0.001, n = 49 patients, Mann-Whitney U test; Additional file 3: Fig. S4). Since most LGGs are IDH1-mutated, we examined MARCO expression in the 4 GBM cases with IDH1 mutations. Of the 281 macrophages from IDH1-mutated GBMs, only one cell had non-zero expression of MARCO. Comparing normalized expression levels, the mean MARCO expression in IDH1-wild-type GBM macrophages was 160 times higher than in IDH1-mutated GBM (p = 0.020, n = 41 patients, Mann-Whitney U test) and 64 times higher than in LGGs (p = 0.0007, Fig. 1e top, n = 49 patients). Meanwhile, the mean expression of other macrophage markers such as CD14 was not significantly different (p = 0.34 and p = 0.30, Fig. 1e bottom, n = 41 and n = 49 patients). Immunofluorescence imaging in IDH1-wild-type cases validates that MARCO co-localizes with CD163 macrophages in tissues from these patients, while IDH1-mutated and LGG tissues show rare, if any, co-localization (Additional file 3: Figs. S5 and S6). Therefore, MARCO expression is almost exclusively found in macrophages of IDH1-wild-type GBM, rather than their less deadly IDH1-mutated or lower-grade counterparts.
MARCO bulk expression associates with poor clinical outcomes and mesenchymal subtype
To understand the clinical consequences of MARCO expression, we examined its bulk expression in two published datasets: TCGA-GBM and the GBM cohort from Wang et al.  (for a total of n = 603 patients). Splitting patients into the MARCO-high and MARCO-low groups down the median, we found a negative association of MARCO expression with both overall survival (OS: p = 0.0046, n = 592 patients, log-rank test, Additional file 3: Fig. S7A) and disease-free survival (DFS: p = 0.018, n = 387 patients, Additional file 3: Fig. S7B). However, given the aforementioned association of MARCO with IDH1 and the strong impact of IDH1 mutations on survival, we repeated the analysis with only the 437 GBMs with known IDH1-wild-type status, and the effect persisted (OS: p = 0.0084, n = 437 patients; DFS: p = 0.035, n = 267 patients, log-rank test, Fig. 2a, b). Continuing this analysis with the IDH1-wild-type cohort, we treated MARCO expression as a continuous variable and found that the effect of MARCO in a univariable Cox model was also significantly detrimental to survival (p = 0.00057 for OS, n = 437 patients, hazard ratio = 1.19; p = 0.0047 for DFS, n = 267 patients, hazard ratio 1.19; Wald test). These effects were most pronounced in terms of long-term survival, with MARCO-high patients having roughly half the overall survival rate at the 2-year (18.0% vs. 31.8%) and 5-year (3.7% vs. 7.8%) time points compared to MARCO-low patients. Similarly, 2-year DFS was over three times lower in the MARCO-high population (5.8% vs. 18.1%; no data for 5-year DFS). The effect of MARCO expression on survival remained significant after controlling for common confounders including age and MGMT methylation status (Additional file 3: Fig. S7C).
Recapitulating the single-cell results, we found an association of IDH1 mutations with decreased MARCO expression in bulk data (p = 0.0039, n = 482 patients, Mann-Whitney U test, Fig. 2c). This was also significant with dichotomizing the population into MARCO-high and MARCO-low patients (p = 0.0016, n = 482 patients, Fisher Exact test, Fig. 2d). We then compared MARCO bulk expression with transcriptomic subtype . Although MARCO is not on the original list of mesenchymal genes , its expression was highly enriched in mesenchymal samples (p < 0.00001 for all pairwise comparisons between mesenchymal and other subtypes, n = 565 patients; Fig. 2e). These bulk transcriptomic results support our single-cell findings of the enrichment of MARCO in IDH1-wild-type tumors and also demonstrate an association of MARCO expression with worse prognosis and the unfavorable mesenchymal subtype.
Single-cell association of MARCO with mesenchymal traits and hypoxia
To uncover the source of this mesenchymal signature, we went back to our single-cell expression data. The mesenchymal expression signature was found primarily in the macrophage population rather than in tumor cells (p < 0.00001, n = 50 patients, Mann-Whitney U test, Fig. 3a, Additional file 3: Fig. S8, Additional file 2: Table S3) and specifically within MARCO-expressing macrophages (p = 0.005). MARCO mean expression in macrophage cells also significantly correlated with the mesenchymal signature score in tumor cells from the same sample (p = 0.0084, n = 41 patients, Spearman correlation; see the “Methods” section; Fig. 3b), while anticorrelating with the proneural signature (p = 0.047). We then performed gene set enrichment analysis (GSEA) based on the Pearson correlation of all other genes with MARCO expression to understand what processes underlie this subpopulation. Out of the 50 hallmark gene sets from MSigDB, the gene sets with the highest enrichment scores among MARCO+ macrophages were epithelial-mesenchymal transition, angiogenesis, glycolysis, and hypoxia (FDR q < 0.001 for all four, permutation test on phenotypes with default weighting, n = 17,132 cells; Fig. 3c, d). This association with hypoxia is supported by single-cell data from Darmanis et al. , in which samples were taken from the core and periphery of the same four tumors. Macrophages from the tumor core had higher MARCO expression than those taken from the tumor periphery (p < 0.001, n = 1846 cells; Mann-Whitney U test; Fig. 3e). Furthermore, we investigated the expression of MARCO in Ivy GAP , a database with laser-microdissected specimens from different anatomic structures. We found significantly different expression of MARCO in various structures (p = 0.014, n = 270 specimens from 37 patients; Kruskal-Wallis test, Additional file 3: Fig. S9), with the highest expression in the perinecrotic zone within the cellular tumor. These spatially informed findings support the GSEA characterization of MARCO+ macrophages as residing in the hypoxic tumor core.
MARCO+ macrophages demonstrate loss of inflammatory pathways and antigen presentation
Since hypoxia can polarize macrophages toward a pro-tumor phenotype , we also investigated if appropriate antitumoral responses were downregulated within MARCO+ macrophages. Also via GSEA, we found that the gene sets with the lowest enrichment scores among the hallmark set were all pro-inflammatory: interferon alpha response, interferon gamma response, allograft rejection, and TNFa signaling via NFKB (FDR q ≤ 0.001 for all four sets, n = 17,132 cells, Fig. 3f). Interestingly, among the eight individual genes with the lowest relative expression in MARCO+ macrophages, four of them were HLA class II genes (HLA-DRB1, DRA, DPA1, and DPB1), and another was CD74 (MHC class II invariant chain). Accordingly, the GO gene set antigen processing and presentation of peptide or polysaccharide antigen via MHC class II was highly downregulated in MARCO+ macrophages (p < 0.001, n = 17,132 cells, Fig. 3g). Similarly, by comparing the differential expression of HLA class II genes in relation to all 17,496 other assayed genes, MARCO was one of the three most negatively associated (Fig. 3h) alongside BNIP3 (an apoptotic Bcl-2 family gene) and VCAN (an extracellular protein implicated in metastasis ). Loss of HLA class II expression in macrophages is generally associated with an anti-inflammatory, inactivated state  and has been previously linked to worse outcomes in melanoma  as well as unfavorable tumor profiles in the context of PD1 immunotherapy in GBM .
Dynamics of MARCO expression under PD1 immunotherapy
To determine if the anti-inflammatory properties of MARCO+ macrophages play a role in PD1 checkpoint inhibitor therapies, we investigated MARCO expression in a longitudinal cohort of 17 PD1-treated GBM patients . Interestingly, there was a decrease in MARCO between pre- and post- immunotherapy recurrences (p = 0.02, n = 37 time points; Mann-Whitney U test; Additional file 3: Fig. S10A), but this was solely found within responders (p = 0.02, n = 19 time points from 10 patients; Fig. 4a). A timeline of two representative cases from this cohort shows a responder that had a strong decrease in MARCO expression following immunotherapy, whereas a non-responder had an increase in MARCO expression after immunotherapy, followed soon after by death (Fig. 4b). Immunofluorescence imaging of tissues from these same two cases visually recapitulates these bulk transcriptomic findings (Fig. 4c, Additional file 3: Fig. S11). While this supports the negative role of MARCO in the long-term following adjuvant therapy, the dynamics may differ in the short-term—in an orthogonal study of neoadjuvant PD1 therapy in GBM, MARCO was reported as one of the top increased genes in patients treated with pembrolizumab approximately 2 weeks prior to sample acquisition . Meanwhile, within a longitudinal cohort of 86 patients treated with standard therapy , we found no significant difference in MARCO expression before and after treatment (p = 0.21, n = 160 time points, Additional file 3: Fig. S10B), suggesting that these changes in MARCO are specific to immunotherapy.
Recruitment of MARCO+ macrophages from the blood by tumor cells
One important question for the potential targeting of MARCO+ macrophages is where they originate from and how they are recruited to the tumor. Based on the gene sets associated with BMDMs versus resident microglia , we found that MARCO+ macrophages more closely resembled BMDMs (p < 0.0001, n = 71 subpopulations; Mann-Whitney U test; Fig. 5a, b). Compared to non-MARCO-expressing macrophages, CD163 was among the highest differentially expressed genes in MARCO+ macrophages (p = 0.0037, n = 17,509 genes; exact permutation test), and TMEM119 was among the lowest (p = 0.0055), which are classic markers characterizing BMDMs and microglia, respectively. Nonetheless, MARCO appears to define a more specific subpopulation than solely BMDM-microglial differences (p = 0.004, n = 17,509 genes; see the “Methods” section; Additional file 3: Fig. S3B).
To understand the factors that recruited these macrophages, we examined the expression of a list of 10 genes known to attract and re-program macrophages in GBM . We found a significant positive correlation of MARCO mean expression in macrophage cells with the corresponding normalized expression of these recruitment factors in tumor cells from the same sample (p = 0.008, n = 41 patients, Spearman correlation; see the “Methods” section; Fig. 5c). These results are supported by correlations in bulk expression data from TCGA-GBM, in which MARCO is positively associated with this same signature of recruitment factors (p < 0.001, n = 528 patients, Spearman correlation; see the “Methods” section; Additional file 3: Fig. S12). Furthermore, we found within our GBM single-cell data that there is significant co-expression of Ki-67 with MARCO (p < 0.00001, n = 17,132 cells, chi-squared test; Additional file 3: Fig. S13), suggesting that these MARCO+ macrophages may also be proliferating within the tumor microenvironment.
In this manuscript, we have characterized the role of MARCO as a marker of pro-tumor macrophages in GBM. In particular, we have found that this MARCO-expressing subpopulation associates with mesenchymal, hypoxic, and anti-inflammatory traits as well as poor clinical prognosis. The mesenchymal nature of MARCO is supported by the role of MARCO in regulating the epithelial-mesenchymal transition outside the context of cancer . Meanwhile, the upregulation of glycolysis and hypoxia gene sets is consistent with the enrichment of MARCO+ macrophages in the tumor core—which is known to be associated with pro-tumor macrophages . Furthermore, the localization of these macrophages within the hypoxic, necrotic core is consistent with the observations that MARCO is upregulated in regions of ischemic brain, suggesting a functional role of this scavenger receptor in clearing cellular debris . Our finding of the opposition of MARCO to inflammatory genes and pathways is also seen in mice, where tolerized BMDMs upregulate MARCO . In fact, the remarkable downregulation of HLA class II genes on MARCO+ macrophages is consistent with murine studies where MARCO expression is associated with a decrease of antigen internalization capacity . Although MARCO has been previously reported as a pro-tumor TAM marker in non-small cell lung cancer , lung adenocarcinoma , and breast cancer , and has been associated with poor prognosis in periampullary adenocarcinoma , its role has not been previously reported in GBM. Our findings of the pro-tumor TAM traits associated with MARCO in GBM adds credence to its importance across cancers.
Here, we show that MARCO+ macrophages appear to be recruited from the blood via the upregulation of a set of factors secreted by tumor cells, including CSF1 and TGF-β. CSF1 expression in tumor cells has been previously shown to be related to higher proportions of TAMs in GBM . Notably, these TAMs were observed to express the cognate receptor CSF1R, which is targetable by existing therapeutics [11, 17]. Meanwhile, TGF-β has been experimentally shown to upregulate MARCO expression in M0 BMDMs , joining a host of other studies implicating TGF-β in glioma progression [28, 55]. While these recruitment factors are potential targets, MARCO itself also presents a promising target—anti-MARCO therapeutic antibodies have demonstrated efficacy in mouse melanoma models [44, 56]. As with other immunotherapies, we do not expect monotherapy with anti-MARCO antibodies to markedly improve outcomes alone. However, including anti-MARCO in combination regimens may be advantageous. While MARCO expression remains unaltered across standard therapy, we found that its expression changes in the course of anti-PD1 immunotherapy, with responders exhibiting decreases in MARCO in the long term following treatment. This result suggests that MARCO+ macrophages may be detrimental to checkpoint immunotherapy, which has been largely unsuccessful in GBM . Simultaneous targeting of MARCO or its associated recruitment factors may provide a promising opportunity to manipulate macrophage polarization and thereby boost the efficacy of checkpoint inhibitor therapy [58, 59]. Furthermore, there is a critical need for markers of response to anti-PD1 therapies. Previously, we had shown that mutations in MAPK genes and loss of PTEN predict response to anti-PD1 immunotherapy . However, there is no current marker of response after treatment with anti-PD1 inhibitors. MARCO may be useful to satisfy this need, and validation of this in prospective trials with anti-PD1 therapy could be beneficial.
In this work, we utilized single-cell RNA sequencing to identify a pro-tumor subpopulation of macrophages characterized by MARCO expression. These macrophages are associated with both worse clinical outcomes as well as with the mesenchymal subtype. Furthermore, these macrophages appear to be derived from the bone marrow and affected by anti-PD1 checkpoint inhibitors. These results highlight a novel macrophage subpopulation that contributes to tumor progression in glioblastomas and suggest potential therapeutic strategies for its mitigation.
Availability of data and materials
All novel sequencing data have been deposited in the Gene Expression Omnibus under accession number GSE141383: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE141383 . All previously published datasets can be found in their respective public repositories [20,21,22,23,24,25]. The processed and source data underlying each figure are provided in Additional file 4. All custom code is available at https://github.com/RabadanLab/MARCOsinglecell .
Bone marrow-derived macrophage
Gene set enrichment analysis
Single-cell RNA sequencing
Uniform manifold approximation and projection
Unique molecular identifier
Gilbert MR, Dignam JJ, Armstrong TS, Wefel JS, Blumenthal DT, Vogelbaum MA, et al. A randomized trial of bevacizumab for newly diagnosed glioblastoma. N Engl J Med. 2014;370(8):699–708. https://doi.org/10.1056/NEJMoa1308573.
Bowman RL, Klemm F, Akkari L, Pyonteck SM, Sevenich L, Quail DF, et al. Macrophage ontogeny underlies differences in tumor-specific education in brain malignancies. Cell Reports. 2016;17(9):2445–59. https://doi.org/10.1016/j.celrep.2016.10.052.
Chen Z, Feng X, Herting CJ, Garcia VA, Nie K, Pong WW, et al. Cellular and molecular identity of tumor-associated macrophages in glioblastoma. Cancer Research. 2017;77(9):2266–78. https://doi.org/10.1158/0008-5472.CAN-16-2310.
Hambardzumyan D, Gutmann DH, Kettenmann H. The role of microglia and macrophages in glioma maintenance and progression. Nature Neuroscience. 2016;19(1):20–7. https://doi.org/10.1038/nn.4185.
Ribas A, Wolchok JD. Cancer immunotherapy using checkpoint blockade. Science. 2018;359(6382):1350–5. https://doi.org/10.1126/science.aar4060.
Dutoit V, Migliorini D, Dietrich PY, Walker PR. Immunotherapy of malignant tumors in the brain: how different from other sites? Front Oncol. 2016;6:256.
Zhao J, Chen AX, Gartrell RD, Silverman AM, Aparicio L, Chu T, et al. Immune and genomic correlates of response to anti-PD-1 immunotherapy in glioblastoma. Nat Med. 2019;25:462–9.
Behnan J, Finocchiaro G, Hanna G. The landscape of the mesenchymal signature in brain tumours. Brain. 2019;142(4):847–66. https://doi.org/10.1093/brain/awz044.
Cassetta L, Pollard JW. Targeting macrophages: therapeutic approaches in cancer. Nat Rev Drug Discov. 2018;17(12):887–904. https://doi.org/10.1038/nrd.2018.169.
Poon CC, Sarkar S, Yong VW, Kelly JJP. Glioblastoma-associated microglia and macrophages: targets for therapies to improve prognosis. Brain. 2017;140(6):1548–60. https://doi.org/10.1093/brain/aww355.
Cannarile MA, Weisser M, Jacob W, Jegg AM, Ries CH, Ruttinger D. Colony-stimulating factor 1 receptor (CSF1R) inhibitors in cancer therapy. J Immunother Cancer. 2017;5(1):53. https://doi.org/10.1186/s40425-017-0257-y.
Müller S, Kohanbash G, Liu SJ, Alvarado B, Carrera D, Bhaduri A, et al. Single-cell profiling of human gliomas reveals macrophage ontogeny as a basis for regional differences in macrophage activation in the tumor microenvironment. Genome Biol. 2017;18(1):234. https://doi.org/10.1186/s13059-017-1362-4.
Lee JK, Wang J, Sa JK, Ladewig E, Lee HO, Lee IH, et al. Spatiotemporal genomic architecture informs precision oncology in glioblastoma. Nat Genet. 2017;49(4):594–9. https://doi.org/10.1038/ng.3806.
Darmanis S, Sloan SA, Croote D, Mignardi M, Chernikova S, Samghababi P, et al. Single-cell RNA-seq analysis of infiltrating neoplastic cells at the migrating front of human glioblastoma. Cell Rep. 2017;21(5):1399–410. https://doi.org/10.1016/j.celrep.2017.10.030.
Tirosh I, Venteicher AS, Hebert C, Escalante LE, Patel AP, Yizhak K, et al. Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. Nature. 2016;539(7628):309–13. https://doi.org/10.1038/nature20123.
Venteicher AS, Tirosh I, Hebert C, Yizhak K, Neftel C, Filbin MG, et al. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science. 2017;355:eaai8478.
Yuan J, Levitin HM, Frattini V, Bush EC, Boyett DM, Samanamud J, et al. Single-cell transcriptome analysis of lineage diversity in high-grade glioma. Genome Med. 2018;10:1–15.
Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014;344(6190):1396–401. https://doi.org/10.1126/science.1254257.
Neftel C, Laffy J, Filbin MG, Hara T, Shore ME, Rahme GJ, et al. An integrative model of cellular states, plasticity, and genetics for glioblastoma. Cell. 2019;178(4):835–49 e21. https://doi.org/10.1016/j.cell.2019.06.024.
Diaz A. UCSF Brain Tumor Bioinformatics DAC. Eur Genome-Phenome Arch. https://ega-archive.org/dacs/EGAC00001000515. Accessed 6 May 2019.
Yuan J, Sims P, Levitin HM. Single-cell transcriptome analysis of lineage diversity and microenvironment in high-grade glioma. Gene Expression Omnibus. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103224. Accessed 23 May 2019.
Darmanis S, Quake SR. Single-cell RNA-seq analysis of diffuse neoplastic infiltrating cells at the migrating front of human glioblastoma. Gene Expression Omnibus. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84465. Accessed 14 May 2019.
Tirosh I, Venteicher AS, Suva ML, Regev A. Single cell RNA-seq analysis of IDH-mutant astrocytoma. Gene Expression Omnibus. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE89567. Accessed 29 Mar 2019.
Tirosh I, Venteicher AS, Suva ML, Regev A. Single cell RNA-seq analysis of oligodendroglioma. Gene Expression Omnibus. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70630. Accessed 2 May 2019.
Laffy J, Tirosh I. Single cell RNA-seq analysis of adult and paediatric IDH-wildtype glioblastomas. Gene Expression Omnibus. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131928. Accessed 3 Sept 2019.
Chen A, Zhou W, Yuan J, Sims P, Rabadan R. Single-cell characterization of macrophages in glioblastoma reveals MARCO as a mesenchymal pro-tumor marker. Gene Expression Omnibus. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE141383. Accessed 3 Dec 2019.
Brennan CW, Verhaak RGW, McKenna A, Campos B, Noushmehr H, Salama SR, et al. The somatic genomic landscape of glioblastoma. Cell. 2014;157(3):753. https://doi.org/10.1016/j.cell.2014.04.004.
Wang J, Cazzato E, Ladewig E, Frattini V, Rosenbloom DI, Zairis S, et al. Clonal evolution of glioblastoma under therapy. Nat Genet. 2016;48(7):768–76. https://doi.org/10.1038/ng.3590.
Ceccarelli M, Barthel FP, Malta TM, Sabedot TS, Salama SR, Murray BA, et al. Molecular profiling reveals biologically discrete subsets and pathways of progression in diffuse glioma. Cell. 2016;164(3):550–63. https://doi.org/10.1016/j.cell.2015.12.028.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Ma G, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102:15545–50.
Griffiths JA, Richard AC, Bach K, Lun ATL, Marioni JC. Detection and removal of barcode swapping in single-cell RNA-seq data. Nat Commun. 2018;9(1):2667. https://doi.org/10.1038/s41467-018-05083-x.
Szabo PA, Levitin HM, Miron M, Snyder ME, Senda T, Yuan J, et al. Single-cell transcriptomics of human T cells reveals tissue and activation signatures in health and disease. Nat Commun. 2019;10(1):4706. https://doi.org/10.1038/s41467-019-12464-3.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7. https://doi.org/10.1038/nmeth.3337.
McInnes L, Healy J, Saul N, Großberger L. UMAP: uniform manifold approximation and projection. J Open Source Software. 2018;3(29):861.
Verhaak RGW, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, et al. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell. 2010;17(1):98–110. https://doi.org/10.1016/j.ccr.2009.12.020.
Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003;34(3):267–73. https://doi.org/10.1038/ng1180.
Russo CD, Cappoli N. Glioma associated microglia/macrophages, a potential pharmacological target to promote antitumor inflammatory immune response in the treatment of glioblastoma. Neuroimmunol Neuroinflammation. 2018;5(9):36. https://doi.org/10.20517/2347-8659.2018.42.
Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 2018;19(1):15. https://doi.org/10.1186/s13059-017-1382-0.
Arredouani MS, Palecanda A, Koziel H, Huang YC, Imrich A, Sulahian TH, et al. MARCO is the major binding receptor for unopsonized particles and bacteria on human alveolar macrophages. J Immunol. 2005;175(9):6058–64. https://doi.org/10.4049/jimmunol.175.9.6058.
Jing J, Yang IV, Hui L, Patel JA, Evans CM, Prikeris R, et al. Role of macrophage receptor with collagenous structure in innate immune tolerance. J Immunol. 2013;190(12):6360–7. https://doi.org/10.4049/jimmunol.1202942.
Yang M, Qian X, Wang N, Ding Y, Li H, Zhao Y, et al. Inhibition of MARCO ameliorates silica-induced pulmonary fibrosis by regulating epithelial-mesenchymal transition. Toxicol Lett. 2019;301:64–72. https://doi.org/10.1016/j.toxlet.2018.10.031.
Granucci F, Petralia F, Urbano M, Citterio S, Di Tota F, Santambrogio L, et al. The scavenger receptor MARCO mediates cytoskeleton rearrangements in dendritic cells and microglia. Blood. 2003;102(8):2940–7. https://doi.org/10.1182/blood-2002-12-3651.
La Fleur L, Boura VF, Alexeyenko A, Berglund A, Pontén V, Mattsson JSM, et al. Expression of scavenger receptor MARCO defines a targetable tumor-associated macrophage subset in non-small cell lung cancer. Int J Cancer. 2018;143(7):1741–52. https://doi.org/10.1002/ijc.31545.
Georgoudaki A-M, Prokopec KE, Boura VF, Hellqvist E, Sohn S, Östling J, et al. Reprogramming tumor-associated macrophages by antibody targeting inhibits cancer progression and metastasis. Cell Rep. 2016;15(9):2000–11. https://doi.org/10.1016/j.celrep.2016.04.084.
Lundgren S, Karnevi E, Elebro J, Nodin B, Karlsson MCI, Eberhard J, et al. The clinical importance of tumour-infiltrating macrophages and dendritic cells in periampullary adenocarcinoma differs by morphological subtype. J Transl Med. 2017;15(1):152. https://doi.org/10.1186/s12967-017-1256-y.
Puchalski RB, Shah N, Miller J, Dalley R, Nomura SR, Yoon JG, et al. An anatomic transcriptional atlas of human glioblastoma. Science. 2018;360(6389):660–3. https://doi.org/10.1126/science.aaf2666.
Leblond MM, Gerault AN, Corroyer-Dulmont A, MacKenzie ET, Petit E, Bernaudin M, et al. Hypoxia induces macrophage polarization and re-education toward an M2 phenotype in U87 and U251 glioblastoma models. Oncoimmunology. 2016;5(1):e1056442. https://doi.org/10.1080/2162402X.2015.1056442.
Gao D, Vahdat LT, Wong S, Chang JC, Mittal V. Microenvironmental regulation of epithelial-mesenchymal transitions in cancer. Cancer Res. 2012;72(19):4883–9. https://doi.org/10.1158/0008-5472.CAN-12-1223.
Chavez-Galan L, Olleros ML, Vesin D, Garcia I. Much more than M1 and M2 macrophages, there are also CD169+ and TCR+ macrophages. Front Immunol. 2015;6:263.
Gartrell RD, Marks DK, Hart TD, Li G, Davari DR, Wu A, et al. Quantitative analysis of immune infiltrates in primary melanoma. Cancer Immunol Res. 2018;6(4):481–93. https://doi.org/10.1158/2326-6066.CIR-17-0360.
Cloughesy TF, Mochizuki AY, Orpilla JR, Hugo W, Lee AH, Davidson TB, et al. Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nat Med. 2019;25(3):477–86. https://doi.org/10.1038/s41591-018-0337-7.
Rivera LB, Bergers G. Location, location, location: macrophage positioning within tumors determines pro- or antitumor activity. Cancer Cell. 2013;24(6):687–9. https://doi.org/10.1016/j.ccr.2013.11.014.
Milne SA, McGregor AL, McCulloch J, Sharkey J. Increased expression of macrophage receptor with collagenous structure (MARCO) in mouse cortex following middle cerebral artery occlusion. Neurosci Lett. 2005;383(1-2):58–62. https://doi.org/10.1016/j.neulet.2005.03.065.
Lavin Y, Kobayashi S, Leader A, Amir E-aD, Elefant N, Bigenwald C, et al. Innate immune landscape in early lung adenocarcinoma by paired single-cell analyses. Cell. 2017;169:750-65.e17.
Han J, Alvarez-Breckenridge CA, Wang QE, Yu J. TGF-β signaling and its targeting for glioma treatment. Am J Cancer Res. 2015;5(3):945–55.
Matsushita N, Komine H, Grolleau-Julius A, Pilon-Thomas S, Mulé JJ. Targeting MARCO can lead to enhanced dendritic cell motility and anti-melanoma activity. Cancer Immunol Immunother. 2010;59(6):875–84. https://doi.org/10.1007/s00262-009-0813-5.
Filley AC, Henriquez M, Dey M. Recurrent glioma clinical trial, CheckMate-143: the game is not over yet. Oncotarget. 2017;8(53):91779–94. https://doi.org/10.18632/oncotarget.21586.
Saha D, Martuza RL, Rabkin SD. Macrophage polarization contributes to glioblastoma eradication by combination immunovirotherapy and immune checkpoint blockade. Cancer Cell. 2017;32(2):253–67 e5. https://doi.org/10.1016/j.ccell.2017.07.006.
Arredouani MS. Is the scavenger receptor MARCO a new immune checkpoint? OncoImmunology. 2014;3:e955709-1-e-3.
Chen A. Single-cell function library. Github. https://github.com/RabadanLab/MARCOsinglecell. Accessed 16 Mar 2021.
We would like to thank Emanuelle M. Rizk and Tejaswi Sudhakar for the assistance with sample collection, Chia-Ing Jan for the support with microscopy, and the Molecular Pathology Core at CUIMC for the help with IF staining.
This work has been funded by NIH grants R01 CA185486 (R.R.), R01 CA179044 (R.R.), R01 NS103473 (P.A.S., J.N.B., P.C.), U54 CA193313, U54 CA209997 (R.R., P.A.S.), KL2TR001874 (R.D.G), NSF/SU2C/V-Foundation Ideas Lab Multidisciplinary Team (PHY-1545805) (R.R.), and 2018 Stand Up To Cancer Phillip A. Sharp Innovation in Collaboration Awards (R.R.). A.X.C. is funded by the Medical Scientist Training Program (T32GM007367). R.D.G. is supported by Swim Across America and Hyundai Hope on Wheels Hope Scholar Award.
Ethics approval and consent to participate
The GBM tissue specimens were collected from surgical resections of de-identified patients at Columbia University Irving Medical Center, who provided informed consent to participate in these studies through a protocol (IRB-AAAJ6163) approved by the Columbia University Institutional Review Board. The research was conducted in accordance with the principles of the Declaration of Helsinki.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Chen, A.X., Gartrell, R.D., Zhao, J. et al. Single-cell characterization of macrophages in glioblastoma reveals MARCO as a mesenchymal pro-tumor marker. Genome Med 13, 88 (2021). https://doi.org/10.1186/s13073-021-00906-x
- Single-cell RNA-seq
- Cancer immunotherapy