Single nucleus multiomics identifies ZEB1 and MAFB as candidate regulators of Alzheimer's disease-specific cis-regulatory elements

Abstract

Cell type-specific transcriptional differences between brain tissues from donors with Alzheimer's disease (AD) and unaffected controls have been well documented, but few studies have rigorously interrogated the regulatory mechanisms responsible for these alterations. We performed single nucleus multiomics (snRNA-seq plus snATAC-seq) on 105,332 nuclei isolated from cortical tissues from 7 AD and 8 unaffected donors to identify candidate cis-regulatory elements (CREs) involved in AD-associated transcriptional changes. We detected 319,861 significant correlations, or links, between gene expression and cell type-specific transposase accessible regions enriched for active CREs. Among these, 40,831 were unique to AD tissues. Validation experiments confirmed the activity of many regions, including several candidate regulators of APP expression. We identified ZEB1 and MAFB as candidate transcription factors playing important roles in AD-specific gene regulation in neurons and microglia, respectively. Microglia links were globally enriched for heritability of AD risk and previously identified active regulatory regions.

Keywords: APP; Alzheimer’s; MAFB; ZEB1; cis-regulatory element; microglia; multiomics; neuron; single cell; transcription factor.

Graphical abstract

Figure 1

Cellular diversity of DLPFC from Alzheimer’s disease and unaffected donors revealed by single cell multiomics (A) Experimental design. (B) Uniform manifold approximation and projection (UMAP) visualization of the weighted nearest neighbor (WNN) clustering of single nuclei colored by cell type and cluster assignment. (C) Total number of cells in each subcluster and the proportion of cells from each individual (red, AD donors; blue, unaffected donors) in the subcluster. (D) Row-normalized gene expression of scREAD cell type markers. (E) Chromatin accessibility across cell types for cell type marker genes (indicated below). (F) Correlation of pseudo-bulked cell type-specific expression profiles between individuals. Colors indicating cell type are consistent throughout the figure.

Figure 2

Cell type-specific transcriptome dysregulation in Alzheimer’s DLPFC (A) MAST log₂(fold change [FC]) of all up- and downregulated genes in AD for each cell type. (B) Number of shared DEGs between cell types in both directions (upper triangle, upregulated in AD; lower triangle, downregulated in AD). (C) Normalized expression of the top DEG in the indicated cell types (log₂[FC] > 1). (D) Overlap of DEGs with agreement on cell type and direction with Morabito et al. and Mathys et al. (E) Heatmap showing the odds ratio of the top enrichR GO terms for up and downregulated DEGs within each cell type (∗adjusted p < 0.01).

Figure 3

Identification of candidate CREs (A) Schematic of gene-peak association (top). Heatmap of row-normalized accessibility and expression for the most correlated peak-gene link for each gene (bottom). Columns are pseudo-bulked on cell type and disease status. (B) Distributions of the number of linked peaks per gene (left) and the number of linked genes per peak (right) for AD (red) and control (blue) samples. (C) Total number of links per cell type for AD and control. Cell type of the link is assigned by the cell type in which the peak was called. (D) ENCODE annotation of linked peaks by cell type. (E) Shared (across cell types) and cell type-specific linked peaks that overlap H3K27ac of the corresponding cell type. (F) Normalized expression of KANSL1 from AD and control samples in each cell type. Expression is significantly different in AD versus control for all cell types. (G) Linkage plot for all links to KANSL1. Top: coverage plot of pseudo-bulked accessibility in excitatory neurons separated by status (red, AD; blue, control). Bottom: significant AD and control peak-gene links. Arc height represents strength and direction of correlation. Arc color indicates if the link was identified in both AD and control (common, gray) or control donors only (blue). A linked peak overlapping a single SNP is highlighted in gray.

Figure 4

Identification of AD-specific TF regulatory networks (A) Strategy for defining peak-gene-TF trios. A linked peak containing a TF motif must be correlated with that TF and the expression of that TF must be correlated with the linked gene for that peak to be considered a part of trio. (B) Genome annotations for location of linked peaks within trios. (C) Heatmap of column-normalized expression of genes within MEF2C trios by cell type. (D) Normalized expression of MEF2C by cell type. (E) Top enriched GO terms for genes within MEF2C trios from excitatory and inhibitory neurons (green, “neuron”) and microglia (purple, “microglia”). (F) Heatmap of correlation values of AD and control-specific trios identified in microglia (left) and excitatory/inhibitory neurons (right) for TFs involved in at least 3 trios. (G) Linkage plot for GABRA5. Top: coverage plot of pseudo-bulked accessibility in indicated cell types. Middle: coverage plot of ZEB1 ChIP-seq signal from NeuN⁺ nuclei isolated from DLPFC tissue from two unaffected donors (1238 and 1242). Bottom: significant peak-gene links; green indicates overlap with ZEB1 motif. Arc height represents strength and direction of correlation. Track of ZEB1 motifs (green) and H3K27ac peaks from neurons (black; Nott et al.¹¹). Linked peak of interest is highlighted in gray. (H) ZEB1 motif from JASPAR 2022 (top). Normalized expression of ZEB1 and GABRA5 in excitatory/inhibitory neurons and microglia.

Figure 5

Validation of candidate CREs (A) sLDSC results using 16 GWAS traits as indicated with our linked peaks stratified by cell type and group (“All” = all links, “Common” = links identified in both AD and control data, “AD” = links specific to AD, “Control” = links specific to control). Heatmap indicates coefficient Z score from running sLDSC with each set of links combined with the 97 baseline features. Feature-trait combinations with a Z score significantly larger than 0 (one-sided Z test with alpha = 0.05, p values corrected within each trait using Benjamini-Hochberg method) are indicated with a numeric value reporting the enrichment score. (B) Bar plot showing enrichment (±95% confidence interval [CI]) of links for previously nominated regulatory regions: active MPRA elements (blue), eQTLs where target gene is same as linked gene (pink), and HiC loops linking region to same target gene (green). MPRA, massively parallel reporter assay; NPC, neural precursor cells; ESC, embryonic stem cells. (C) Boxplots showing statistically significant (∗p < 0.05, ANOVA with Fisher’s LSD) elements representing links tested in luciferase assays. Luciferase elements are denoted by the linked gene for the nominated region. (D) Boxplots showing comparison of rs12445022 with its corresponding reference element linked to JPH3 (∗p < 0.05, ANOVA with Fisher’s LSD). (E) Top: normalized expression of APP in each cell type. Middle: coverage plot of accessibility in indicated cell types. Bottom: significant control (blue) and common (gray) peak-gene links to APP tested in luciferase assays. Arc height represents strength and direction of correlation. Links that contained CREs that increased expression of the luciferase reporter are highlighted in gray. (F) Boxplots showing all tested luciferase elements representing APP-peak links. Elements highlighted in gray are located within the APP gene body (∗p < 0.05, ANOVA with Fisher’s LSD).