Single-Cell Features of Glioblastoma Macrophages Reveal MARCO as a Stroma-Related Tumor Marker

/Abstract/

Macrophages are the most common infiltrating immune cells in gliomas, exhibiting various pro-tumor and anti-tumor effects. However, the different subpopulations of macrophages and their impact on the tumor microenvironment remain unclear.

We combined new and previously published single-cell RNA sequencing data from 66 gliomas, analyzing 19,331 macrophages from a total of 98,015 single cells.

Unsupervised clustering revealed a pro-tumor subpopulation of bone marrow-derived macrophages characterized by the scavenger receptor MARCO, which was found almost exclusively in IDH1 wild-type glioblastomas. Previous studies indicated that MARCO is a marker of tumor malignancy in melanoma and non-small cell lung cancer; in this study, we found that high MARCO expression correlates with poor prognosis and mesenchymal subtypes. Moreover, MARCO expression levels significantly changed in a treatment-dependent manner during anti-PD1 checkpoint inhibitor therapy, which we confirmed through immunofluorescence imaging.

These findings indicate the presence of a novel macrophage subpopulation in brain gliomas that drives the tumor progression of glioblastoma and may serve as a potential therapeutic target.

/Background/

Glioblastoma (GBM) is a refractory primary brain tumor. Despite standard treatments involving surgery, chemotherapy, and radiotherapy, GBM ultimately recurs, with a median survival of only about 16 months. One reason for the poor treatment outcomes is the complexity of the role of macrophages in the tumor microenvironment. While immunotherapy has been successful in various other cancers, the immunosuppressive microenvironment in GBM, including tumor-associated macrophages (TAMs), hinders the efficacy of immunotherapy. A recent study on checkpoint inhibitor therapy for GBM found an association between the infiltration of HLA class II-deficient macrophages and poor treatment response to immunotherapy.

Targeting these TAMs (e.g., using CSF1R inhibitors) is a potentially effective treatment strategy, but it requires a better understanding of the functional specificity and necessary markers of TAMs. Our traditional understanding of macrophage polarization (M1 vs. M2) is a generalized statement about macrophages in different polarized states; the actual situation is much more complex. Furthermore, in GBM, there are various monocyte lineage cells (commonly referred to as “macrophages”), including bone marrow-derived macrophages (BMDMs) recruited from the blood and tissue-resident microglia. Therefore, the specific markers and pathways involved in pro-tumor macrophages in GBM remain unclear.

Due to the diversity of cell types in the tumor microenvironment, tumor bulk expression profiles are not ideal for studying cell subpopulations. In contrast, single-cell RNA sequencing (scRNA-seq) has proven helpful in understanding the heterogeneity within GBM. In recent years, several scRNA-seq datasets have been published in gliomas, aiding in elucidating the general differences between BMDMs and microglial populations in the tumor microenvironment. However, each study had a limited number of patients. In this paper, we combined nine new scRNA-seq samples from glioblastomas with 57 previously published cases, exploring the distribution of macrophages in gliomas on an unprecedented scale. In the process, we identified a novel pro-tumor macrophage marker (MARCO) in GBM through immunofluorescence imaging. In addition to studying the clinical impact of MARCO expression levels on survival and immunotherapy, we also elucidated the relationship between MARCO expression levels and mesenchymal, hypoxic, and anti-inflammatory features.

/Methods/

Published single-cell RNA sequence count matrices were obtained from multiple databases and combined with nine previously unpublished cases, forming a cohort of 50 GBM and 16 LGG samples (Supplementary File 1: Table S1). Tumor bulk expression profile data and survival data were obtained from 528 GBM patients in TCGA-GBM and 75 GBM cases from Wang et al. Expression subtypes and IDH1 mutation status of TCGA were obtained from Ceccarelli et al. Other gene sets were obtained from MSigDB v6.2, in addition to the BMDM microglial cell gene set from Yuan et al.

First, all datasets were filtered to remove mitochondrial and ribosomal proteins. The datasets were then merged (for GBM and LGG separately) while retaining the intersection of genes and discarding genes with zero total counts. The raw counts were normalized to log2(1+TPK), as described by Yuan et al. The expressed genes were used as a filtering list to reduce batch effects, intersecting with LM22 from CIBERSORT. For visualization, principal component analysis (PCA) was first applied to reduce the total dimensions to 5% of the number of genes, and then UMAP clustering was performed with default parameters to non-linearly reduce to two-dimensional data. Three cell types were identified using standard markers for macrophages, T lymphocytes, and tumor cells: CD14, CD3, and SOX2. Macrophages were isolated, and the same dimensionality reduction procedure was applied to generate macrophages. K-means clustering analysis was used to separate macrophage populations, with scores from k=2 to k=6, where k=2 provided the highest score. This filtering list was used to evaluate the enrichment gene sets characterizing each of the two macrophage clusters.

To validate the presence of MARCO in GBM patient tissues, we used double-staining immunofluorescence (IF) to co-stain MARCO with CD163. Tissue specimens from one patient of each tumor subtype (including IDH1 wild-type GBM, IDH1 mutant GBM, and grade III anaplastic astrocytoma (LGG)) were evaluated. Additionally, two patients receiving anti-PD1 inhibitor treatment were compared before and after treatment, including one responder and one non-responder. Anti-MARCO (1:100, PA5-64134, Invitrogen (Carlsbad, CA)) was used, followed by biotinylated goat anti-rabbit IgG (1:200, BA1000, Vector Laboratories) and streptavidin conjugated 594 (1:1000, S11127, Invitrogen) along with anti-CD163 (1:100, 10D6, Biocare Medical (Pacheco, CA)), biotinylated horse anti-mouse IgG (1:200, BA2000, Vector Laboratories), and streptavidin conjugated 488 (1:300, Invitrogen) on 5μm sections for IF staining. The stained sections were examined and photographed using a Nikon Eclipse TE2000-E. Synthetic images were captured from each area of interest at emission wavelengths of 395, 488, and 568. All images in Figure 4c were thresholded for the MARCO channel using a fixed value. For single-staining controls, all channels were thresholded using a fixed value.

TCGA-GBM U133 microarray data and preprocessed expression data from Wang et al. were used to determine MARCO expression. Each cohort was first independently log and Z-score normalized before merging. These normalized values were combined with survival data, including overall survival and disease-free survival. Survival differences were assessed in two ways: (1) dividing MARCO expression at the median into high and low expression groups and comparing survival curves; (2) directly generating a univariate Cox model based on MARCO expression levels and using the Wald p-value of MARCO as a covariate. In Figures 2a and b, we only used cases known to be IDH1 wild-type.

Patel et al.’s single-cell subtype scoring method was used to calculate characteristic scores for each Verhaak expression type. The same method was used to obtain BMDM and microglial cell characteristics from Yuan et al., as well as to assess the comprehensive score of tumor-macrophage interaction genes.

GSEA was performed using the official Linux client v4.0.1. The normalized gene expression of all macrophages was used as input data, with MARCO expression levels serving as phenotype markers, and Pearson correlation used as the metric for gene ranking. Default weights (p=1) were applied. P-values were calculated based on 1000 permutations of the phenotype, and the final gene sets were ranked according to enrichment scores.

The average normalized expression level of MARCO in the macrophage population of each GBM sample was calculated and compared to a selected list of recruiting factors in the tumor population of the same samples: macrophage colony-stimulating factor (CSF1), granulocyte/macrophage colony-stimulating factor (CSF2), hepatocyte growth factor (HGF), monocyte chemotactic protein 1 (MCP-1), macrophage migration inhibitory factor (MIF), stromal-derived factor 1 (SDF-1), transforming growth factor beta (TGF-β, including TGFB1, 2, and 3), interleukin 10 (IL-10), osteopontin (SPP1), and lactoferrin (MFGE8). Gene set expression subtype scores were calculated using this gene set, and the tumor cells of each case were averaged. The score was then compared with the corresponding MARCO expression level in the macrophage population using Spearman correlation. The same procedure was used to assess the correlation between macrophage MARCO expression levels and paired tumor cell expression subtype scores.

To quantify the utility of MARCO as a marker compared to traditional myeloid markers, BMDM markers, and microglial markers, we assessed the standardized dispersion coefficient of each gene. We used the Highly_Variable_Genes function in SCANPY with default parameters and batch correction. We determined p-values using the exact permutation test.

All statistical analyses were conducted in Python 3.6. In all box plots, the center line represents the median, with lower and upper limits representing the first and third quartiles, respectively. The shaded area in survival curves represents the confidence interval. Violin plots used a Gaussian kernel to estimate density, with the center line indicating the median. Non-parametric comparisons were made between two groups. Statistical significance was evaluated at an adjusted p-value threshold of 0.05.

/Results/

To understand the heterogeneity of cell populations in GBM, we collected single-cell RNA-seq data from nine GBM cases and combined it with previously published studies from 41 GBM patients, resulting in a total of 79,968 single-cell transcriptomes (Supplementary File 1: Table S1). To reduce batch effects, we filtered the gene list to 499 genes overlapping with LM22, which is the reference matrix used by CIBERSORT to differentiate immune cells (Supplementary File 2: Table S2). After filtering, log normalization, batch effect reduction, and dimensionality reduction, these cells visualized to produce a macrophage population characterized by CD14 (17,132 cells; Figure 1a, Supplementary File 3: Figure S1). Unsupervised k-means clustering on this macrophage population revealed two subpopulations (Supplementary File 3: Figure S2). One subpopulation was enriched with inflammatory-related genes, with CCL4 and IL1B being the top two genes on the filtered gene list (p=0.004 and p=0.008, exact permutation test, n=499 genes; Figure 1b). The other subpopulation, in contrast to the inflammatory side, had MARCO (a macrophage receptor with collagen structure) as a highly expressed gene (p=0.004, exact permutation test, n=499 genes; Figure 1b). Even ignored during dimensionality reduction, MARCO remained the most abundant gene on this side, proving its representativeness (p=0.004, exact permutation test, n=499 genes). Observing the expression levels of MARCO across all cell populations, we found it specifically expressed in the macrophages of this subpopulation (Supplementary File 3: Figure S3A). MARCO is a scavenger receptor typically present on alveolar macrophages, with various immunomodulatory functions. Given the negative role of MARCO in other cancers, we chose to focus on these MARCO+ macrophages and their impact on the tumor microenvironment.

Next, we investigated whether this macrophage subpopulation could also be observed in low-grade gliomas (LGG), including grade II and III astrocytomas and oligodendrogliomas. We collected single-cell expression profiles from 18,047 cells from 16 LGG cases and processed them in the same manner as the GBM above (Figure 1c). Although 2,199 macrophages were identified, only 1% (23 cells) showed non-zero MARCO expression (Figure 1d), whereas in GBM, 12% (2,092 out of 17,132 cells, p<0.001; chi-squared test). Although GBM samples had lower library complexity (p<0.001, n=49 patients, Mann-Whitney U test; Supplementary File 3: Figure S4), this situation still existed. Since most LGGs are IDH1 mutant, we assessed the MARCO expression levels in four IDH1 mutant GBM cases. Among the 281 macrophages from IDH1 mutant GBM, only one cell exhibited non-zero MARCO expression. Comparing standardized expression levels, the average MARCO expression level in IDH1 wild-type GBM macrophages was 160 times that of IDH1 mutant GBM (p=0.020, n=41 patients, Mann-Whitney U test), and 64 times that of LGG (p=0.0007, top of Figure 1e, n=49 patients). Meanwhile, the average expression levels of other macrophage markers such as CD14 showed no significant differences (p=0.34 and p=0.30, bottom of Figure 1e, n=41 and n=49 patients). Immunofluorescence imaging in IDH1 wild-type cases confirmed the co-localization of MARCO with CD163 macrophages in these patient tissues, while IDH1 mutant and LGG tissues showed rare co-localization (if any) (Supplementary File 3: Figures S5 and S6). Therefore, MARCO expression was found almost exclusively in the macrophages of IDH1 wild-type GBM, rather than in the less lethal IDH1 mutant or low-grade gliomas.

To understand the clinical outcomes of MARCO expression levels, we assessed its tumor bulk expression in two published datasets: TCGA-GBM and the GBM cohort from Wang et al. (a total of 603 patients). Dividing patients into high and low MARCO expression groups, we found that MARCO expression levels negatively correlated with overall survival (OS: p=0.0046, n=592 patients, log-rank test, Supplementary File 3: Figure S7A) and disease-free survival (DFS: p=0.018, n=387 patients). However, given the association of MARCO with IDH1 and the significant impact of IDH1 mutation on survival, we repeated the analysis solely on 437 cases with known IDH1 wild-type status, yielding similar results (OS: p=0.0084, n=437 patients; DFS: p=0.035, n=267 patients, log-rank test, Figures 2a, b). Continuing the analysis on the IDH1 wild-type cohort, we treated MARCO expression levels as a continuous variable and found that the impact of MARCO on survival was also significantly adverse in the univariate Cox model (OS patients p=0.00057, n=437, hazard ratio 1.19; DFS patients p=0.0047, n=267, hazard ratio 11.9; Wald test). These effects were most pronounced in long-term survival rates, with two-year (18.0% vs. 31.8%) and five-year (3.7% vs. 7.8%) overall survival rates for MARCO high patients being only half that of MARCO low patients. Similarly, the two-year DFS in the MARCO high population decreased more than threefold (5.8% vs. 18.1%; no five-year DFS data). The impact of MARCO expression levels on survival remained significant after controlling for common factors including age and MGMT methylation status (Supplementary File 3: Figure S7C).

Consistent with the single-cell sequencing results, we found in tumor bulk sequencing data that IDH1 mutations were associated with decreased MARCO expression levels (p=0.0039, n=482 patients, Mann-Whitney U test, Figure 2c). This was also important for dividing the population into high and low MARCO groups (p=0.0016, n=482, Fisher’s exact test, Figure 2d). We then compared MARCO expression with transcriptomic subtypes. Although MARCO is not listed among the mesenchymal gene set, its expression was highly enriched in mesenchymal samples (all paired comparisons between mesenchymal and other subtypes, p<0.00001, n=565 patients; Figure 2e). These tumor bulk transcriptomic results support our single-cell findings that MARCO expression is abundant in IDH1 wild-type tumors and also indicate that MARCO expression levels are associated with poor prognosis and unfavorable mesenchymal subtypes.

To reveal the source of this mesenchymal feature, we returned to the single-cell expression data, where mesenchymal expression features were initially found in the macrophage population rather than in tumor cells (p<0.00001, n=50 patients, Mann-Whitney U test, Figure 3a, Supplementary File 3: Figure S8, Supplementary File 2: Table S3), particularly in macrophages expressing MARCO (p=0.005). The average expression level of MARCO in macrophages was also significantly correlated with the mesenchymal feature scores in tumor cells from the same samples (p=0.0084, n=41 patients, Spearman correlation; see “Methods” section; Figure 3b). We then performed gene set enrichment analysis (GSEA) based on the Pearson correlation of MARCO expression levels with all other genes to understand the underlying processes in this subpopulation. Among the 50 hallmark gene sets in MSigDB, the gene set with the highest enrichment score in MARCO+ macrophages was epithelial-mesenchymal transition, angiogenesis, glycolysis, and hypoxia (all four gene sets FDR Q < 0.001, default weighted phenotype permutation test, n=17,132 cells; Figure 3c, d). Single-cell data from Darmanis et al. supported the association with hypoxia, where samples were taken from the core and periphery of the same four tumors. Macrophages from the tumor core had higher MARCO expression levels than those from the tumor periphery (p<0.001, n=1,846 cells; Mann-Whitney U test; Figure 3e). Furthermore, we examined the expression levels of MARCO in the Ivy GAP database, which includes laser microdissection samples from different anatomical structures. We found that MARCO expression levels varied significantly across various structures (p=0.014, n=270 samples from 37 patients; Kruskal-Wallis test, Supplementary File 3: Figure S9), with the highest expression level in the perivascular area of the tumor. These findings support that MARCO+ macrophages reside in the hypoxic tumor core with GSEA characteristics.