Molecular and genetic landscapes of retina and brain microglia in neurodegenerative diseases

Retina and brain microglia share transcriptional signatures

We analyzed data from 97 independent human samples from retina (Fig. 1A) across 11 different studies and analyzed their transcriptome (single-nucleus RNA sequencing [snRNA-seq]) (Fig. 1B) and epigenome (single-nucleus ATAC sequencing [snATAC-seq]) (Fig. 1C). Finally, we explored the disease implications of retinal microglia in AD and brain microglia in AMD using respective GWAS data from genetic studies (Fig. 1D; Fritsche et al. 2016; Bellenguez et al. 2022).

Figure 1.

Overview of the study. (A) Retina and brain microglia share a common origin from the yolk sac. (B,C) Integration of single-nucleus RNA sequencing (snRNA-seq) and single-nucleus ATAC sequencing (snATAC-seq) data from 97 human retina and brain microglia samples. (D) Integration of genome-wide association study (GWAS) data for age-related macular degeneration (AMD) and Alzheimer's disease (AD) with transcriptomic and epigenetic annotations using various enrichment methods.

Open in new tabDownload PowerPointLink to figure

Microglia constitute ∼0.2%–0.4% of total retinal cells (Li et al. 2019; Wang et al. 2023). To ensure an adequate number of cells for analysis, we integrated the snRNA-seq of 28 human retinas from five different data sets (Supplemental Table S1; Orozco et al. 2020; Wang et al. 2022; Kuchroo et al. 2023; Yu et al. 2024) and implemented fast integration pipeline in Seurat 4.4 (Stuart et al. 2019) to harmonize the data (Supplemental Fig. S1). In total, we analyzed 214,635 cells and identified 11 cell types using cell type–specific markers (Supplemental Table S2; Macosko et al. 2015; Yan et al. 2020), including amacrine, astrocyte, cones, horizontal, microglia, Müller glia, OFF cone bipolar, ON cone bipolar, retinal ganglion, rod bipolar, and rod cells (Supplemental Table S3; Supplemental Fig. S2). For the brain, we analyzed 125,343 cells from three data sets (Velmeshev et al. 2023; Haney et al. 2024; Sepp et al. 2024), comprising samples from 24 individuals (Supplemental Table S1). We annotated six distinct brain cell types, including astrocytes, endothelial cells, excitatory neurons, inhibitory neurons, microglia, and oligodendrocytes using cell type–specific markers from published studies (Supplemental Tables S2, S3; Supplemental Fig. S3; Haney et al. 2024). Next, we integrated and clustered cells from the retina and brain, applying batch effect correction, and visualized them using uniform manifold approximation and projection (UMAP) plots generated with the Seurat package (Hafemeister and Satija 2019). Our analysis revealed overlapping microglia cell clusters between the retina and brain, indicating similar transcriptional signatures (Fig. 2A). To identify molecular similarities among different cell types in the retina and brain, we performed a cross-tissue analysis for both integrated data sets, using a tree-based method, SpeciesTree (Bakken et al. 2021). The resulting hierarchical clustering revealed that retina and brain microglia were closely related compared with other cell types (Fig. 2B).

Figure 2.

Shared transcriptome signature of retina and brain microglia. (A) The snRNA-seq analysis of retina and brain. The plot shows uniform manifold approximation and projection layouts (UMAP) for the 11 retina cell types and six brain cell types, with overlapping microglia cluster highlighted in a red circle. (B) Hierarchical clustering of cell type clusters from the retina and brain reveals distinct yet overlapping patterns of similarity. Notably, retina and brain microglia cluster together on the same branch, indicating a close relationship between these cell types across both tissues. (C) Dot plot of module scores displaying the z-scores for transcriptional signatures based on differentially expressed genes and curated markers, visualizing the expression profiles of annotated microglia subclusters. The color and size of each dot correspond to the average expression level and fraction of expressing cells, respectively. (D) High-resolution UMAP plots showing the distinct clustering of microglia, PVMs, and monocytes in the retina. (E) High-resolution UMAP plots of microglia, PVMs, and monocytes in the brain. (F) The UMAP plots show that the clusters of microglia, PVMs, and monocytes from both the retina and brain overlap perfectly, indicating a high degree of similarity in their subcluster organization across these two tissues. (RGC) Retinal ganglion cell, (hMG1) human retinal microglia clusters marker 1, (hMG2) human retinal microglia clusters marker 2 (hMG3) human retinal microglia clusters marker 3, and (PVM) perivascular macrophages.

Open in new tabDownload PowerPointLink to figure

Microglia and monocytes are both members of the mononuclear phagocyte system and share essential tissue-specific functions that are critical for homeostasis and the response against pathogen evasion (Guilliams et al. 2014). Thus, we further refined the microglia cluster to distinguish yolk sac–derived microglia (Ginhoux et al. 2010) from perivascular macrophages (PVMs) and monocyte-derived myeloid precursors (Guilliams et al. 2014). We annotated these cell types based on the module scores (Fig. 2C) calculated using the top differentially expressed genes identified in human donor retina data from retina (Yu et al. 2024) and brain microglia (Böttcher et al. 2019) as well as cell type–specific marker genes (Supplemental Table S2; Butovsky and Weiner 2018; O'Koren et al. 2019; Hao et al. 2021; Bloomfield et al. 2024; Yu et al. 2024) for the retina (Fig. 2D) and brain (Fig. 2E). Upon merging, 2064 retina microglia cells clustered with 4927 cells from brain microglia, forming a distinct and unified microglial group. Transcriptomic signatures from these microglia subcluster in the retina and brain were used in all subsequent analyses to ensure greater specificity and resolution. Additionally, retina PVMs clustered with brain PVMs, as well as retina monocytes with brain monocytes, each forming separate and distinguishable groups (Fig. 2F). Taken together, our analyses shows that the transcriptome profile of retina microglia closely resembles that of brain microglia.

Chromatin accessibility is similar in retina and brain microglia

RNA expression is regulated by the changes in chromatin accessibility at gene promoters and other distal regulatory elements, such as enhancers. Furthermore, snRNA may not fully capture a cell's transcriptional profile, as it primarily detects mRNA levels. Thus, to understand a comprehensive picture of transcriptional regulation that dictates microglia identity, we compared global chromatin accessibility profiles of the retina and brain using the assay for transposase-accessible chromatin using sequencing (ATAC-seq) (Buenrostro et al. 2013). We integrated chromatin accessibility data from 38 human retinas from three different studies (Orozco et al. 2020; Wang et al. 2022, 2023) and eight samples from three human brain data sets (Supplemental Table S1; Supplemental Fig. S4; Corces et al. 2020; Morabito et al. 2021; Zhu et al. 2023). After performing initial analyses with SEURAT (Hafemeister and Satija 2019; Stuart et al. 2021), we annotated 11 cell types in the retina (Supplemental Table S2; Supplemental Fig. S5) and six cell types in the brain (Supplemental Table S2; Supplemental Fig. S6), using cell type–specific markers for the retina (Macosko et al. 2015; Yan et al. 2020) and brain (Haney et al. 2024). This analysis identified 275,579 and 366,041 peaks in the retina and brain, respectively. Next, we integrated the data by clustering 255,747 retina and brain cells (Supplemental Table S3), correcting for batch effects, and visualized the results with UMAP plots generated using Seurat (Hafemeister and Satija 2019). This analysis revealed overlapping microglia cell clusters between the retina and brain, indicating shared epigenetic signatures (Fig. 3A). Hierarchical clustering further supported this, with microglia from both retina and brain clustering in a distinct branch separate from other cell types in these tissues (Fig. 3B).

Figure 3.

snATAC-seq analysis of the retina and brain. (A) UMAP of annotated clusters based on chromatin accessibility data from retina and brain cell types, with a red circle highlighting the overlap between retina and brain microglia clusters. (B) Hierarchical clustering of chromatin accessibility data shows that retina and brain microglia form a separate, closely related cluster, highlighted in a red box, distinct from other cell types. (C) Module scores based on differentially expressed gene lists and our curated marker list are visualized using a dot plot for the annotated microglia, PVMs, and monocyte clusters. The color and size of each dot correspond to the average expression level and fraction of expressing cells, respectively. (D) High-resolution UMAP plots of chromatin accessibility data highlight the distinct clustering of microglia, PVMs, and monocytes within the retina. (E) Higher-resolution UMAP plots of brain ATAC-seq data shows distinct clustering of microglia, PVMs, and monocytes. (F) UMAP plots reveal that the chromatin accessibility profiles of microglia and PVMs overlap closely between the retina and brain, suggesting conserved subcluster organization across these tissues. (RGC) Retinal ganglion cell, (hMG1) human retinal microglia clusters marker 1, (hMG2) human retinal microglia clusters marker 2, (hMG3) human retinal microglia clusters marker 3, and (PVM) perivascular macrophages.

Open in new tabDownload PowerPointLink to figure

Similar to transcriptome, we next performed unsupervised clustering at higher resolution followed by annotation based on module scores derived from retina (Yu et al. 2024) and brain microglia (Böttcher et al. 2019), as well as cell type–specific marker genes (Fig. 3C; Supplemental Table S3; Butovsky and Weiner 2018; O'Koren et al. 2019; Hao et al. 2021; Bloomfield et al. 2024; Yu et al. 2024). In our analysis of the retina, we were able to identify distinct immune cell populations, including microglia, monocytes, and PVMs as distinct clusters (Fig. 3D). A comparable assessment in the brain revealed the presence of these immune cell populations (Fig. 3E). A combined projection of the two tissues onto the UMAP showed significant overlap of 702 retina microglia with the 1887 brain microglia, with a clear distinction from the other cell types (Fig. 3F), confirming the similarity seen in the transcriptome level. For all subsequent analyses, microglia were referred to the annotations derived from these subclusters. We also observed that unlike in the transcriptomic data, in which monocytes from the retina and brain clustered together (Fig. 3F), in the ATAC-seq data, they were distinctly separated. This separation was further supported by differential analysis, which identified 99 differentially expressed genes between retina and brain monocytes compared with 1421 differentially accessible regions (DARs) in the snATAC-seq data (Supplemental Fig. S7). Notably, in the transcriptomic data set, the number of retina and brain samples was comparable, with 28 and 24, respectively. However, in the ATAC-seq data, we had a much larger number of from the retina—38 in total—compared with only eight samples from the brain. Therefore, the observed separation of monocytes in the ATAC-seq data may be influenced by sample imbalance or may reflect tissue-specific epigenetic programming in response to the distinct local environments of the retina and brain. These findings demonstrate that the gene expression and chromatin accessibility of the retina and brain microglia cells are similar. Furthermore, despite their functional similarities, macrophages and microglia have distinct expression profiles, indicating cell type–specific functions.

Molecular signatures of retina and brain microglia

To gain deeper insight into the similarities and differences between microglia in the retina and brain at the molecular level, we examined global gene expression profiles using correlation heatmaps. For these analyses, we obtained the average gene expression for all cell types in the retina and brain and computed the Spearman's rank correlation for all possible pairs. Correlation between microglia, monocytes, and PVMs was significantly higher between the brain and retina compared with other cell types (Fig. 4A; Supplemental Fig. S8). We observed a clear separation between this cluster in which correlations within microglia and related-cell types were significantly higher than correlations between all other cell types of the retina and brain regardless of the tissues (P = 1.06 × 10⁻¹⁰). Specifically, the correlation of retina and brain microglia reached the value of 0.88 (Fig. 4B; Supplemental Fig. S8). For the ATAC-seq data, we observed the same pattern with brain and retina microglia, monocytes, and PVMs, showing higher (Fig. 4C) and significantly different correlations compared with other cell types (P = 3.5 × 10⁻¹⁶). The correlation between the brain and retina microglia was 0.96 (Fig. 4D).

Figure 4.

Molecular similarities and differences in retina and brain microglia. (A) Correlation heatmap of gene expression between retina and brain cell types from snRNA-seq data. The color scale ranges from blue (low correlation) to red (high correlation). Column labels correspond to retina cell types, and row labels represent brain cell types. (B) XY plot depicting the log₂-transformed gene expression profiles of microglia from the retina (y-axis) and brain (x-axis). The correlation coefficient, displayed at the top of the plot. (C) Heatmap of correlations from the inferred gene activity matrix in snATAC-seq data, with a color scale from blue (low) to red (high). (D) XY plot of log₂-inferred gene activity from ATAC-seq chromatin accessibility profiles for the retina (y-axis) and brain (x-axis) microglia, with the correlation coefficient at the top denoting their similarity. (E) Venn diagram illustrating the overlap between 596 retina microglia–specific genes (compared with other retina cell types) and 470 brain microglia–specific genes (compared with other brain cell types). (F) Venn diagram illustrating the overlap between 596 retina microglia–specific genes (vs. other retina cell types) and 374 genes upregulated in retina microglia compared with brain microglia (differentially expressed genes between brain and retina microglia). (G) Venn diagram showing the overlap between 470 brain microglia–specific genes (vs. other brain cell types) and 547 genes upregulated in brain microglia relative to retina microglia (differentially expressed genes between brain and retina microglia). (H) Dot plot showing scaled gene expression of the top 20 genes upregulated in retina microglia compared with other retina cell types and brain microglia. (I) Dot plot of the top 20 genes upregulated in brain microglia compared with other brain cell types and retina microglia. (J) Dot plot displaying marker genes for microglia, with the top panel representing retina microglia and the bottom panel representing brain microglia. (K) Bar plot illustrating the number of shared and unique ligands and receptors between retina and brain microglia. (RGC) Retinal ganglion cell.

Open in new tabDownload PowerPointLink to figure

A differential expression (DE) analysis of gene expression profiles identified 596 genes that were upregulated (more than twofold) in retina microglia compared with other retinal cell types (Fig. 4E; Supplemental Table S4). Similarly in the brain, 470 genes were upregulated in microglia compared with other brain cell types (Fig. 4E; Supplemental Table S4). One hundred ninety-seven of these genes were shared between the retina and brain microglia (Fig. 4E). Next, we performed DE analyses between retina and brain microglia that identified 374 genes upregulated in retina microglia and 547 in brain microglia (Supplemental Table S4). We then compared these two gene lists, identifying 311 genes upregulated in retina microglia relative to other cell types in the retina and brain microglia (Fig. 4F). In the brain microglia, we identified 218 genes that were upregulated compared with other brain cell types and retina microglia (Fig. 4G). Among the top 20 differentially regulated genes (Fig. 4H,I), many were immune-related and included several known microglia marker genes (Supplemental Fig. S9). In the retina microglia, we identified upregulation of CD74, C1QA, and MHC class II molecules (HLA-DRA, HLA-DPA1, and HLA-DPB1) whereas P2RY12 and CX3CR1 were upregulated in brain microglia (Fig. 4J; Supplemental Table S4).

Considering the dynamic role of transcription factors (TFs) in determining transcriptional program and regulation, we next compared the TF motifs of retina and brain microglia using Signac (Stuart et al. 2021). Among these, 140 motifs were shared between the two tissues, each showing greater than twofold enrichment (Supplemental Table S5). These included motifs for TFs previously associated with microglial function, such as HINFP, EGR2, HES1, and PLAG2 (Manouchehri et al. 2023). These findings support our characterization of both the similarities and differences between retina and brain microglia, with shared motifs reflecting substantial commonalities. Finally, we compared the intracellular communications between retinal and brain microglia using LRLoop (Xin et al. 2022), which utilizes transcriptome data from both ligand sender and receiver cells. In the retina, we identified 110 ligands and 102 receptors specific to microglia compared with other retinal cell types. In the brain, microglia exhibited 74 ligands and 88 receptors (Supplemental Table S6). Of these, 55 ligands and 57 receptors were shared between retinal and brain microglia (Fig. 4K). These shared interactions were enriched in Gene Ontology (GO) terms related to immune regulation and phagocytosis (Supplemental Table S6). Taken together, our data identify unique signatures in retina and brain microglia and highlight the molecular diversity despite the common origin.

Retina and brain microglia are dominant pathogenic cell type for AMD and AD

As our understanding of immune system dysfunction in relation to AMD (Nussenblatt et al. 2014) and AD (Bettcher et al. 2021) pathology grows, the role of microglia in the pathogenesis of AMD (Yu et al. 2024) and AD (Hansen et al. 2018) has garnered considerable attention. However, the molecular similarities between microglia's roles in AMD and AD are still not fully elucidated. Therefore, we set out to examine the molecular-genetic mechanisms of AMD and AD in the context of retina and brain microglia. We examined the relationship between the microglia-enriched/specific gene expression signatures and the disease risk for AMD and AD. We generated the microglia-enriched/specific signature using the “FindMarkers” (Butler et al. 2018) function, which compared the microglia cluster against all other clusters to identify potential marker genes (Supplemental Table S4). We also included three gene sets: 197 genes common to both retina and brain microglia (Fig. 4E), as well as genes upregulated in either retina or brain microglia compared with the other, with 374 upregulated in retina and 547 in brain microglia (Fig. 4F,G).

Next, we employed two enrichment methods to evaluate whether cell type–specific or enriched gene sets from the retina and brain show significant association with AMD and AD. First, we used multimarker analysis of GenoMic annotation (MAGMA) (de Leeuw et al. 2015), which uses a multiple regression approach to properly incorporate LD between markers and to detect multimarker effects. For this analysis, we used summary statistics from two GWAS studies: one with 16,144 AMD patients and 17,832 controls, identifying 34 AMD loci (Fritsche et al. 2016), and another with 111,326 AD cases and 677,663 controls, identifying 75 AD loci (Bellenguez et al. 2022). Using the MAGMA baseline model, we examined the effects of variant-to-gene mapping augmentations on gene-set scores and identified enrichment for retina and brain microglia in AMD (Fig. 5A; Supplemental Table S7). Similarly, the gene sets in both retina microglia and brain microglia were also significantly enriched in AD (Fig. 5B). Next, we implemented a LD score regression method, LDSC-SEG (Finucane et al. 2018), to quantify the heritability enrichment (Supplemental Table S8). Consistent with the results of MAGMA, LDSC-SEG identified retina and brain microglia to be significantly enrichment for both AMD (Fig. 5C) and AD (Fig. 5D). In addition, the 197 genes that were shared between retina and brain microglia also showed significant enrichment for AMD and AD in MAGMA and for AD in LDSC-SEG (Fig. 5A–D).

Figure 5.

Integration of gene expression and GWAS data identified microglia as disease-critical cell type in AMD and AD. (A) MAGMA enrichment analysis of gene expression data from retina and brain cell types for AMD-associated variants. The plot shows the P-value for MAGMA enrichment of AMD-associated variants in gene expression across retina cell types. The x-axis represents the mean strength of association (–log₁₀(P-value)), and the y-axis displays the different cell types. (B) Enrichment of AD–GWAS variants within retina and brain cell types using MAGMA. (C) Retina and brain cell type–specific enrichment of SNP heritability for AMD estimated using LDSC-SEG. The y-axis displays negative log₁₀ P-values of enrichment P-value for each individual test. (D) Enrichment estimates for retina and brain cell types for AD–GWAS variants using LDSC-SEG. The dashed line indicates the P-value cutoff for individual annotations, and the asterisks represent statistical significance at a false-discovery rate (FDR) of >5% across all four plots. (E,F) Heatmap of AMD–GWAS genes across retina and brain cell types. This highlights that AMD–GWAS genes that are enriched in retina microglia (E) are also predominantly expressed in brain microglia (F). (G,H) Heatmap of AD–GWAS genes across retina and brain cell types. This reveals that AD–GWAS genes enriched in brain microglia (H) are also highly and specifically expressed in retina microglia (G). (RGC) Retinal ganglion cell, (Common) genes common to both retina and brain microglia, (Retina-specific) genes upregulated in retina microglia compared with brain microglia, and (Brain-specific) genes upregulated in brain microglia compared with retina microglia.

Open in new tabDownload PowerPointLink to figure

Next, we investigated how GWAS signals are enriched within genes differentially expressed between brain and retinal microglia. For AMD, signatures from both the retina and brain showed enrichment, whereas for brain microglia, enrichment was observed only with AD in the MAGMA analyses (Fig. 5A,B; Supplemental Table S7). In LDSC-SEG analyses (Supplemental Table S8), retina-specific and brain microglia–specific signatures showed enrichment exclusively for AMD and AD, respectively (Fig. 5C,D). These findings indicate that although microglia in both the retina and brain play roles in AMD and AD, the distinct molecular signatures are likely specific to each disease.

Finally, we examined the expression of genes (40 for AMD [Fritsche et al. 2016] and 76 for AD [Bellenguez et al. 2022]) located within GWAS-associated loci (Supplemental Fig. S10). We selected genes with normalized expression values greater than 1.5 in retina or brain microglia. Among these, several AMD-associated genes showed selective or preferential expression in retina microglia (Fig. 5E), whereas those linked to AD exhibited similar patterns in brain microglia (Fig. 5H; Hansen et al. 2018). We observed that AMD genes that were highly expressed in retina microglia also showed specific expression in brain microglia (Fig. 5F), as well as AD genes in retina microglia (Fig. 5G), consistent with shared transcriptome profiles of retina and brain microglia. Collectively, these results implicate microglial dysfunction as a contributing factor in AMD and AD pathogenesis.

Retina and brain microglia open chromatin is enriched in disease heritability

Because majority of AMD- and AD-associated variants are noncoding (Fritsche et al. 2016; Bellenguez et al. 2022), they are likely to impact the function of regulatory elements. Thus, we leveraged snATAC-seq data to determine the relevance of retina and brain microglia in AMD and AD using MAGMA (de Leeuw et al. 2015). First, we identified the DARs for each cell type within the retina and brain, using a log-fold-change threshold of 0.25 and Bonferroni-adjusted P-values (P_adj < 0.05) (Supplemental Table S9). DARs represent the peaks that were either exclusively or more open in one cell type compared with others. The DARs were then assigned to the nearest genes using Signac (Stuart et al. 2021), and gene-set enrichment analysis was conducted using MAGMA (de Leeuw et al. 2015) to determine whether these genes were enriched in genetic associations with AMD and AD. We observed enrichment of AMD–GWAS variants in retina microglia compared with other retina cell types (Fig. 6A; Supplemental Table S10). Although this did not reach significance after correction for multiple testing, likely owing to the relatively small sample size of the AMD–GWAS data set. Among brain cell types, microglia and excitatory neurons showed significant enrichment (Fig. 6A). For AD–GWAS, significant enrichment was observed only in microglia from both the retina (Fig. 6B) and brain (Fig. 6B). We next applied S-LDSC (Finucane et al. 2015) for partitioning heritability analyses of chromatin accessibility. We estimated heritability enrichment (defined as the proportion of heritability explained by SNPs in an annotation divided by the proportion of SNPs in the annotation) for each disease–cell type pair conditional on annotations from the baseline-LD model for each. For AMD, enrichment was not significant for retina or brain microglia (Fig. 6C), possibly because of small sample size in the AMD–GWAS data. For AD, both the retina (enrichment score = 31.00, SNP heritability explained = 27.06%) and brain microglia (enrichment score = 25.91, SNP heritability explained = 36.33%) showed significant enrichment (Fig. 6D; Supplemental Table S11). Next, we investigated how GWAS signals are enriched within DARs identified between brain and retinal microglia. In contrast to gene expression analyses, the DAR analysis revealed only a limited number of regions: 58 upregulated in retinal microglia and 45 in brain microglia. Subsequently, no significant enrichment was observed for AMD or AD, likely owing to limited statistical power (Fig. 6A–D; Supplemental Table S11).

Figure 6.

Integration of epigenetic annotations and GWAS shows enrichment for disease heritability in retina and brain microglia. (A) MAGMA enrichment analysis of open chromatin data. Bar plots display the mean strength of association (–log₁₀ P-values) for MAGMA enrichment in retina and brain cell types, highlighting significant enrichment in retina and brain microglia for AMD–GWAS variants. (B) MAGMA enrichment analysis of AD–GWAS variants in open chromatin data across retina and brain cell types, showing significant enrichment specifically in retina and brain microglia. (C) Heritability enrichment estimates from retina and brain annotations for AMD–GWAS variants using S-LDSC. (D) Enrichment estimates with S-LDSC for epigenetic annotations from retina and brain cell types for AD. In all four plots, the dashed line marks the P-value significance for individual annotations, and asterisks represent the significant at an FDR <5%. (E,F) GARFIELD enrichment analysis of open chromatin data for AMD–GWAS (E) and AD–GWAS (F) for retina and brain cell types. Radial lines show odds ratio (OR) values at five GWAS P-value thresholds (T) for ATAC-seq peaks from retina and brain microglia cell types. Colored dots in the inner ring of the outer circle show where GARFIELD enrichment is significant, with stronger significance closer to the center (from T < 10⁻⁵ to T < 10⁻⁸). The dot's color represents a specific cell type. (G,H) The bar plot shows enrichment for AMD-associated (G) and AD-associated (H) variants in ATAC-seq data from ENCODE and Roadmap Epigenomics data sets across different tissues. ORs for enrichment are shown for variants at GWAS threshold of P < 1 × 10⁻⁸ after multiple-testing correction for the number of effective annotations. The dashed lines represent the P-value significance threshold for enrichment. (RGC) Retinal ganglion cell, (Retina-specific MG) DARs identified between in retina microglia compared with brain microglia, and (Brain-specific MG) DARs identified between in brain microglia compared with retina microglia.

Open in new tabDownload PowerPointLink to figure

Finally, we applied GWAS analysis of regulatory or functional information enrichment with LD correction (GARFIELD) (Iotchkova et al. 2019) to estimate the enrichment of AMD–GWAS and AD–GWAS associations in chromatin accessibility peaks of 11 retina and six brain cell types (Supplemental Table S12). Significant enrichment of AMD-associated loci was observed in retina (OR = 2.68, P = 5.54 × 10⁻⁶) (Fig. 6E). Within the brain cell types, AMD–GWAS signals were most enriched in microglia (OR = 2.32, P = 2.91 × 10⁻⁵) compared with other cell types (Fig. 6E). Similarly, AD–GWAS signals were significantly enriched in both brain (OR = 2.53, P = 1.42 × 10⁻²⁰) and retina microglia (OR = 2.67, P = 9.23 × 10⁻¹⁹) at different GWAS P-value thresholds (Fig. 6F). The findings reveal enrichment of both AMD–GWAS and AD–GWAS signals in microglia, consistent with the shared profiles. To better understand the specificity of AMD–GWAS variants in brain microglia and AD–GWAS variants in retinal microglia, we aimed to extend our analysis by incorporating the regulatory annotation of open chromatin regions, as assessed through ATAC-seq, across 63 tissues available in the Encyclopedia of DNA Elements (ENCODE) (The ENCODE Project Consortium et al. 2020). The ENCODE data exhibited wide variation in sample sizes across tissues, leading to differences in the number of open chromatin regions, which could introduce bias in comparative enrichment analyses. To minimize this, we focused on regions that were uniquely open or shared with no more than 10% of the tissues. Along with microglia, AMD–GWAS variants were significantly enriched in the liver (Fig. 6G; Supplemental Table S12), consistent with the liver's role in complement pathways (Qin and Gao 2006), which accounts for significant proportion of AMD heritability (Fritsche et al. 2014). Similar analysis of AD–GWAS variants identified the most significant enrichment retina microglia (Fig. 6H; Supplemental Table S12). These findings emphasize the specificity of brain and retinal microglia compared with other tissues in explaining the heritability of AMD and AD, underscoring their critical relevance to both conditions.

Data acquisition and primary analysis

Data were downloaded from multiple published NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) accession numbers along with one data set on Zenodo (https://zenodo.org) and two data sets from the CZI public website (https://chanzuckerberg.com) (Supplemental Table S1). This study was performed in accordance with the ethical standards of the Declaration of Helsinki, and informed consent was taken from all the participants in original studies. The study was approved by the Baylor College of Medicine review board. FASTQ files were processed using Cell Ranger version 8.0.0 for RNA-seq data and Cell Ranger ATAC version 1.2.0 for ATAC-seq data. Multiomic single-nucleus data were processed using Cell Ranger Arc 2.0.2. Other data sets from Zenodo and CZI included preprocessed count matrices, which were directly loaded without the need for further primary analysis. GWAS summary statistics were obtained from the GWAS catalog (Supplemental Table S1).

snRNA-seq analyses

Single-nuclei data from 28 retina samples from five data sets (Orozco et al. 2020; Wang et al. 2022; Kuchroo et al. 2023; Yu et al. 2024) and 24 brain samples from three data sets (Velmeshev et al. 2023; Haney et al. 2024; Sepp et al. 2024) were integrated for the analysis (Supplemental Table S1). Counts data from Cell Ranger was analyzed using Seurat (Butler et al. 2018). Quality control was performed on each sample to remove low-quality cells and doublets with counts of fewer than 200 and more than 7500, respectively. Cells with mitochondrial reads >10% were also removed. We implemented the fast integration pipeline in Seurat 4.4 (Stuart et al. 2019) for performing normalization with dimension reduction using SCTransform V2 (Hafemeister and Satija 2019) and RunPCA(seed.use = 1) on data set–level objects. The lists for retina and brain data sets were fed into SelectIntegrationFeatures (nFeatures = 3000), PrepSCTIntegration, and FindIntegrationAnchors (normalization method = “SCT,” dims = 1:30, reduction = “rpca,” k.anchors = 20). Using the anchors, we ran IntegrateData (normalization method = “SCT,” dims = 1:30). Cells embeddings were regenerated again after integration and used for RunUMAP (reduction = “pca,” dims = 1:30, min.dist = 0.2, seed.use = 1). Unsupervised clustering and cluster identification were performed using FindNeighbors and FindClusters (algorithm = 1, resolution = 0.5, random.seed = 1). We annotated the cells using Seurat's addModuleScore to calculate scores for cell type–specific marker genes (Menon et al. 2019; Peng et al. 2019; Yan et al. 2020). Module scores were visualized using the DotPlot function from Seurat. Cells that could not be confidently assigned to specific cell types were excluded in downstream analysis. To combine data sets for the retina and brain, the two Seurat objects were merged and regenerated the cell embeddings using RunPCA, RunHarmony (group.by.vars = “sample”), and RunUMAP (reduction = “harmony,” dims = 1:30, seed.use = 1, min.dist = 0.2) (Korsunsky et al. 2019).

Subclustering analyses of snRNA-seq

We selected the clusters with the highest module score for microglia and regenerated the cell embeddings and UMAP plot. We then run unsupervised clustering and cluster identification and annotated them, including the sets of markers for microglia, PVMs, and monocytes (Supplemental Table S2; Butovsky and Weiner 2018; O'Koren et al. 2019; Hao et al. 2021; Bloomfield et al. 2024; Yu et al. 2024). We generated module scores based on two sets of lists of marker genes. The first list contains sets of markers for microglia, PVMs, and monocytes (Butovsky and Weiner 2018; O'Koren et al. 2019; Hao et al. 2021; Bloomfield et al. 2024; Yu et al. 2024). The second list includes the top differentially expressed genes, ranked by fold change, for hMG1, hMG2, and hMG3—corresponding to human microglia clusters 1, 2, and 3—as well as for PVMs and monocytes from a recently published single-cell atlas of the human retina (Yu et al. 2024) and brain (Böttcher et al. 2019). After annotation, cells identified as microglia, PVMs, and monocytes were extracted in the retina and brain for performing another round of integration (Stuart et al. 2019). The lists for retina and brain data sets were fed into SelectIntegrationFeatures (nFeatures = 3000), PrepSCTIntegration, and FindIntegrationAnchors (normalization method = “SCT,” dims = 1:30, reduction = “rpca,” k.anchors = 20). Using the anchors, we ran IntegrateData (normalization method = “SCT,” dims = 1:30). Cell embeddings were regenerated again after integration and used for RunUMAP (reduction = “pca,” dims = 1:30, min.dist = 0.2, seed.use = 1) (Stuart et al. 2019).

DE analysis was performed using FindMarkers (logfc_threshold = 0.25, min.pct = 0.2, test.use = “wilcox,” assay = “RNA”) from Seurat (Butler et al. 2018). Pairwise DE analysis was done for all annotated cell types to produce lists of DE genes for use in downstream processes. We also performed a separate DE analysis between the identified microglia from the retina and microglia from the brain. DE genes were filtered for genes with twofold change or higher.

snATAC-seq analysis for the brain and retina

We integrated chromatin accessibility data from 38 human retinas from three different studies (GEO: GSE196235, GSE247157) (Orozco et al. 2020; Wang et al. 2022, 2023) and eight samples from three human brain data sets (GEO: GSE174367, GSE204682, GSE147672) (Supplemental Table S1; Corces et al. 2020; Morabito et al. 2021; Zhu et al. 2023). Count matrices in form of .h5 files were generated for each sample using Cell Ranger ATAC 1.2.0 imported into Seurat (Butler et al. 2018) and Signac (Stuart et al. 2021). Peaks in each sample were merged, and the counts were regenerated for the combined data set. Gene annotation for Homo sapiens hg38 was added with annotation package EnsDb.Hsapiens.v86 (https://bioconductor.org/packages/EnsDb.Hsapiens.v86/). Quality control was performed by removing cells with counts lower than 200 and higher than 35,000, nucleosome signal more than five, TSS enrichment less than two, blacklist ratio higher than 0.05, and percentage reads in peaks less than 15. We retained 224,131 and 31,616 cells for the retina and brain, respectively. Integration of all data sets within each tissue was performed by normalization and dimension reduction for each data set. First, features, peaks in this case, were ranked based on their observed frequency observed using FindTopFeatures (min.cutoff = 10). Normalization through term frequency-inverse document frequency (TF-IDF) was done to correct for sequencing depth, and dimension reduction was performed by running singular value decomposition (SVD) on the corrected count matrix. Anchors were identified using FindIntegrationAnchors (reduction = “rlsi,” dims = 2:30) (Stuart et al. 2019). The IntegrateEmbeddings function created a new cell embeddings layer in the Seurat object of each tissue for use with clustering and visualization.

Next, a UMAP layer was generated using the RunUMAP function (min.dist = 0.2, dims = 2:30, seed.use = 1) after performing unsupervised clustering and identification, as described above, with FindNeighbors and FindClusters (algorithm = 1, resolution = 0.5, random.seed = 1). Cells were annotated by generating gene activity matrices with the GeneActivity function, followed by the creation of an RNA assay. The count matrix in this assay was normalized using NormalizeData with default parameters. Finally, module scores based on our curated list of marker genes were computed and visualized using Seurat's DotPlot function. To maintain high specificity, cells lacking a definitive cell type label were removed from further analysis.

To combine the retina and brain Seurat objects, the peaks from both objects were first merged, and then the peaks were requantified in each tissue's respective object before combining them. A new gene activity matrix was calculated using the updated count matrix. The cell embeddings were recalculated following the steps outlined above, with the addition of RunHarmony (group.by.vars = ‘sample’) for batch correction (Korsunsky et al. 2019) prior to running RunUMAP.

Subclustering of retina and brain microglia from the snTAC-seq data

To further refine the annotation of microglia cells from the retina and brain, we extracted the microglia clusters from both tissues and applied the same procedures to generate UMAP plots of their subclusters. The Seurat object for the retina and brain initially labeled as microglia were split into lists of smaller objects for each data set. Each data set then undergoes cell embeddings generated using FindTopFeatures (min.cutoff = 10), RunTFIDF, and RunSVD. Anchors were created using FindIntegrationAnchors(reduction = “rlsi,” dims = 2:30) and additional parameters that are different: k.anchor = 100 for the brain and default for the retina. Finally, the data sets from the retina and brain separately were integrated with IntegrateEmbeddings (dims.to.integrate = 1:30) and k.weight = 50 for the retina and default for the brain. After integration, UMAP coordinates were generated by RunUMAP (reduction = “integrated_lsi,” dims = 2:30, seed.use = 1, min.dist = 0.2) (Stuart et al. 2019). Module scores were calculated based on differentially expressed genes and our curated list of marker genes, as described above, and were then visualized using the DotPlot function to annotate microglia, PVMs, and monocytes. After annotating the cells, the CallPeaks function was run to generate a set of peaks for each cell type in both tissues. These peaks were then filtered to include only those with conventional chromosome notations and to exclude any that overlapped with a blacklist provided in Signac (Stuart et al. 2021). To visualize the identified microglia, PVMs, and monocytes from the retina, brain, and combined, we extracted the three cell types from each Seurat object and then reintegrated them using the same steps as above, with some differences in parameters. For Retina, we changed k.weight = 40 in IntegrateEmbeddings. For the brain, k.anchor was changed to default, and k.weight was kept at default. For integrating both the retina and brain data sets, k.anchor was kept at default and k.weight = 20. RunUMAP was performed after these integrations with dims = 2:30, min.dist = 0.2, and seed.use = 05192025 (Stuart et al. 2019).

Pairwise differential accessibility region and overrepresented TF motif analyses

Pairwise differential accessibility region (DAR) analysis was performed using FindMarkers (logfc.threshold = 0.25, min.pct = 0.2, random.seed = 1, test = ‘wilcox’) to compare the cell type of interest against other cell types (Butler et al. 2018). DAR peaks were linked to genes using the ClosestRegion function from Signac (Stuart et al. 2021) to identify genes that overlapped with or were in close proximity to these peaks. Motif analysis on the retina and brain microglia data used the Signac package (Stuart et al. 2021). We acquired the DNA sequence motif information from JASPAR2022 and added into snATAC-seq seu.obj with AddMotifs function in Signac (Stuart et al. 2021). Finally, the FindMotifs function was used to identify overrepresented motifs within DAR peaks.

Hierarchical clustering of microglia subclusters

Hierarchical clustering for snRNA-seq and snATAC-seq data was performed using speciesTree (Bakken et al. 2021) to quantify the similarity of the microglia cluster in each tissue. Counts from either the gene expression matrix or chromatin accessibility peaks were used as input. The distance matrix was calculated using the cor function with method = “Spearman,” and the matrix was then used for hierarchical clustering in speciesTree. The resulting plot was generated using the default functions provided in the package.

Correlation heatmap for all cell types in the retina and brain

The average expression data from the RNA assay in the combined retina and brain Seurat object were generated using the AverageExpression function from Seurat. This matrix was then input into the cor function with method = “Spearman” to account for differences in gene expression scaling between the retina and brain. A heatmap was generated using ComplexHeatmap (Gu et al. 2016), and clustering was done using the hclust default method. The significance of the difference between correlations of immune cells with other cell types versus within immune cells was determined using the Mann–Whitney U test via wilcox.test function.

Ligand receptor identification

We predicted the ligand–receptor interaction between retina microglia and other retinal cell types, as well as brain microglia with other brain cell types, using LRLoop (Xin et al. 2022). Interaction strengths were first calculated based on the expression levels of ligands in sender cells and receptors in receiver cells. We then integrated cell type–specific signaling and regulatory networks by overlaying known interactions with gene expression data. The get_LRscores function (default parameters) was used to compute interaction scores between each pair of cell types in both directions. To refine these scores, we applied the updateLRscores function (default parameters), using microglia interaction scores and those from the partner cell type to account for regulatory loops. Finally, we filtered for interactions in which both ligand and receptor genes were expressed in at least 10% of their respective cell types and had a score of 0.5 or more across all time points.

Gene expression heatmap for all cell types in retina and brain using AMD– and AD–GWAS genes

Relevant gene lists for AMD and AD were obtained to create heatmaps of their expression in each cell type in the retina and brain. For AMD, the list comprised most the likely target genes at those loci identified through GWAS (Fritsche et al. 2016), eQTL (Orozco et al. 2020), and rare variant (Ratnapriya et al. 2020) studies. For AD, we gathered the reported nearest protein-coding genes to the associated GWAS variants (Bellenguez et al. 2022). Heatmaps were generated using ComplexHeatmap (Gu et al. 2016). Genes with scaled expression greater than 1.5 in microglia, relative to other cell types in the retina or brain, were further selected for comparison.

Pathway enrichment analysis

DEGs were filtered based on average fold change greater than two and an adjusted P-value 0.05. Venn diagrams were made using DeepVenn (Hulsen 2022). Pathway enrichment analysis was performed using a 5% FDR cutoff with Ensembl IDs as the key type and the biological process subontology (Yu et al. 2012). Top five GO terms with the lowest adjusted P-values were selected for plotting. Dot Plots were generated using dotplot function from enrichplot package.

MAGMA

The significance of microglia in the retina and brain in neurodegenerative diseases was evaluated. For snRNA-seq data, genes most specific to each cell type were selected (those showing at least a twofold higher expression compared with other cell types) and analyzed using MAGMA (de Leeuw et al. 2015). GWAS data were obtained from IAMDGC consortium for AMD (Fritsche et al. 2016) and the GWAS catalog for AD (Bellenguez et al. 2022). We computed the full gene analysis using all gene locations and the GWAS SNPs with default parameters. This step requires the LD reference, which we obtained from The 1000 Genomes Project Consortium (2015); the gene locations in hg19 provided in the MAGMA repository; and the summary statistics from the GWAS. We used the hg19 reference genome for our analysis, as the AMD–GWAS and AD–GWAS data sets were originally aligned to this build. Although newer reference genomes such as GRCh38 offer improved annotations, we do not expect them to significantly impact our conclusions, as the study focuses on common variants that are well represented across reference builds. Gene-set analysis was performed using the result from the full gene analysis along with the list of genes from our DE analysis, including genes that were present in >20% of the cells. For ATAC-seq data, genes overlapping or closest to the DAR peaks with a lower threshold of twofold change were selected for gene-set enrichment analysis using MAGMA (de Leeuw et al. 2015). We applied the Benjamini–Hochberg correction for multiple testing in R (R Core Team 2025) using p.adjust (method = “fdr”). The BETA value from the MAGMA output was used to represent the enrichment, and results were visualized in the form of a bar plot using log₁₀ of the adjusted P-value.

LDSC-SEG and S-LDSC

We used the same set of genes as in the MAGMA analysis, selecting cell type–specific genes with twofold or more expression. LD reference was acquired through the 1000 Genomes Project (The 1000 Genomes Project Consortium 2015) of European ancestry. For LDSC-SEG (Finucane et al. 2018), we added a 100 kb window to the hg19 gene body coordinates to create the binary annotation for LD score calculation for each set of DE genes. A separate LD score for all protein-coding genes was calculated. Partitioned heritability was estimated by stratified LD score regression using GWAS summary statistics for AMD and AD, the baseline model, LD score for all protein-coding genes, and individual LD score sets for each cell type in each tissue. For S-LDSC (Finucane et al. 2015), the annotation was generated by overlapping the peaks of interest with the SNPs’ position located in the LD reference to create a single column file, containing 1 sec and 0 sec for hits and misses, respectively. Peaks called for cells in microglia in the retina and brain are used for creating annotations. These peaks were first filtered for conventional chromosomes and through a blacklist provided by Signac, which were then lifted to hg19 coordinates. The baseline-LD score was provided by LDSC group (Finucane et al. 2015). After LD score computations, we ran stratified partitioned heritability calculations for AMD–GWAS or AD–GWAS summary statistics using the baseline-LD scores and each of the LD scores for those each annotation separately. P-values were adjusted for multiple comparisons in R (R Core Team 2025) using the Benjamini–Hochberg method via p.adjust (method = “fdr”). The results were visualized through bar plots of the log₁₀ of the adjusted P-value.

GARFIELD

We performed the functional enrichment analysis using GARFIELD (Iotchkova et al. 2019) of the open chromatin annotation in retina and brain by leveraging GWAS summary statistics of AMD and AD. We ran GARFIELD with default parameters and our annotations, and the method is described in detailed elsewhere (Iotchkova et al. 2019). Briefly, GARFIELD first performs the LD pruning to extract independent SNPs from genome-wide genetic associations, followed by LD tagging annotation to annotate each variant with regulatory annotations. Then, utilizing a logistic regression model, GARFIELD calculates enrichment statistics (odds ratios [OR]) and P-values at user-defined GWAS P-value thresholds for each annotation. Enrichment of significant SNPs was assessed for each tier of GWAS P-value thresholds: 1 × 10⁻⁵, 1 × 10⁻⁶, 1 × 10⁻⁷, and 1 × 10⁻⁸. The results are presented in the radial plots. For the ENCODE data, we downloaded narrowPeak files from 61 tissues (The ENCODE Project Consortium et al. 2020). Differences in the number of open chromatin regions across tissues could introduce potential bias in comparative heritability analyses. To minimize this effect, we initially focused on uniquely open regions, defined as those active in ≤10% of the tissues. However, for the AD–GWAS comparison of retina versus all other tissues, we increased this threshold to 20% for retina microglia, as the 10% threshold resulted in a very small number of peaks (fewer than 3000) owing to the smaller sample size.