Transcriptional diversity across human and iPSC-derived microglia atlas
To provide a compressive repertoire of microglia diversity in postmortem human brains and human cell based experimental models we integrated single-cell and SRT dataset from multiple cohorts and brain regions. Table 1 provides a summary of the clinical, neuropathological and experimental data of these datasets.
We harmonized microglia from two in vitro models: The first experimental dataset involves Xenografted-microglia (xMGs) that represents seven distinct clusters of approximately ~26K single cell nuclei(11). In this experiment, wild-type and R47H human microglial progenitors were transplanted into 5x-hCSF1 mice, to investigate how this mutation affects lipid droplet accumulation, and the microglial response to Aβ plaques. After seven months, GFP-expressing TREM2 and isogenic TREM2-R47H xMGs were extracted from the brains of 5x-hCSF1 mice and subjected to scRNA-seq(11). Each cluster represents a unique microglial state for both TREM2 and TREM2-R47H xMGs(11) (Fig 1a, Table 1a). We also included in this atlas multimodal CRISPR-based genetic screens in human iPSC-derived microglia(12) (Fig 1a). Unsupervised clustering and UMAP dimensional reduction of iTF-Microglia revealed nine transcriptionally distinct clusters for ~18K nuclei(12). Our atlas also integrated ex vivo data. We included microglia from postmortem human microglia from parietal cortex from our previous work(13). We identified nine microglia transcriptional states for about 16K nuclei from 57 donors, including 15 carriers of pathogenic mutations in APP and PSEN1 (autosomal dominant AD; ADAD), 25 sAD non-carriers. In addition, we included three individuals who matched AD-neuropathological criteria but without clinical cognitive impairment at death (presymptomatic), 7 individuals who matched non-AD neurodegenerative pathologies criteria (Other), and 7 individuals who exhibited neither neurodegenerative pathology nor evidence of dementia (Control) (Fig 1a, Table 1b). The last dataset includes ~166K single microglia nuclei from aged postmortem brain samples of 427 individuals from the ROSMAP study(14). This includes 207 AD donors (133 early and 74 late) and 220 control individuals across six brain regions: prefrontal cortex (PFC), hippocampus, mid-temporal cortex, angular gyrus, entorhinal cortex, and thalamus. Here, we identified 12 distinct microglia transcriptional states (Fig 1a, Table 1c).
We selected clusters that were sufficiently represented in each of the datasets (Fig 1b), considering both the number of cells contained within each cluster and the extent to which these clusters fairly represented the diverse range of donors included in the study. This selection process ensured that each cluster provided a comprehensive overview of the biological diversity present across the samples. By maintaining broad biological representation, we were able to maximize the significance and robustness of the integration, thereby enhancing the reliability and validity of our findings (Additional file 1). After comparing methods and optimizing parameters, we chose the Canonical Correlation Analysis (CCA) method implemented in the SeuratV5 package that high accuracy and exhibited minimal batch effects (Fig S1a-h). We harmonized a total of 222,822 nuclei into 14 clusters that are represented in the four datasets (Fig 1c). We confirmed that all samples and batches were adequately represented across the different clusters (Fig S2a-b).
We examined the contribution of each cohort’s transcriptional states in these integrated clusters. Cells from all cohorts were reorganized while preserving their biological relevance. (Fig 1d, see NMI in Methods: Fig S1). For instance, the majority of homeostatic cells from both in vitro and ex vivo were grouped into Cluster 0, stress-related cells from both ex vivo were found in Cluster 6 while activated states from in vitro and the ex vivo parietal cortex cohort converged into Cluster 8. Although no activated cluster was clearly identified in the Human multiple region (ROSMAP) samples(14), by virtue of aggregating multiple dataset we were able to reassign few nuclei to the activated cluster (Fig 1d). Finally, to evaluate the association of transcriptional states to amyloid plaques and tau tangles, we selected SRT datasets from human Middle temporal gyrus (MTG) AD samples (N=3) generated using the 10X genomics Visium platform.
Characterization of integrated transcriptional states
We annotated the 14 integrated clusters based on their molecular signatures and functions, utilizing well-known canonical standard markers for Homeostatic, Activated, IFN, and MHCII (22)(Additional file 2). We expanded this list by incorporating additional representative markers from Sun et al. (2023)(14) for newly reported transcriptional states. We identified an average of 941 upregulated genes (FDR <0.05) for each transcriptional, revealing a rich transcriptional diversity within all integrated transcriptional states (Additional file 3).
Homeostatic microglia (Cluster 0) were characterized by high expression of well-known homeostatic markers such as P2RY12 (log2FC = 0.28, adj. p = 7.73 × 10-504), P2RY13 (log2FC = 0.21, adj. p = 5.24 × 10-131), and CX3CR1 (log2FC = 0.19, adj. p = 9.23 × 10-165) (Fig 1e, Additional file 3). Activated microglia (Cluster 8) exhibited higher expression of SPP1 (log2FC = 1.45, adj. p = 7.97 × 10-3290), CD9 (log2FC = 1.83, adj. p = 2.27 × 10-1331), and TREM2 (log2FC = 0.59, adj. p = 1.50 × 10-511) (Fig 1e, Additional file 3).
In MHCII (Cluster 3), we observed higher expression of HLA-A (log2FC = 0.22, adj. p = 1.52 × 10-182), HLA-E (log2FC = 0.08, adj. p = 4.57 × 10-42), and HLA-DPB1 (log2FC = 0.64, adj. p = 6.17 × 10-529) and marker genes associated with “positive regulation of lymphocyte proliferation” (adj. p = 1.75 × 10−2) (Fig 1e, Additional file 3 and 4). In MHCII (Cluster 3), we observed higher expression of HLA-A (log2FC = 0.22, adj. p = 1.52 × 10-182), HLA-E (log2FC = 0.08, adj. p = 4.57 × 10-42), and HLA-DPB1 (log2FC = 0.64, adj. p = 6.17 × 10-529) and marker genes associated with “positive regulation of lymphocyte proliferation” (adj. p = 1.75 × 10−2) (Fig 1e, Additional file 3 and 4).
We identified three inflammatory states showing high expression of IFN-induced gene (Clusters1, 5, and 7). Cluster 7 showed overexpression of IFIT3 (log2FC = 1.57, adj. p = 1.37× 10−1418) (Fig 1e, Additional file 3), and showed strong functional enrichment in immune responses, including defense response to virus (adj. p = 7.45 × 10−20), and cytokine-mediated signaling pathway (adj. p= 2.77× 10−11) (Additional file 4). We also found TMEM163 (log2FC = 0.52, adj. p = 9.82 × 10−483) in Cluster 1 and LRRK2 (log2FC = 0.36, adj. p = 1.45 × 10−125) in Cluster 5 (Fig 1f, Additional file 3).
We observed significant enrichment of stress markers such as HSP90AA1 (log2FC = 1.73, adj. p = 5.92 × 10-6228) and HSPH1 (log2FC = 1.56, adj. p = 3.21 × 10-3248) related to the regulation of cellular response to heat (adj. p = 3.23 × 10-8) in Cluster 6 (Fig 1f, Supplementary Tables 3 and 4). Interleukin IL1B (log2FC = 3.70, adj. p = 8.41 × 10-2568) and CCL4 (log2FC = 4.82, adj. p = 1.31 × 10-1480) and additional genes related to the “cytokine-mediated signaling” (adj. p = 2.29 × × 10-10) were overexpressed in Cluster 11 (Fig 1f, Supplementary Tables 3 and 4).
Cluster 4 showed expression of SPARC (log2FC = 0.51, adj. p = 1.27× 10−55) (Fig 1e, Additional file 3), and strong functional enrichment of extracellular structure organization (adj. p = 2.68 ×10-6), external encapsulating structure organization (adj. p = 2.68 ×10-6), extracellular matrix organization (adj. p = 5.64 ×10-4), and ECM-receptor interaction (adj. p = 8.17×10-4), henceforth designated as Extracellular Vesicular (EV) microglia (Additional file 4). And Cluster 2 exhibited high expression of neurotransmitter receptors related to cell-cell junction organization (adj. p = 2.73 × 10−3) (Fig 1f, Additional file 4), while we identified Cluster 12 as a cycling state due to the enrichment in DNA metabolic process (adj. p = 2.42 × 10−31) (Fig 1f, Additional file 4). We annotated Cluster 13 as a lipid-processing state based on the enrichment of representative marker TPRG1 (log2FC = 1.20, adj. p = 9.96 × 10−19) (Fig 1f, Additional file 3). We did not find any specific functional enrichment in clusters 9 and 10, neither were those nuclei flagged by our quality control process. Therefore, we labeled them as Unknown (Fig 1e-f, Additional file 3).
Furthermore, we examined the expression of each marker gene in each dataset (Fig 2a and Fig S3), to verify they were consistently expressed. Our analysis revealed consistency in gene expression patterns across all datasets, further enhancing the reliability and generalizability of our insights. In summary, Clusters 0, 3, 8, 6, and 11 were classified as Homeostatic, MHCII, Activated, Stress, and IL1B. Based on strong evidence of gene expression and proportional analysis, Clusters 0, 3, 8, 6, and 11 were classified as Homeostatic, MHCII, Activated, Stress, and IL1B. IFN-I, IFN-II, and IFN-III cells were identified in clusters 7, 5, and 1. Clusters 2, 4, 12, and 13 were annotated as Neuronal Surveillance, extracellular vesicular, Cycling, and Lipid Processing respectively (Fig 1 c). The biological pathway analysis of each annotated transcriptional state is shown in Additional file 5.
Relative proportions of microglia transcriptional states across AD groups
To evaluate whether the relative proportions of these transcriptional states change in AD, we tested the statistical significance of cell fractional differences between control and other individuals, including those in ex vivo datasets and in Xenografted-mic with AD risk variants in TREM2 and wildtype (Fig 2b, Additional file 6).
A prominent Stress cell state was enriched in ADAD samples from the Human parietal cortex (WASHU) cohort (β = 0.30, p = 1.67 × 10−2; Fig 2b, Additional file 6). The Gene Ontology analysis revealed its association with “regulation of cellular response to heat” (adj. p = 3.29 × 10-8), “chaperone cofactor-dependent protein refolding” (adj. p = 1.15 × 10-7) and “regulation of inclusion body assembly” (adj. p = 6.19× 10-4) (Fig 2c, Additional file 7). These nuclei overexpressed HSP90AA1 and HSPH1 genes, suggesting these play a role in synaptic homeostasis and Aβ pathology, as previously observed in the entorhinal cortex in AD(23). In the case of Human multiple regions (ROSMAP) we did not find significant enrichment of Stress transcriptional state in AD groups comparison to control.
Both ex vivo datasets exhibited enrichment of Activated cell states in AD groups, but such association was not significant in the Xenografted TREM2-R47H samples. The Activated transcriptional state was specific to ADAD (β = 0.23, p = 1.78 × 10−3; Fig 2b, Additional file 6), and sporadic AD (β = 0.19, p = 5.53 × 10−3; Fig 2b, Additional file 6) in the Human parietal cortex (WASHU), similarly in the earlyAD (β = 0.12, p = 9.44 × 10−12; Fig 2b, Additional file 6) and late AD (β = 0.19, p = 9.81 × 10−17; Fig 2b, Additional file 6) from the Human multiple region (ROSMAP) cohort. Several experimental models support that in response to Aβ pathology, microglia cells proliferate, and acquire an activated transcriptional profile(24–26). However, excessive microglia activation may exacerbate neuronal damage at later stages of amyloid pathology(24,27–29). A gene ontology analysis reveals that activated cells are associated with the regulation of cell migration (p = 1.67 × 10-4), regulation of ERK1 and ERK2 cascade (p = 1.18 × 10-2), and positive regulation of stress fiber assembly (p = 2.61 × 10-3) (Fig 2c, Additional file 7). Additional studies support that the ERK1/2 pathway governs inflammatory cytokine production and iNOS expression in activated microglia, highlighting its importance in neuroinflammatory processes and as a potential therapeutic target (30). Additionally, activation of the Rho/ROCK signaling pathway boosts stress fiber assembly in microglia, facilitating essential cytoskeletal rearrangements for phagocytosis and activation (31,32).
The Cycling transcriptional state was enriched in Presymptomatic individuals (general linear regression β = 0.19, p = 1.86 × 10−2; Fig 2b, Additional file 6) from the Human parietal cortex (WASHU). Biological pathway analysis revealed associations of the Cycling transcriptional state with the DNA metabolic process (p = 2.42 × 10-31), microtubule cytoskeleton organization involved in mitosis (p = 1.12 × 10-19), and mRNA splicing via spliceosome (p = 1.32 × 10-16) (Fig 2c, Additional file 7). These findings suggest a potential role for these pathways in driving microglial proliferation, implicating their involvement in the early stages of neuroinflammatory processes associated with AD.
Neuronal surveillance cells were enriched in late onset AD (general linear regression β = 0.06, p = 2.06 × 10−2; Fig 2b, Additional file 6) samples from multiple region cohorts, showing it significant involvement in “cell-cell adhesion via plasma-membrane adhesion molecules” (p = 9.00 × 10-4), “cell-cell junction organization” (GO:0045216) (p = 2.73 × 10-3), and “regulation of cell migration” (p = 4.52 × 10-2) (Fig 2c, Additional file 7). These processes suggest a potential role in modulating immune cell infiltration, microglial activation and possibly involvement in the blood-brain barrier integrity, highlighting them as potential therapeutic targets for mitigating neuroinflammation in AD(33–36).
Microglial transcriptional shift in response to AD pathology
We integrated SRT datasets from human MTG AD samples (N=3) generated using 10X genomics Visium platform(21). We selected those spots characterizing the cortical layers II to VI obtaining a total of 7,884 spots. As few spots were captured for layer I, we removed those from our analyses. Thus, we integrated this data with the scRNA-seq datasets(11–14) to investigate the distribution of microglia transcriptional states among cortical layers (Fig 3a) and their proximity to AD pathological lesions (Aβ or neurofibrillary tangles (NFT)). To achieve this, we conducted a reference-based integration to deconvolve and assign scores to each spot based on the resemblance of microglial transcriptional states. In our analysis, we focused on the highly predicted (top 25%) Visium spots for each AD pathology, examining the distribution of Homeostatic and Activated transcriptional states across cortical layers (Fig. 3b), and dataset (Fig S4).
We observed a significant increase in spots capturing Activated microglia in the external cortical layers II-III (p=2.20e-16), while internal layers IV-VI showed an increase in Homeostatic signatures (p<2.2e-16) (Fig. 3c, Additional file 8). We discarded that this was not driven by Aβ or NFT. However, in the case of the Human multiple regions (ROSMAP) reference, we were unable to capture significant Activated spots (p=5.24e-01).
Furthermore, we analyzed the relationship of Activated and Homeostatic transcriptional states in relation to the proximity of amyloid plaques (identified using Rabbit anti-Aβ antibody) and NFT (identified using Human/murine phospho-tau pSer202/Thr205 antibody) (Fig. 3d). We identified that spots proximal to plaques (<55 microns) tended to show Activated microglia (p=1.06e-08), while distal spots showed a predominance of Homeostatic microglial expression (p=1.36e-06). The presence of the microglial Activation state in proximity to amyloid plaques was predominantly observed in external cortical layers II and III (p=2.56e-03), while the Homeostatic state was predominantly observed distal to amyloid plaques in external cortical layers II and III (p=6.95e-03) (Fig. 3e, Additional file 9). Notably, no such distinction was observed in relation to NFT (Fig S5). To ensure that our observations were not driven by a single AD sample and/or limited to integrated data, we examined the distribution of transcriptional states across all AD samples using each cohort as a reference (Fig S6, Additional file 9). We found a prevalence of microglial Activation states near amyloid plaques, and the Homeostatic state was predominantly observed distal to amyloid plaques in each sample. Importantly, this distribution pattern was not driven by differences in the burden of amyloid plaques among cortical layers.
Validation of homeostatic microglia distributions in the cortical layers of human MTG
We performed immunochemistry experiments to validate the finding from the transcriptional studies. We stained for P2RY12, a unique marker for homeostatic microglia to map and analyze the distribution and activity of microglia in the adjacent 10 µm section of human MTG, which is 10 µm interval from the Visium gene expression section(21). This technique leverages antibodies that specifically bind to the P2RY12 protein, a purinergic receptor predominantly expressed in homeostatic microglia. We performed P2RY12 and Aβ co-immunostaining to assess the spatial relationship of homeostatic microglia with Aβ plaques, although the immunostaining of P2RY12 and Aβ plaque was not on the same section with Visium gene expression assay.
Our registration of the immunostaining image with H&E staining on the Visium gene expression sections shows that the distribution of P2RY12-positive (+) counts (p = 1.54e-04) and the percentage of the area with P2RY12+ signals (p = 4.41e-02) across the internal cortical layers IV-VI, as observed through single-cell and SRT data (Fig. 4b, Additional file 10). The distribution of Aβ+ spots (p = 1.06e-01), as well as the overlap of P2RY12+ with Aβ+ spots (p = 1.95e-03) among cortical layers, further confirmed that this distribution was not driven by amyloid plaques (Fig. 4b, Additional file 10).