Transcriptional landscape of human microglia implicates age, sex, and APOE‐related immunometabolic pathway perturbations

Abstract Microglia have fundamental roles in health and disease; however, effects of age, sex, and genetic factors on human microglia have not been fully explored. We applied bulk and single‐cell approaches to comprehensively characterize human microglia transcriptomes and their associations with age, sex, and APOE. We identified a novel microglial signature, characterized its expression in bulk tissue and single‐cell microglia transcriptomes. We discovered microglial co‐expression network modules associated with age, sex, and APOE‐ε4 that are enriched for lipid and carbohydrate metabolism genes. Integrated analyses of modules with single‐cell transcriptomes revealed significant overlap between age‐associated module genes and both pro‐inflammatory and disease‐associated microglial clusters. These modules and clusters harbor known neurodegenerative disease genes including APOE, PLCG2, and BIN1. Meta‐analyses with published bulk and single‐cell microglial datasets further supported our findings. Thus, these data represent a well‐characterized human microglial transcriptome resource and highlight age, sex, and APOE‐related microglial immunometabolism perturbations with potential relevance in neurodegeneration.


| INTRODUC TI ON
Microglia are the resident macrophages of the central nervous system (CNS), responsible for clearance of cellular debris and pathological protein aggregates. In the healthy brain, they exist in a homeostatic state and can be induced to a reactive state in response to changes in the CNS microenvironment, such as inflammation and neuronal damage (Masuda et al., 2020). They are fundamental to maintaining brain homeostasis during development, aging, and disease; therefore, microglial dysfunction could ultimately lead to neurodegeneration (Li & Barres, 2018). Microglia are integral to the pathophysiology of neurodegenerative diseases, including Alzheimer's disease (AD) and multiple sclerosis, with chronic inflammation implicated as a contributing factor (Hickman et al., 2018;Keren-Shaul et al., 2017;Krasemann et al., 2017).
Fresh human brain tissue studies are imperative to the characterization of the microglial transcriptome in health and disease; however, accessibility is limited. Although single nuclei studies using frozen tissue provide an easier alternative, recent studies have demonstrated limitations in detecting substantial populations of less abundant cell types (Del-Aguila et al., 2019;Mathys et al., 2019). Additionally, it was recently reported that many microglial activation genes are expressed in the cytosol and therefore are likely to be missed by single nuclei RNA sequencing (snRNAseq; Thrupp et al., 2020). Recent single-cell studies aiming to characterize microglial gene expression using fresh tissue have highlighted the heterogeneity in microglial phenotypes Olah et al., 2020;Sankowski et al., 2019). This has revealed that phenotypic changes are not binary but rather a spectrum of states in which microglia can simultaneously co-exist during transition from homeostatic to more reactive states. Additionally, these different subsets could have specialized functions in brain homeostasis and dysfunction. Thus, it is increasingly important to characterize these heterogeneous subpopulations to understand their roles in health and disease. This could also help facilitate the design of novel therapeutic approaches to target-specific subpopulations of cells and modulate their activity (Li & Barres, 2018).
Obtaining fresh human tissue from neurosurgeries allows us to study the mechanisms of microglial function in living cells.
Unfortunately, this tissue is usually excised from surgical procedures for tumor resection or relieving temporal lobe epilepsy, rendering it difficult to distinguish between normal and disease-affected tissue. Darmanis et al. (2017) investigated the effect of GBM tumors on CNS cell types and surrounding tissue, revealing that myeloid cells are greatly affected by the tumor microenvironment. They found that peri-tumor myeloid populations were primarily pro-inflammatory microglia compared to macrophages within the tumor core (Darmanis et al., 2017). In temporal lobe epilepsy, two distinct microglial phenotypes have been identified with microglia present in sclerotic areas with few neurons expressing markers of activation, including anti-inflammatory cytokine IL10 (Kinoshita & Koyama, 2021;Morin-Brureau et al., 2018). The other phenotype occurs transiently following a seizure, with secretion of interleukins CXCL8 and IL1B mediated by the NLRP3 inflammasome (Morin-Brureau et al., 2018).
Microglial expression has also been shown to be affected by aging (Galatro et al., 2017;Olah et al., 2018); however, few studies have investigated the effects of sex and genetic factors on human microglia. Sex differences in microglia have been previously reported in mice, with females being predisposed to harboring more activated microglia than males (Frigerio et al., 2019;Nelson et al., 2017;Stephen et al., 2019). APOE, a lipoprotein of which the ε4 allele (APOE-ε4) is a major risk factor for AD and also implicated in other neurodegenerative diseases (Yamazaki et al., 2019]), is upregulated in disease-associated microglia (DAM) in mice and humans, but downregulated in astrocyte and oligodendrocyte subpopulations (Grubman et al., 2019;Hammond et al., 2019;Keren-Shaul et al., 2017;Mathys et al., 2019). In microglia and neurons, APOE interacts with LDL receptors to facilitate endocytosis of cholesterol and phospholipids and modulate lipid homeostasis in the brain (Gamache et al., 2020). Such studies provide growing support for cell typespecific functions of APOE; however, its effects on microglia remain to be fully elucidated. Thereby, identifying age, sex, and APOEassociated pathways in microglia will provide greater insight into the functions of specific microglial subsets in relation to these risk factors. Interindividual variability and diversity in functional states makes targeting specific microglial subsets in disease challenging for modulating these cells (Li & Barres, 2018). Identifying the mechanisms regulating microglial homeostasis and activation can allow us to manipulate these cells for therapeutic purposes.
In this study, we leveraged both bulk and single-cell approaches to provide a comprehensive characterization of the adult human microglial transcriptome. We obtained fresh intraoperative neurosur-| 3 of 21 PATEL ET AL.

| RE SULTS
To uncover microglial transcriptional profiles and their associations with age, sex, and APOE, we performed microglial cell typespecific and single-cell RNA sequencing (scRNAseq) studies in fresh human brain tissue. We obtained neurosurgical tissue unaffected by the primary disease process from 19 human donors ( Figure S1). Microglia were isolated by CD11b + microbead selection followed by FACS sorting of cells expressing the CD11b + / CD45 intermediate microglial signature to acquire a more purified population. These samples underwent bulk microglia RNAseq, with subsets of these and additional samples also undergoing 10x scRNAseq (n = 5) and bulk tissue RNAseq (n = 9; Figure 1a; Table 1). Validation of sorted microglia using qPCR showed the expected CD11b + /CD45 intermediate /P2RY12 + microglial signature (Li & Barres, 2018) with no expression of other cell type markers, indicating that we isolated a highly enriched microglial population ( Figure S2).

F I G U R E 1
Characterization of our core human microglial signature. (a) Schematic illustrating our experimental approach for isolating microglial populations from fresh brain tissue and data analyses. [Created with BioRender.com] (b) MSigDB GO terms enriched in our microglial signature genes and top 25 genes for each. (c) Venn diagram showing number of overlapping genes between our microglial signature and those previously reported from Galatro et al. (2017), Gosselin et al. (2017) and Olah et al. (2018) To define a core human microglial signature, we calculated log 2 fold change and q-values of differential expression for each gene between bulk microglia RNAseq data in our study and bulk brain RNAseq data from 7 AMP-AD datasets provided by Mayo Clinic (Allen et al., 2016]), Mount Sinai Brain Bank ,  (Table S7). GO enrichment using MSigDB showed that this signature was enriched for genes involved in immune-related and inflammatory response pathways as would be expected, and leukocyte-mediated immunity ( Figure 1b).
To determine the ability of bulk brain tissue data to capture microglial genes, we assessed the expression levels of our microglial signature genes in each of the 7 AMP-AD bulk brain RNAseq datasets. Of the 1971 microglial signature genes in our study, 37-47% were captured in these bulk brain datasets ( Figure S3A-B). Our microglial signature genes comprised 3.6-4.5% of the expressed bulk brain transcriptome, consistent with prior estimations (Mathys et al., 2019;Wang et al., 2020). We next compared bulk microglia RNAseq transcript levels to that obtained from bulk tissue RNAseq of neurosurgical fresh brain tissue samples. Bulk fresh brain tissue does not capture all microglial marker genes, as demonstrated by the low correlation between bulk tissue and bulk microglia data (R = 0.46, Table S6). This reiterates the need for complementary single-cell type data to deconvolute cell typespecific expression. We provide the list of microglial signature genes that are also expressed at high levels in bulk brain tissue data (Table   S5), which can serve as a validated resource for microglial signature gene markers in bulk RNAseq datasets.
To determine how the microglial signature in this study compared where 49.7% of their genes were also present in our signature, and 22% of ours in their signature. Gosselin et al. (2017) samples were also obtained from neurosurgical tissue resections like our cohort and are unlike Galatro et al. (2017) and Olah et al. (2018) samples that were harvested during autopsy. Although there appears to be a common set of microglial genes consistent across signatures, each also harbors many unique genes, which could be due to study or individual specific differences.
2.2 | Transcriptional profiling of microglia discovers co-expression networks and implicates lipid and carbohydrate metabolism pathways associated with age, sex, and APOE We generated gene co-expression networks using WGCNA (Langfelder & Horvath, 2008) to reduce number of tests and increase power to detect genetic associations with age, sex, and APOE.
We identified 7 modules with significant associations (Figure 2; Figure S5;   Rangaraju et al., 2018b), CTSL, (Cermak et al., 2016) and UNC5B (Ahn et al., 2020;Xu et al., 2016). Due to the nature of the tissue, we compared expression levels of the key module genes within tumor and epilepsy samples ( Figure S6). We observed significant differences due to diagnosis in NPC2, PLAU, APOC1, and IKBKE; therefore, it is important to note that some of their associations may be confounded due to the disease state.

| Single-cell transcriptome reveals specific subtypes of microglia
To uncover distinct microglial subtypes, a subset of sorted microglial samples from neurosurgical brain tissue underwent single-cell expression profiling. We obtained 26,558 cells from 5 unique individuals, including one individual who underwent epilepsy surgery and had samples from two brain regions (Table 1). Analysis of the scRNAseq data from these samples revealed 13 distinct cell clusters, which were annotated using established neuronal and glial marker genes from the literature ( Table S11). For these samples, the most marked difference was observed for macrophages (cluster 9) and homeostatic microglia (cluster 2), which had greater contributions from the mesiotemporal and anterior temporal regions, respectively. This could be due to the proximity of the mesiotemporal sample to the disease-affected region.
We characterized the microglial clusters by their expression of es-  Figure S8; Table S12). Most of these markers are distinct to a single cluster, although some clusters appeared to have similarities in their marker expressions. To define the proximity of their transcriptional profiles, we performed hierarchical clustering of the microglial clusters ( Figure 3d). We determined that the homeostatic microglia cluster 2 was transcriptionally closest to clusters 7 and 11, which may represent subtypes of homeostatic microglia. Cluster 11 is enriched for markers of cell proliferation (STMN1, H2AFZ, PCNA, and MKI67), some of which were also observed by Olah et al. (2020), suggesting that these could be proliferating microglia. Clusters 1 and 6 both expressed inflammatory chemokines CCL2 and CCL4 and antiinflammatory molecule EGR2; however, cluster 6 was more closely related to DAM, whereas cluster 1 (named as Leukocyte-recruiting cluster) represented a more pro-inflammatory signature with greater expression of IL6 and TNFα. Cluster 6 (named as interferon-response cluster) also highly expressed interferon-related marker IFITM3 and ISG15, also observed in a cluster by Olah et al. (2020), which they defined as an interferon response-enriched subset. The upregulation of chemokines and interleukins in these clusters suggest that they could be involved in recruitment of other immune cells. These findings highlight different transcriptional profiles for these two inflammatory clusters that may represent distinct activated microglia subtypes. Cluster 3 highly expressed heat shock protein HSPA1A, an immediate early gene (Schmunk et al., 2020) reportedly involved in antigen processing (Aung et al., 2012), response to stress and injury and exhibiting decreased gene expression in multiple sclerosis patients (Gandhi et al., 2010;Satoh et al., 2005). Several were upregulated in this cluster, suggesting that this may represent cells that underwent dissociation-induced stress .
Several of the clusters did not express well known existing cell type markers. Clusters 5/8 and FOS + 0/4 were transcriptionally closest to one another (Figure 3d). Cluster 5 has distinct expression of immunoreactive marker CD163, which was not observed in other subsets except macrophages. Several HLA genes are also highly expressed in this cluster, suggesting that these may be immunoreactive microglia (Hendrickx et al., 2017). Cluster 8 marker FTL has recently been used to characterize iron-accumulating microglia (Kenkhuis et al., 2021). Our findings highlight transcriptional profiles for known microglial clusters, describe the transcriptional proximity of these clusters and suggest that less well-defined clusters could potentially represent novel or intermediate transcriptional states of microglia.
Our microglial signature was significantly enriched in more clusters expressing more activated markers (Table S13) full enrichment statistics provided in Table S15). Age-associated co-expression network ME14, implicated in lipid metabolism, was significantly enriched in interferon-response (cluster 6) and DAM (cluster 10) clusters. Genes within module 28, which was significantly upregulated with APOE-ε4, had statistically significant enrichment in all clusters except cluster 7. There was no statistically significant enrichment for any of the other microglial modules that had significant age, sex, or APOE associations, suggesting that these factors may have ubiquitous effects on most microglial subtypes. Some of the remaining microglial co-expression networks had distinct patterns of cluster enrichment ( Figure S9), suggesting that some but not all networks could be representative of distinct microglial subtypes.

| Meta-analyses with published datasets support age, sex, and APOE associations
Each dataset was individually processed through the same MAPR-seq pipeline for quality control to minimize variability due to data processing. Meta-analysis of WGCNA results was performed by coercing the external dataset co-expression networks onto our existing co-expression networks. The forest plot in Figure 4a highlights the individual and combined associations for networks of interest across the datasets where age, sex, or APOE information was available (Table   S17). For module ME14 associations with age, Galatro et al. (2017) also exhibited a similar direction of effect (R = 0.30, p = 0.06), whereas Srinivasan et al. (2020) did not show any association (R = −0.003, p = 0.99), likely due to a higher median age of individuals in the study.
When meta-analyzed, sex-associated ME26 genes were significantly downregulated in males in all datasets, supporting our finding. This module was also significantly associated with APOE e4 carriers in the Srinivasan dataset (R = 0.42, p = 0.04). However, Srinivasan samples were inversely correlated with APOE for ME23 in comparison to our data (R = 0.61, p = 0.07). Srinivasan et al. (2020) was inversely correlated with age and APOE for ME23, but not with sex or APOE for ME26. Overall, it appears that the direction of effects was similar for correlation with these traits in most datasets ( Figure S10).  (Table S18). This is similar to what we observed in our single-cell data (Table S13), supporting the notion that our microglial signature is enriched for activated cells.
In this study, we sought to overcome these knowledge gaps by characterizing the transcriptome of sorted bulk and single-cell microglial populations isolated from fresh human brain tissue. We iden- Our findings in sorted cells from fresh human brain tissue provide transcriptional evidence for immunometabolism changes and proinflammatory phenotype with microglial aging, thereby contributing essential complementary data from humans for this cell type.
Besides module ME14, we determined that ME26 is also enriched for lipid metabolism genes. ME26 module expression is higher in both APOE-ε4 and female sex; however, we note that in our sorted bulk microglia RNAseq samples, there were no male APOE-ε4 carriers.
Therefore, the distinct influence of sex and APOE on the expression of this module remains to be established. APOE-ε4, a major risk factor for Alzheimer's disease, has the lowest lipid binding efficiency compared with other APOE isoforms (Chausse et al., 2020). Increased cholesterol accumulation has been reported in both iPSC-driven astrocytes from APOE-ε4 carriers (Lin et al., 2018) and in Apoe-deficient microglia (Nugent et al., 2020). These findings collectively support a role for APOE-ε4 associated microglial transcriptional changes and disrupted cholesterol metabolism. Using our sorted microglia RNAseq data, we identified five additional modules that associate with APOE-ε4, one in a positive direction (ME28) and four negatively (ME4, ME23, ME34, and ME36). Of these, module ME23 had the second most significant APOE-ε4 association after ME26. Interestingly, ME23 was enriched for carbohydrate metabolism biological processes, which are also tightly regulated in microglia (Bernier et al., 2020). Module ME23 harbors known AD risk genes BIN1 and PLCG2, where the latter is a microglial gene that modulates signaling through TREM2  and also a hub gene in this module. ME23 genes BIN1, JUN, and TGFBR2 were found to be reduced in a mouse microglial neurodegenerative phenotype gene signature (Krasemann et al., 2017). These findings further demonstrate the consistency of our human microglial data with that from mouse models and supports perturbed microglial immunometabolism as a potential pathogenic mechanism in neurode-  (Table   S11). As there were only 5 cells in this cluster, this is likely to denote contamination from the sorting procedures. We identified microglial clusters that share characteristics of those previously reported in mice (Keren-Shaul et al., 2017) and humans Olah et al., 2020), such as DAM. We also uncovered clusters that had not been previously characterized in the literature, including cluster 7, exhibiting high expression levels of astrocytic marker gene SLC1A3. Microglial expression of SLC1A3 was previously shown to occur in mice and humans especially in disease states (Chrétien et al., 2004;Grassivaro et al., 2020;Wilhelmsson et al., 2017). We also leveraged these scRNAseq data to further characterize the sorted bulk microglial expression modules. Integration of our samples with Olah et al. allowed us to explore microglial subtypes at single-cell resolution with increased power. We found that even in this larger dataset, our microglial signature was similarly enriched in the more activated subtype clusters, demonstrating the robustness of this signature across datasets. The homogeneity in clustering highlighted shared expression patterns between cells from both datasets, providing further support for the putative microglial subtypes we identified here. Hence, our microglial scRNAseq data contribute further to the emerging single-cell landscape of this cell type.
Activated microglia are a hallmark in both brain cancers and temporal lobe epilepsy in response to the inflammation (Darmanis et al., 2017;Kaminska et al., 2021;Kinoshita & Koyama, 2021;Morin-Brureau et al., 2018;Ochocka et al., 2021). We acknowledge that the tissue used in this study is sourced from tumor or epileptogenic tissue adjacent regions in individuals with diagnoses of brain tumors or epilepsy, meaning it is not entirely healthy. Therefore, it is possible that the activated subtypes we identify here may be due to the nature of the tissue. Additionally, sex-specific expression differences have been observed in glioma-activated microglia . Nonetheless, we did not observe significant differences in most genes we highlighted here due to diagnosis ( Figure S5). There is also a large cluster of homeostatic microglia present in these samples, indicating that not all cells are activated. While the focus of our work is AD, these findings also have potential relevance in the fields of cancers, epilepsy, and other diseases.
We recognize that our study has several limitations, primarily owing to the difficulty in obtaining high-quality neurosurgical brain tissue, which leads to limited sample size and variability in tissue, diagnoses, and patient demographics. Even though we have utilized tissue surgically separated from disease tissue, the samples are from epilepsy and various brain tumor patients representing multiple diagnoses. We isolated microglia using an approach which should minimize activation; however, we cannot definitively rule out stress-induced transcriptomic changes during isolation. Despite these caveats, we could identify microglial co-expression modules and subclusters with multiple features that are consistent with prior publications from model systems (Hammond et al., 2019;Keren-Shaul et al., 2017;Krasemann et al., 2017;Marschallinger et al., 2020). Our scRNAseq clusters have contributions from both tumor and epilepsy samples, suggesting that our findings are unlikely to be driven by any one diagnoses. Furthermore, there are few published studies using fresh brain tissue to study microglia and thus our integrated single-cell data showing a homogenous cluster of microglia highlights the robustness of the methodology.
In summary, our study on sorted bulk microglia RNAseq and scRNAseq from fresh brain tissue yield several key findings. We identify a microglial gene signature from sorted bulk microglia, characterize its expression in bulk brain RNAseq across 7 datasets comprising 6 regions, in bulk fresh brain RNAseq and in microglial scRNAseq

| Patient samples
Fresh human brain tissue was obtained from patients undergoing neurosurgical procedures for epilepsy or tumor resection. Tissues determined to be grossly unaffected by the primary disease process were utilized for the present study ( Figure 1). Patient samples were transported from the operating room to the laboratory in 1X DPBS (Thermo Fisher; 14287080) for processing within 1-2 h of resection.
Human tissue was collected with informed consent prior to surgery and all procedures were approved by the Mayo Clinic Institutional Review Board and are HIPAA compliant.

| Tissue dissociation
Tissue was dissected to remove necrotic tissue, white matter and excess vascular tissue, to retain only cortical gray matter. The remaining tissue was cut into sagittal slices and weighed before being processed using the Adult Brain Dissociation Kit (Miltenyi; 130-107-677) as per the manufacturer's protocol. Debris removal (Miltenyi; 130-109-398) and red blood cell lysis (Miltenyi; 130-094-183) were also performed. All procedures were carried out on ice. The resulting homogenate was filtered through a 70µm filter before proceeding.

| Magnetic-activated cell sorting (MACS)
The cell suspension was first enriched for CD11b + cells by incubat-

| Identification of a core microglial signature from bulk microglia data
To define a core microglial signature, we compared our bulk microglia data to cognitively normal control samples from the AMP-AD bulk tissue transcriptome data from 7 different datasets representing 6 Genes that passed these criteria and were significant in comparisons with all 7 bulk brain datasets determined the microglial signature.
These signature genes were assessed for GO term enrichment with biological pathways using MsigDB. REViGO (Supek et al., 2011) tree plots were generated in R using GO terms obtained from MsigDB.

| Weighted gene co-expression network analysis
The CQN normalized expression values from bulk microglia were input to R WGCNA (Langfelder & Horvath, 2008)

| Over-representation and correlation analyses
Hypergeometric testing was performed in R to determine the enrichment of a select set of genes in previously reported signatures, bulk tissue expressed genes, WGCNA modules or 10X single-cell clusters. Correlation between bulk tissue and bulk microglial normalized CQN data was calculated using Spearman's rank correlation.
Concordant and discordantly correlated genes were determined using the upper and lower quartiles from each dataset. Quality control, normalization, clustering, and marker gene identification were performed with Seurat v3 (Stuart et al., 2019), followed by annotation of clusters using established cell type markers. We kept

| Meta-analysis of bulk microglia RNAseq datasets
We obtained sets of co-expression genes, that is, WGCNA modules from our own data and identified several modules whose eigengenes were correlated with age/sex/APOE. To check whether these correlations hold for external datasets, we downloaded raw reads in FASTQ format for three external bulk microglia datasets listed in Table S17. Reads were mapped to human genome build hg38 and were counted and normalized in the same fashion as described previously. For each WGCNA module, we identified the central genes (i.e., genes whose Pearson correlation with module eigengene > 0.75). Using these central genes, module eigengenes were calculated in all datasets-ours and three external ones. We correlated module eigengenes with traits. Finally, we performed meta-analysis to combine correlations from multiple datasets using metacor function (random effect model) in R meta package.

| Meta-analysis of single-cell RNAseq datasets
Raw single-cell RNAseq data from Olah et al. (2020) was downloaded through Synapse (syn21438358) and processed through Cell Ranger. Quality control was performed as described prior to integration with our samples using Seurat v4.0.4. Integration was performed to combine the datasets by individual rather than sample. A total of 9 samples and 15,819 cells were retained from the Olah data.
Cells were annotated to include dataset of origin. Hypergeometric tests were performed to determine the enrichment of our microglial signature within each integrated cluster. Additionally, we also looked at overlap between integrated cluster genes and marker genes from our original single-cell dataset. To identify where our original cells localized within the integrated dataset, we mapped the original cell IDs to the integrated clusters and calculated percentage overlap.

ACK N OWLED G M ENTS
The authors thank the patients and their families for their participa- The results published here are in whole or in part based on data ob- Eric Schadt from Mount Sinai School of Medicine.

CO N FLI C T O F I NTE R E S T
None.

AUTH O R CO NTR I B UTI O N S
TP and NET wrote the manuscript; NET and MA designed the study; TP, XW, and ZQ performed data analysis; JC consulted on statistical methods; TP, TPC, XW, YM, and RMA generated tables and figures; EM, CAG, SG, KC, RW, HGC, and AQH provided neurosurgical tissue samples; TP, TPC, LJLT, SJL, SL, FQTN, CCGH, KGM, and TN performed experimental procedures from blood and tissue samples. All authors read the manuscript and provided input and consultation. NET oversaw the study and provided direction, funding, and resources.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data in this manuscript are available via the AD Knowledge Portal