Single-cell multiregion dissection of Alzheimer’s disease

1Alzheimer’s disease (AD) is characterized by pathological protein aggregation in a stereotyped pattern across multiple brain regions,. Post-mortem diagnosis of AD is staged by the severity and distribution of these pathological hallmarks: extracellular amyloid-β deposits and intracellular neurofibrillary tangles (NFTs) in neurons. Tangles are first seen in the entorhinal cortex (EC) (Braak stages I–II), then the hippocampus and thalamus (Braak stages III–IV) and finally the neocortex (Braak stages V–VI), a sequence that is typically synchronous with cognitive decline from mild cognitive impairment to severe dementia,,–. Understanding the cellular architecture of affected brain regions has important implications for early and region-specific therapeutic interventions and may shed light on the molecular mechanisms underlying the regional progression of pathology. Although some brain regions relevant to AD have been studied individually at scale or jointly in samples from a few individuals–, a comprehensive molecular characterization of region-specific differences in AD is currently lacking and could capture differences in regional molecular architecture– and region-specific neuronal and glial subtype alterations in AD and in cognitive resilience to AD pathology,–.

2Here we present a transcriptomic atlas of the human brain spanning six distinct anatomical regions from persons with and without Alzheimer’s dementia as a basis for understanding disease-associated differences. We profile the transcriptomes of over 1.3 million nuclei from the EC, hippocampus (HC), anterior thalamus (TH), angular gyrus (AG), midtemporal cortex (MT) and prefrontal cortex (PFC) from 48 individuals, 26 of whom have a pathologic diagnosis of AD. We annotate region-specific neuronal and glial subtype diversity, present an online resource for navigating this atlas (http://compbio.mit.edu/ad_multiregion) and provide mechanistic insights into cellular vulnerability, response and resilience to AD.

3To characterize cellular diversity in the human brain, and the genes, pathways and cell types that underlie AD progression across brain regions, we performed single-nucleus RNA-sequencing (snRNA-seq) analysis of nuclei isolated from 283 post-mortem brain samples across six brain regions from 48 participants in the Religious Order Study (ROS) or the Rush Memory and Aging Project (MAP) (together, ROSMAP; Fig. 1a). We selected 48 participants on the basis of pathologic diagnosis of AD (stratified by NIA-Reagan score of 26 (with AD) and 22 (without AD; labelled non-AD)) and on the basis of clinical diagnosis of Alzheimer’s dementia (n = 16) versus non-dementia (n = 32), (Fig. 1a, Extended Data Fig. 1a and Supplementary Table 1). From these 48 individuals, we profiled six brain regions: the EC (221,493 cells), which is affected in early AD (stages I–II); the HC (221,415) and TH (207,625), which are affected in mid-AD (stages III–IV); and the AG (220,409), MT (227,412) and PFC (254,721), which are affected in late AD (stages V–VI), for a total of 1.35 million transcriptomes of independent nuclei after removing doublets, low-quality cells and highly sample-specific clusters. We annotated 76 high-resolution cell types in 14 major cell type groups, including 32 excitatory neuron subtypes (436,014 nuclei, 32.2% of total) and 23 inhibitory subtypes (159,838 nuclei, 11.8% of total) (Extended Data Fig. 1b–d, Supplementary Figs. 1 and 2 and Supplementary Table 2). We characterized these cell types in terms of their transcriptome size and proliferative status, compared our atlas with previously published data across species– (Extended Data Fig. 1e,f and Supplementary Figs. 3–5) and identified broad cell type identity programs using non-negative matrix factorization (NMF) and transcriptional regulons using SCENIC, (Extended Data Figs. 2 and 3 and Supplementary Tables 3 and 4).

4To gain insights into the cellular architecture of the human brain, we investigated differences in the composition of major cell types between the six brain regions. The fraction of neurons increased significantly from the TH (14.4% neurons) to the three-layer allocortical HC (32.2%), the entorhinal periallocortex (36.6%) and the six-layered neocortical regions (AG, MT and PFC, 58.9%) (Fig. 1b–e and Supplementary Fig. 6). Glia, including astrocytes, oligodendrocytes, oligodendrocyte precursor cells (OPCs) and microglia/immune cells, tended to be less abundant in neocortical samples (Fig. 1b–e), in agreement with previous studies in humans, and mice, (Supplementary Fig. 7a–d). Differences in the composition of major cell types between regions were reproducibly observed across study participants, irrespective of the individual’s disease status (Supplementary Fig. 7e–h), suggesting that variability in the major cell type composition between regions is a fundamental characteristic of the human brain and is not affected by AD pathology.

5We first characterized the regional diversity of excitatory neuron subtypes, which were consistent across individuals and were either highly region-specific to the HC, EC and TH (7, 9 and 2 subtypes, respectively) or were predominantly shared across neocortical regions (12 subtypes) (Fig. 1f and Supplementary Fig. 8–12). Hippocampal subtypes included neurons from the highly structured CA1 and CA2/CA3 subfields and dentate gyrus and the more entorhinal-proximal subiculum and para/presubiculum areas. EC-specific subtypes that clustered separately from neocortical subtypes for the same layers were often marked by expression of RELN, TOX3 and GPC5, and contained subtypes from both the lateral (L2 RELN+GPC5+) and medial (L2 TOX3+POSTN+) EC– (Supplementary Fig. 10). Excitatory neurons in the TH were predominantly composed (74%) of a subtype (NXPH1+RNF220+) that was not observed in the neocortex and is predicted to be regulated by LHX9, SOX2, SHOX2 and TCF7L2, (Extended Data Fig. 2a,b and Supplementary Figs. 9e–i and 11n,o). We found that the thalamic–neocortex separation is conserved in mice and recapitulated both this divide and thalamic marker genes in independent single-cell, bulk and microarray data in both mice and humans,,, (Supplementary Fig. 12).

6In contrast to excitatory neuron subtypes, the majority of inhibitory neuron subtypes (22 out of 23 subtypes) were observed in all five cortical regions (Fig. 1g and Supplementary Figs. 13–17), although some inhibitory subtypes had regional biases, including PVALB+HTR4+ and CUX2+MSR1+ (enriched in neocortex), layer 6 SST+NPY+ (EC and HC) and GPC5+RIT2+ (EC), suggesting that there are significant differences in inhibitory neuron composition between the neocortex and allocortex (Fig. 1g and Supplementary Fig. 14). Moreover, in the HC, EC and MT, caudal ganglionic eminence-derived GABAergic neurons (VIP+LAMP5+) were significantly more abundant than medial ganglionic eminence-derived neurons (SST+PVALB+), but these two major clades were not significantly different in the PFC (Extended Data Fig. 1g). By contrast, the TH contained a single unique, thalamus-specific inhibitory subtype (MEIS2+FOXP2+) marked by genes that are involved in neurite outgrowth, such as the semaphorins SEMA3C and SEMA3E, DISC1 and SPON1, and receptors for serotonin (HTR2A), acetylcholine (CHRM2, CHRNA3) and glutamate (GRM3) (Fig. 1g,h and Supplementary Figs. 14 and 15). These genes were in a single inhibitory program (Inh-22, from NMF) that included the SCENIC-predicted subtype regulators FOXP2 and LEF1, (Extended Data Figs. 2c,d and 3b). We recapitulated this thalamic difference and program genes in the mouse thalamus and human lateral geniculate nucleus (dLGN) using previously published single-cell data (Supplementary Figs. 16 and 17). To validate the localization and specificity of markers of the thalamic inhibitory neuron subtype, we performed in situ hybridization for both FOXP2 and MEIS2 with GAD2 on TH and PFC post-mortem brain samples from four individuals, and found significant thalamus-specific co-localization of both marker genes with GAD2 (Fig. 1i).

7As thalamic MEIS2 neurons expressed several typically glutamatergic neuron genes, we determined glutamatergic and GABAergic module scores for every neuronal cell to further examine the chimeric nature of this subtype (Supplementary Fig. 15g–k and Supplementary Table 5). These scores matched the cortical excitatory and inhibitory split and were negatively correlated both across and within broad neuronal classes (Fig. 1j and Supplementary Fig. 15g,h). Both thalamic MEIS2+ inhibitory and NXPH1+ excitatory neurons had intermediate scores, placing them between the cortical excitatory and inhibitory clusters, suggesting that they are less polarized with regard to the expression of cortical glutamatergic versus GABAergic programs (Fig. 1j and Supplementary Fig. 15i–k). Predicted cell–cell communication interactions were mostly shared across multiple regions, but the thalamus had multiple differential interactions (Supplementary Figs. 18 and 19). The top thalamus-specific interactions were between excitatory NXPH1 and neuronal NRXN1 or NRXN3, whereas inhibitory neurons expressed NXPH1 in the other regions, suggesting that neurexophilin signalling swaps from excitatory neurons in the thalamus to inhibitory neurons in cortical brain regions (Extended Data Fig. 4).

8We next tested whether glial cells also had transcriptional differences between brain regions. We identified multiple transcriptionally distinct subsets for each major glial cell type and determined their characteristic marker genes (Fig. 2a and Supplementary Fig. 20–25). Among glial cell types, astrocytes had the highest regional heterogeneity, containing both highly neocortex-enriched (GRM3+DPP10+) and thalamus-enriched (LUZP2+DCLK1+) subtypes (Fig. 2a–c and Supplementary Fig. 20). Region-specific astrocyte subtypes were experimentally validated using RNA in situ hybridization (Fig. 2d) and confirmed by analysing a separate snRNA-seq dataset (Supplementary Fig. 23m). Cortical astrocytes were enriched for markers involved in glutamate processing and transport, whereas hippocampus- and anterior thalamus-enriched DCLK1 astrocytes had lower glutamate transporter activity and were enriched instead for focal-adhesion-related genes (Fig. 2c and Supplementary Fig. 25a). Thalamic astrocytes (LUZP2+) expressed GABA-uptake genes SLC6A1 and SLC6A11 at much higher levels compared with other subtypes, even though the proportion of inhibitory neurons was not markedly higher in the thalamus (Fig. 2c). Notably, the thalamic MEIS2+FOXP2+ interneurons shared multiple markers with neocortical GRM3 astrocytes, including GRM3, MEIS2 and VAV3,, (Supplementary Fig. 23n), suggesting that astrocytes in evolutionarily newer regions may share some functions with inhibitory neurons in older regions.

9We developed a method, single-cell decorrelated module networks (scdemon), to identify gene expression modules from highly correlated sets of genes in atlas-scale snRNA-seq datasets (Fig. 2e). Highly imbalanced cell type composition in single-cell datasets, in which rare cellular states are outnumbered by common cell types, can lead to under-recovery of gene–gene interactions, especially for genes that are expressed at low levels. To account for these issues, our method estimates a sample-decorrelated gene–gene correlation matrix, thresholds gene–gene pairs based on their sparsity and uses the adjusted matrix to identify modules of highly correlated genes (Methods). We used our method to identify modules both across all cells in the atlas and for each major cell type independently, and recovered modules expressed to varying degrees, ranging from identity modules for each glial cell type to a cell cycle module found in just 0.7% of microglia (Fig. 2f, Extended Data Fig. 5, Supplementary Figs. 26–36 and Supplementary Table 6). Cells expressing these modules were enriched for diverse aspects of our dataset, including cellular subtype identity (205 modules), brain region (156, with 77 thalamus specific and 34 EC specific), AD status (73), APOE genotype (78) and sex (24) (Fig. 2f and Supplementary Table 7). We hierarchically clustered modules across the cell types and found that many cell types expressed gene programs for cholesterol biosynthesis (C10), chaperones (C5), ribosomes (C1 and C2), ER protein processing (C7), oxidative phosphorylation (C18), synapse interaction (C16), and glycolysis and response to hypoxia (C20) (Extended Data Fig. 6a,b).

10Using this approach, we identified 32 modules in astrocytes, including an astrocyte-wide program (M9, expressed in >99% of astrocytes) marked by GPM6A and GPC5 and enriched for cell junction assembly, and subtype- and region-specific identity programs such as thalamus-associated M19 (SLC6A11, LGR6, MRAS), which were enriched for sonic hedgehog signalling, M12 (GRM3: forebrain neuron development) and M7 (DCLK1: synaptic membrane) (Fig. 2g,h and Supplementary Fig. 31). Other modules spanned a diverse set of functions, including metallostasis, RNA splicing, glycolysis, oxidative phosphorylation and cholesterol biosynthesis and were shared by multiple subtypes (Fig. 2i–k and Extended Data Figs. 5b and 6a–c). For example, chaperone-enriched and APOE-ε4-associated M8 was expressed in multiple different astrocyte subtypes and regions, and expression of AD-associated M28 (metallostasis) overlapped with expression of both APOE+ (M0) and reactive (M3) astrocytes (Fig. 2i–k). Module–module correlations across samples revealed co-expressed programs, such as reactive astrocytes (M3, marked by TPST1, CLIC4 and EMP1) with cholesterol biosynthesis module M17 (r = 0.60), and glycolysis (M6) with AP-1 module M13 (r = 0.39, including FOS/JUN and ubiquitin), a pair that is potentially co-expressed in astrocytes under metabolic stress (Extended Data Fig. 6d,e).

11In contrast to astrocytes, immune cells showed little regional specificity and oligodendrocyte-lineage cells had thalamus-enriched subtypes with minor transcriptomic differences to neocortex-enriched subtypes (Supplementary Figs. 21–25). Immune modules included identity programs, such as for T cells (M6), macrophages (M7) and cycling microglia (MKI67+, M5) as well as modules found across immune cells and enriched for genes involved in NF-κB signalling (M18), interferon (M20), p53 and DNA damage response (M22) and TGFβ signalling (M14) (Extended Data Fig. 6f and Supplementary Fig. 32). Oligodendrocyte-lineage modules showed high regional specificity, and two OPC modules—thalamus-enriched M11 and EC-enriched M25—were marked by synapse-associated genes such as neural adhesion-related SEMA3D, SEMA6D and CNTN5, and glutamate receptor GRIA4, suggesting a role for OPC sensation and response to neuronal activity in specific brain regions (Extended Data Fig. 5c,e and Supplementary Figs. 33, 34 and 36).

12After constructing our atlas across AD-affected brain regions, we examined how AD affects the cellular composition. At the level of major cell types, we observed slight, non-significant decreases in the number of both excitatory neurons (odds ratio (OR) = 0.94, individuals stratified by pathologic diagnosis of AD), inhibitory neurons (OR = 0.93) and OPCs (OR = 0.85), as well as an increase in the number of oligodendrocytes (OR = 1.14, adjusted P (Padj) = 0.01) and vascular cells (OR = 1.24), mostly driven by differences in the EC, HC and PFC regions, especially in late AD (Extended Data Fig. 7a,b). We next tested whether the fractions of region-specific neuronal subtypes were significantly altered relative to both individual-level pathologic and clinical diagnoses of AD and region-level NFT and plaque accumulation (Fig. 3a and Extended Data Fig. 7c,d). Among excitatory neurons, we identified one HC-specific (CA1 pyramidal neurons) and four EC-specific (L2 RELN+ lateral EC, L3 RELN+, L5 and L2/3 TOX3+TTC6+ neurons) subtypes that were significantly less abundant (OR = 0.38–0.66) in individuals with a pathologic diagnosis of AD (Fig. 3a and Supplementary Fig. 37a–c). Neocortical L2–3 neurons were also significantly less abundant in samples with high NFT levels and in individuals with neocortical NFT involvement (Fig. 3a). Individuals with lower percentages of these vulnerable excitatory neuron subtypes performed significantly worse on cognitive tasks, with the strongest observed impacts on episodic memory and global cognitive function for subtypes marked by GPC5 (EC L5 and L2 RELN+) (Extended Data Fig. 7e,f). Notably, while the overall excitatory fraction was not associated with cognition, lower OPC fraction across regions and, in particular, in non-neocortex regions was significantly associated with impaired cognition (Supplementary Fig. 37d).

13Given that these neuronal subtypes lie in highly interconnected regions, we next examined whether neuronal subtypes connected across regions were coordinately depleted. We found that vulnerable neuronal subtypes were co-depleted specifically in individuals with AD, with some of the strongest effects observed in established connections between the CA1, subiculum, EC–L3 and EC–L5 (Fig. 3b,c and Extended Data Fig. 7g). These included co-depletion for entorhinal L5 versus L5-projecting subiculum (Kendall’s τ = 0.37 (AD); −0.1 (non-AD)) or CA1 (τ = 0.42 (AD) and −0.16 (non-AD)); and for CA1 versus L2-lateral EC (LEC, τ = 0.26 (AD) and −0.07 (non-AD)) and L3 RELN+ (τ = 0.24 (AD) and −0.13 (non-AD)) EC neurons, both of which project in part to the CA1 subfield, (Fig. 3b,c and Extended Data Fig. 7g).

14We next investigate whether vulnerable subtypes share marker genes that might mediate their vulnerability, and identified 391 genes with significantly higher baseline (non-AD) expression in vulnerable subtypes (Fig. 3d and Supplementary Table 8). These included Reelin signalling pathway genes RELN and DAB1; kinase-associated genes MAP2K5, PRKCA and SPHKAP; and multiple genes associated with heparan sulfate proteoglycan biosynthesis (including HS6ST3, XYLT1 and NDST3) (Fig. 3d and Extended Data Fig. 7h,i). Notably, while RELN expression, which is typically restricted to inhibitory neurons, was highly specific to two EC excitatory subtypes, its downstream partner DAB1 was present across subtypes (Extended Data Fig. 7h,i and Supplementary Fig. 37e,f).

15We next examined whether vulnerable inhibitory neuron subtypes in the PFC share characteristics with vulnerable excitatory neuron subtypes across our brain regions using single-cell transcriptomes from 621 ROSMAP study participants,. We identified specific inhibitory neuron subtypes that are depleted in individuals with a high tangle density burden, consistent with our previous findings (Extended Data Fig. 7j). Vulnerable and non-vulnerable inhibitory neuron subtypes differed in the expression of genes involved in neuron projection morphogenesis (ROBO2, SEMA6A and EPHB6), enzyme-linked receptor protein signalling pathways (FGFR2, TGFBR1 and PLCE1) and heparan sulfate proteoglycan biosynthesis (Extended Data Fig. 7k and Supplementary Table 8). Notably, vulnerable inhibitory neuron subtypes expressed significantly higher levels of the Reelin signalling pathway components RELN and DAB1, mirroring the observed higher expression of these two genes in vulnerable excitatory neuron subtypes (Fig. 3e). Furthermore, the Reelin receptors LRP8 (also known as ApoER2) and NRP1 exhibited significantly different baseline expression in vulnerable compared with non-vulnerable inhibitory neuron subtypes (Fig. 3f).

16To test the selective vulnerability of Reelin-expressing excitatory neurons in AD, we performed in situ hybridization (RNAscope) analysis of Reelin and vGlut (excitatory neuron marker) in EC tissue samples from both patients with AD and healthy individuals without AD. We found a significant decrease in the percentage of Reelin-expressing excitatory neurons in the EC of individuals with AD (Fig. 3g). To determine whether this finding was conserved in animal models of AD, we used immunohistochemistry to assess the expression of Reelin in the EC of both 12-month-old App knock-in (KI) mice and 9-month-old Tau(P301S) transgenic mice. We found that, relative to wild-type littermate controls, App-KI mice and Tau(P301S) mice had a significantly decreased percentage of Reelin-positive neurons in the EC (Fig. 3h,i), in agreement with our human transcriptomic data suggesting a selective vulnerability of Reelin-expressing neurons (Fig. 3d–f).

17To understand how vulnerable subtypes are altered in AD, we computed differentially expressed genes (DEGs) for each excitatory neuron subtype (Methods and Supplementary Fig. 38a–c). We partitioned DEGs into sets associated with either vulnerable or non-vulnerable subtypes according to their expression levels in individuals with late AD (Extended Data Fig. 7l,m). DEGs linked to non-vulnerable subtypes were enriched for a diverse set of functions, including ubiquitin-ligase binding, heat-shock-family chaperones, ER protein processing and mediators of neuronal death, whereas vulnerability-associated DEGs were highly enriched only for mitochondrial oxidative phosphorylation but included CRK and NEUROD2, which are both associated with Reelin signalling, (Extended Data Fig. 7l,m and Supplementary Fig. 38d–f). Some DEGs associated with non-vulnerable subtypes had higher differential effect sizes in the vulnerable subtypes, and showed additional enrichment for aerobic glycolysis (including PGK1, LDHB and SLC2A3) and clathrin-mediated endocytosis (including AP2M1/AP2S1, OCRL and COPS8) (Extended Data Fig. 7m).

18To identify regional differences in cellular expression and function specific to individuals with pathologic AD, we computed DEGs for each major cell type in every region alone and across regions using a negative binomial linear mixed model framework, adjusting for both known covariates and potential unknown batch effects (Methods) (Extended Data Fig. 8a and Supplementary Table 9). Astrocytes and inhibitory and excitatory neurons showed the highest number of DEGs over all of the regions, with the largest number of changes in the EC (Extended Data Fig. 8a). Notably, neuronal DEGs showed little overlap across regions, indicating that neuronal differences in AD are primarily determined by subtype or region of origin (Extended Data Fig. 8b). By contrast, microglia and OPC DEGs overlapped within the non-neocortex regions, and astrocyte and oligodendrocyte changes were more consistent across all regions (Extended Data Fig. 8b). AD DEGs were consistent with published results both for region-specific DEGs and for DEGs computed jointly over all regions for multiple AD variables, and were further corroborated by comparisons with various independent studies,,,,– (Supplementary Fig. 39).

19Excitatory DEGs were strongly enriched for electron-transport functional terms across regions and showed weak region-specific enrichments for protein-folding-, ubiquitination- and synapse-associated terms (Extended Data Fig. 8c). Inhibitory DEGs were also broadly enriched for protein-folding- and synapse-associated terms and for oxidative phosphorylation uniquely in the thalamus (Extended Data Fig. 8c). While microglia DEGs were broadly enriched for clathrin-coated endocytosis (up) and viral response (down), they also had diverse region-specific enrichments, including upregulation of major histocompatibility complex type II (MHC-II) binding in the EC and HC, RNA processing in thalamus and glycolysis in the PFC and EC; and HC-driven downregulation of phagocytosis, phospholipase signalling and protein kinase activity (Extended Data Fig. 8c).

20The majority of region-specific DEGs was either broadly shared (on average, 11% of genes were differentially expressed in 3+ cell types in a region) or were in cell-type-specific programs (40% of DEGs were in 3+ regions for a cell type) (Extended Data Fig. 8d,e). Such genes included SLC38A2 and EIF4G2 (broadly shared across regions) and PRDX5, HLA-DRA or CD44, upregulated DEGs in excitatory neurons, microglia and astrocytes, respectively (Extended Data Fig. 8f and Supplementary Fig. 40). Broadly shared genes across cell types showed region-specific enrichment, including for DNA damage (EC), amyloid-β binding and iron transport (HC) and glycolysis (thalamus) in upregulated genes as well as for phospholipid biosynthesis and autophagy in downregulated genes (Extended Data Fig. 8e). Gene sets based on DEGs for global AD pathology burden in the PFC across 427 individuals changed consistently in each region and glial cell type across global pathology, indicating that a significant component of the glial AD response is consistent across regions (Extended Data Figs. 8b–e and 9a,b). The remaining regional DEGs (on average, 48% of DEGs) highlighted region- and cell-type-specific changes. In microglia, these included PPARG and MSR1, upregulated in the HC, each associated with microglia polarization, as well as upregulation of lipoprotein modifier APOC1 and downregulation of transcription factor FOXP2 in the EC (Extended Data Fig. 8f and Supplementary Fig. 40).

21We next examined which cell types and regions were most enriched for genes identified in genome-wide association studies (GWASs) of AD by computing GWAS scores for each cell using single-cell disease-relevance score (scDRS),. Microglia and immune cells showed consistently high scores across regions, with the top scores for the microglia TPT1+ subtype and macrophages in the HC, thalamus and AG (Extended Data Fig. 9c). We examined whether GWAS genes showed region-specific differences in expression that might be linked to the region specificity of AD progression. We identified eight GWAS genes with region-specific expression in microglia, including PLCG2 (EC), APOE and SORL1 (thalamus), and MS4A4A (midtemporal cortex) (Extended Data Fig. 9c–f).

22To determine whether GWAS-identified genes have regional associations with Alzheimer’s pathology, we intersected DEGs for regional pathology measurements with 149 identified familial AD and GWAS locus genes– (Extended Data Figs. 10 and 11a). We found that 74 genes (49%) were differentially expressed for at least one cell type, and multiple genes showed region-specific expression, including the lipid transporter ABCA7 (enriched in thalamus), the zinc-finger protein ZNF655 (EC) and the complement receptor CR1 (neocortex) (Extended Data Fig. 10). GWAS and familial AD genes were maximally expressed (75 genes) and differentially expressed (30 genes) in microglia, and 25 genes were differentially expressed in at least three cell types, including upregulated CLU, PLCG2 and SORT1, and downregulated DENND6A (Extended Data Fig. 10). Among all of the cell types, astrocytes and microglia showed the largest differential changes for these genes in regions with high neuritic plaque density, for example, for APOE, HLA-DRA, PILRA and SORT1, and showed the most response to diffuse plaque. Neurons and oligodendrocyte-lineage cells showed stronger differences for these genes, including for PLCG2, CLU and MAF, in regions with high NFT density (Extended Data Fig. 10).

23To determine whether different pathologies induce distinct transcriptional responses, we computed DEGs for region-specific measurements of NFT and neuritic amyloid-β plaque burden (measured in each region except thalamus) (Fig. 4a, Extended Data Fig. 11a,b and Supplementary Fig. 41). DEGs for AD pathology showed a high overlap with DEGs for pathologic diagnosis (NFT: 45% and plaque: 53% on average) (Fig. 4a). Agreement between NFT and plaque DEGs was highest in the EC and HC for all cell types (average adjusted R2 of 67% in both) and lowest in the PFC (43%) and AG (21%), consistent with late-AD NFT appearance in the neocortex (Fig. 4b).

24We next identified genes with higher differential effects in either NFTs or neuritic plaques (Fig. 4c, Extended Data Fig. 11c–f and Supplementary Table 10). Consistently, NFT-associated genes (374 genes, differentially expressed in 2+ cell types) included PLCG2, CLU and CTNNA2 (in oligodendrocytes and OPCs) and mitochondrial subunits, and were enriched for ER protein processing, electron transport and cadherin binding (Fig. 4d,e). Neuritic-plaque-associated genes (190 genes) included the energy-homeostasis-regulating genes IRS2, PDK4 and HIF3A, and genes enriched for immune response, chromatin regulation and lipid droplets. Notably, in excitatory neurons, plaque-associated and upregulated DEGs were strongly enriched for aerobic transport chain components (including NDUFA4 and COX6B1) (Fig. 4f,g). On the other hand, NFT-associated and downregulated DEGs were enriched for TCA cycle genes, whereas upregulated DEGs were enriched for unfolded protein response and lysosome-linked genes. Finally, astrocytes contained more plaque-associated DEGs compared with other cell types, and their pathology-associated DEGs were enriched in our expression modules, including in metallostasis (M28) for plaque-associated DEGs and oxidative phosphorylation (M27) and chaperones (M8) for NFT-associated DEGs (Fig. 4c,h–j). Interestingly, a reactive astrocyte module (M3) was enriched in upregulated genes for plaques but in downregulated genes for NFTs (Fig. 4j).

25Given the enrichment of NFT-associated or plaque-associated DEGs in expression modules, we next examined whether gene modules were enriched for AD DEGs (for AD pathology or for AD diagnosis) (Fig. 4k, Extended Data Fig. 11g and Supplementary Fig. 42). The same modules enriched for pathology-associated astrocyte DEGs were also enriched for the full sets of DEGs, including metallostasis (M28) with neuritic plaque DEGs and oxidative phosphorylation (M27) with NFT DEGs (Fig. 4k). Modules including ECM, adhesion and neurogenesis-related genes were much lower in AD (M1 and M11), while the modules for specific astrocyte subtypes (M7, DCLK1+; and M24, DPP10+) were enriched for upregulated DEGs (Fig. 4k).

26We independently identified modules for heat-shock chaperones, glycolysis and oxidative phosphorylation in multiple cell types, which were correlated across cell types and were enriched for upregulated DEGs (Fig. 4k, Extended Data Figs. 6a,b and 11g and Supplementary Fig. 42a,b). The glycolysis modules were enriched among diffuse plaque DEGs in microglia and astrocytes and shared a set of genes that included canonical glycolysis genes (PDK1/3/4, PFKL/P, PKM and PGK1), anaerobic glycolysis enzymes (TPI1 and LDHA) and stress-induced genes (EGLN1, DDIT4, VEGFA and BNIP3L) (Fig. 4k–m and Extended Data Fig. 6a,b). All glial types upregulated core glycolysis driver GAPDH and mitophagy-regulating BNIP3L in response to NFT burden and in individuals with cognitive impairment (Fig. 4m). In regions with high diffuse plaque, astrocytes upregulated glycolysis enzymes converting glucose-6-phosphate to pyruvate, while downregulating MPC1, the mitochondrial pyruvate transporter (Fig. 4n). In parallel, astrocytes uniquely upregulated DDIT4, PFKP and ADCY8, along with genes that suppress fatty acid metabolism (ANGPTL4) and promote lipid droplet storage of fatty acids (HILPDA), while microglia upregulated multiple glycogen-related genes (GBE1, UGP2 and PYGL) (Fig. 4m).

27To validate differential expression of ADCY8 and PFKP, we performed in situ hybridization (RNAscope) in AG tissue samples from patients with AD and control individuals without AD and found a significant increase in transcripts of both genes in AQP4+ astrocytes (Fig. 4o,p). Finally, we noticed that the glycolysis pathway genes were maximally expressed at different points in global AD progression for each region (pathology diagnosis by ABC score),. The pathway peaked very early in the EC (ABC score of 1, low levels of AD pathology), later in the HC and midtemporal cortex (intermediate levels), and very late in the PFC (high levels) (Supplementary Fig. 42c), suggesting that the glial metabolic response to AD may not be coordinated globally.

28In addition to understanding cellular alterations associated with specific pathological measures in AD, we investigated what transcriptional changes are associated with cognitive resilience (CR) in AD, cases in which individuals with AD brain pathology display much less cognitive impairment than expected,–. To identify potential molecular mediators that confer CR to AD pathology, we defined CR either categorically as the absence of cognitive impairment despite a pathologic diagnosis of AD (clinical diagnosis condition), or continuously, as the difference between observed cognition and the cognition expected on the basis of pathology level (Fig. 5a). We computed both scores for CR based on global cognitive function and for cognitive decline resilience (CDR) based on the rate of change of global cognitive function over time, and used four different measures of AD pathology: global AD pathology, neuritic plaque burden, NFT burden and tangle density (Fig. 5a).

29We calculated DEGs for both CR and CDR in each major cell type in the PFC (snRNA-seq from 427 ROSMAP study participants). Astrocytes were the only cell type with a consistently high number of genes associated with CR across all of the measures tested (Fig. 5b). To identify specific molecular pathways within astrocytes that may contribute to CR, we focused on genes that are consistently associated with multiple measures of CR in astrocytes (termed CR-associated genes). Several CR-associated genes, including GPX3, HMGN2, NQO1 and ODC1 (encoding a rate-limiting enzyme of polyamine biosynthesis), possess or promote antioxidant activities– (Fig. 5c–f and Supplementary Fig. 43a–d). The expression of HMGN2, NQO1, ODC1 and GPX3 in astrocytes was also positively associated with cognitive function (Fig. 5g and Supplementary Fig. 43e), and these genes exhibited the highest expression in astrocytes isolated from those individuals with the least cognitive decline over time (Fig. 5h and Supplementary Fig. 43f). Analysis of bulk RNA-seq data from the ROSMAP cohorts (n = 638) confirmed a significant positive association between the expression level of HMGN2, ODC1 and GPX3 and multiple measures of cognitive function and CR to AD pathology (Supplementary Fig. 44a–d).

30Furthermore, we noticed that several CR-associated genes within astrocytes encode enzymes that catalyse metabolic reactions that are involved in choline formation and breakdown. The expression of GPCPD1, which encodes glycerophosphocholine phosphodiesterase 1, an enzyme that is critical for cleaving glycerophosphocholine (GPC) to produce choline, was positively associated with measures of CR in astrocytes (Fig. 5c–f and Supplementary Fig. 43a–d). Conversely, PNPLA6, which encodes a phospholipase that catalyses the hydrolysis of intracellular phosphatidylcholine, a major membrane lipid, generating GPC, and CHDH, which encodes choline dehydrogenase, an enzyme that catalyses the conversion of choline to betaine aldehyde, were both negatively associated with multiple measures of CR in astrocytes (Fig. 5c–f and Supplementary Fig. 43a–d). Many of the CR-associated genes identified in PFC astrocytes were also associated with CR in astrocytes from other regions of the human brain (Fig. 5i and Supplementary Fig. 45), corroborating a link between astrocytes and CR beyond the PFC.

31To validate the choline pathway genes PNPLA6, GPCPD1 and CHDH, we selected PFC samples from individuals with high amyloid and tau pathology and compared transcript levels between individuals with intact cognition (that is, cognitively resilient) to those with cognitive impairment, and performed in situ hybridization (RNAscope) analysis of these genes with AQP4 as a marker for astrocytes. We found a significant decrease in PNPLA6 and CHDH transcripts and a significant increase in GPCPD1 transcripts in cognitively resilient individuals, in agreement with the differential expression results (Fig. 5j–l). Notably, choline oxidation to betaine generates a labile methyl group that can be used for homocysteine remethylation, resulting in methionine formation, which is subsequently transformed into the universal methyl donor S-adenosylmethionine. S-adenosylmethionine is involved in the biosynthesis of spermidine, linking choline metabolism and polyamine biosynthesis in astrocytes in CR to AD pathology (Fig. 5m).

32Here we present a transcriptomic atlas of the aged human brain—spanning six brain regions from 48 individuals with and without a diagnosis of AD—that we used to annotate regional cellular diversity, identify gene expression programs and differences in AD across cell types, and pinpoint region-specific cell populations that are vulnerable to AD. We provide an interactive website for exploring the atlas and these annotations, markers, functional modules and differences in AD at both the single-cell and pseudo-bulk levels (http://compbio.mit.edu/ad_multiregion).

33By annotating neuronal and glial subtypes by brain region, we found significant compositional differences between regions, including a subtype of thalamic GABAergic neurons (MEIS2+FOXP2+) that is molecularly distinct from the canonical subclasses of inhibitory neurons in the neocortex. We used region-specific measurements of AD pathology to identify changes in gene expression associated with neurofibrillary tau tangle or amyloid-β plaque burden, including plaque-associated upregulation of metallostasis in astrocytes and of the electron-transport chain in excitatory neurons. We found that AD-risk genes were highly perturbed in AD—in particular for microglia, consistent with their enrichment for GWAS signal—but few risk genes showed region-specific expression. To further examine the cellular and regional heterogeneity of the human brain, we developed a scalable method, scdemon, which uses sample decorrelation to annotate both ubiquitous and rare gene expression programs in each major cell type, and used annotated modules to identify functional programs associated with specific pathological variables, including a glycolysis- and energy-metabolism-linked program in glia, associated with diffuse plaque burden.

34We identified five excitatory neuron subtypes that were reduced in patients with AD (vulnerable subtypes) in the early-affected EC and HC,,, including EC layer II (L2), RORB-positive L5 (AGBL1+GPC5+) and hippocampal CA1 subfield neurons–. Notably, vulnerable excitatory neurons shared expression of genes involved in Reelin signalling and heparan sulfate proteoglycan biosynthesis, both of which were also predictive of inhibitory neuron vulnerability to AD. Recent case studies have identified variants in RELN and APOE as potentially mediating CR to autosomal-dominant AD. Notably, the RELN variant enhanced its binding to glycosaminoglycans (GAGs) and NRP1, and the APOE variant decreased binding to GAGs, potentially affecting their ability to compete for receptor binding,. Thus, our findings suggest a convergence of factors associated with cellular vulnerability in sporadic AD, and resilience to autosomal-dominant AD.

35Finally, we analysed the transcriptomic correlates of cognition and pathology in AD, and identified a set of astrocyte genes linked to CR to AD pathology. Notably, these genes converged on the pathways of choline metabolism and polyamine biosynthesis. This finding aligns with studies showing benefits of dietary choline intake and supplementation on cognitive performance in human individuals and in animal models–. Similarly, dietary supplementation with the polyamine spermidine prolongs life span and health span in several animal models, and spermidine has also been shown to enhance memory performance and counteract age-related cognitive decline–. Our findings support choline metabolism and polyamine synthesis as attractive targets for promoting CR in AD.

36Our study has several limitations: isotropic fractionation and read depth cut-offs may bias cell recovery based on their nuclear content; nuclear RNA may not fully capture microglial states or localized transcriptomic changes; and pathology burden is based on per-sample averages instead of on the spatial context of each cell. Additional individuals and data modalities will strengthen future analyses of region-specific alterations in AD, and spatial data may help to further separate pathology-associated changes.

37No statistical methods were used to predetermine sample size. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.

38We selected 48 individuals from ROSMAP, both ongoing longitudinal clinical–pathological cohort studies of ageing and dementia, in which all of the participants are brain donors. The studies include clinical data collected annually, detailed post-mortem pathological evaluations, and extensive genetic, epigenomic, transcriptomic, proteomic and metabolomic bulk-tissue profiling. For the purpose of this study, individuals were selected based on the modified NIA-Reagan diagnosis of AD and the Braak stage score (Braak stages 0, 1 and 2, n = 20; Braak stages 3 and 4, n = 14; Braak stages 5 and 6, n = 14), with 26 individuals having a positive pathologic diagnosis of AD and 22 individuals having a negative pathologic diagnosis of AD. Details of the clinical and pathological data collection methods have been previously reported,,,,. Individuals were balanced between sexes (male:female ratios 13:13 in AD, 11:11 in NoAD), matched for age (median, 86.6 years (AD) and 86.0 years (no AD)) and post-mortem interval (median, 5.9 h (AD) and 6.3 h (no AD)). Informed consent was obtained from each participant, and the Religious Orders Study and Rush Memory and Aging Project were each approved by an Institutional Review Board (IRB) of Rush University Medical Center. The participants also signed an Anatomic Gift Act, and a repository consent to allow their data to be repurposed.

39All dissections were done on a bed of dry ice using either a fine-toothed razor saw (for cortical regions) or a jewellers saw with diamond wire (for subcortical regions). Region-specific descriptions are as follows. (1) AG: full thickness cortex from the AG (Brodmann area: BA 39/40); take from the first or second slab posterior to the end of the HC. Minimize white matter. (2) MT: full thickness cortex from the middle temporal gyrus (BA 22); take as close to the level of the anterior commissure as possible. Minimize white matter. (3) PFC: full thickness cortex from the frontal pole (BA 10); take from the lateral side of the first or second slab. Minimize white matter. (4) EC: full thickness cortex from the EC (BA 28); take at the level of the amygdala. Avoid amygdala. Minimize white matter. (5) Posterior HC: take from the last slab containing HC. If the last slab has less than 5 mm of HC, take from the next slab anterior. Collect a full cross section. (6) TH: take from the first slab with thalamus. Include the most medial aspect.

40The protocol for the isolation of nuclei from frozen post-mortem brain tissue was adapted from a previous study. All of the procedures were performed on ice or at 4 °C. In brief, post-mortem brain tissue was homogenized in 700 µl homogenization buffer (320 mM sucrose, 5 mM CaCl2, 3 mM Mg(CH3COO)2, 10 mM Tris HCl pH 7.8, 0.1 mM EDTA pH 8.0, 0.1% IGEPAL CA-630, 1 mM β-mercaptoethanol and 0.4 U µl−1 recombinant RNase inhibitor (Clontech)) using a Wheaton Dounce tissue grinder (15 strokes with the loose pestle). The homogenized tissue was then filtered through a 40 µm cell strainer, mixed with an equal volume of working solution (50% OptiPrep density gradient medium (Sigma-Aldrich), 5 mM CaCl2, 3 mM Mg(CH3COO)2, 10 mM Tris HCl pH 7.8, 0.1 mM EDTA pH 8.0 and 1 mM β-mercaptoethanol) and loaded on top of an OptiPrep density gradient (750 µl 30% OptiPrep solution (30% OptiPrep density gradient medium, 134 mM sucrose, 5 mM CaCl2, 3 mM Mg(CH3COO)2, 10 mM Tris HCl pH 7.8, 0.1 mM EDTA pH 8.0, 1 mM β-mercaptoethanol, 0.04% IGEPAL CA-630 and 0.17 U µl−1 recombinant RNase inhibitor)) on top of 300 µl 40% OptiPrep solution (40% OptiPrep density gradient medium, 96 mM sucrose, 5 mM CaCl2, 3 mM Mg(CH3COO)2, 10 mM Tris HCl pH 7.8, 0.1 mM EDTA pH 8.0, 1 mM β-mercaptoethanol, 0.03% IGEPAL CA-630 and 0.12 U µl−1 recombinant RNase inhibitor). The nuclei were separated by centrifugation (5 min, 10,000g, 4 °C). A total of 100 µl of nuclei was collected from the 30%/40% interphase and washed with 1 ml of PBS containing 0.04% BSA. The nuclei were centrifuged at 300g for 3 min (4 °C) and washed with 1 ml of PBS containing 0.04% BSA. The nuclei were then centrifuged at 300g for 3 min (4 °C) and resuspended in 100 µl PBS containing 0.04% BSA. The nuclei were counted and diluted to a concentration of 1,000 nuclei per μl in PBS containing 0.04% BSA.

41For droplet-based snRNA-seq, libraries were prepared using the Chromium Single Cell 3′ Reagent Kits v3 according to the manufacturer’s protocol (10x Genomics). The generated snRNA-seq libraries were sequenced using NextSeq 500/550 High Output v2 kits (150 cycles) or NovaSeq 6000 S2 reagent kits.

42Gene counts were obtained by aligning reads to the GRCh38 genome using Cell Ranger software (v.3.0.2) (10x Genomics). To account for unspliced nuclear transcripts, reads mapping to pre-mRNA were counted. After quantification of pre-mRNA using the Cell Ranger count pipeline, the Cell Ranger aggr pipeline was used to aggregate all libraries (without equalizing the read depth between groups) to generate a gene–count matrix. The Cell Ranger v.3.0 default parameters were used to call cell barcodes. We used SCANPY to process and cluster the expression profiles and infer cell identities. We retained only protein-coding genes and filtered out cells with over 20% mitochondrial or 5% ribosomal RNA, leaving 1.47 million cells over 48 individuals and 283 samples across all regions. We further filtered the dataset to the top 5,000 most variable genes and used them to calculate the low dimensional embedding of the cells (UMAP) (default parameters, using 50 principal components and 15 nearest neighbours) and clusters using the Leiden clustering algorithm at a high resolution (15), giving 337 preliminary clusters. We separately called doublets using DoubletFinder and flagged and removed clusters with strong doublet profiles and clusters showing strong individual-specific batch effects, leaving a final dataset of 1.35 million cells.

43For the UMAP visualization of individual major cell type classes (excitatory neurons, inhibitory neurons, astrocytes, oligodendrocytes, OPCs, immune cells), the SCTransform-based integration workflow of Seurat was used to align data from individual samples, using the default settings,. We selected the set of relevant principal components on the basis of Elbow plots. We annotated cell types using previously published marker genes and single-cell RNA-seq data,,,–. To annotate cell types on the basis of previously published single-cell RNA-sequencing data (Allen Institute’s cell types database; https://portal.brain-map.org/atlases-and-data/rnaseq/human-multiple-cortical-areas-smart-seq), we used three separate approaches. First, Spearman rank correlation coefficients between the average expression profiles of neuronal subpopulations previously defined by the Allen Brain Institute and the neuronal subtypes identified in this study were computed using the cor function in R. Second, to project annotations of neuronal subpopulations previously defined by the Allen Brain Institute onto the neuronal cells analysed in this study, we followed the integration and label transfer workflow of Seurat. Third, we determined cell type marker genes based on data published by the Allen Brain Institute using the FindAllMarkers function from Seurat (Wilcoxon rank-sum test with Bonferroni correction for multiple testing; Padj < 0.05) and computed module scores for each cell type marker gene set across all neuronal cells analysed in this study using the AddModuleScore function of Seurat. To further annotate cell types, we determined marker genes using the FindAllMarkers function from Seurat (Wilcoxon rank-sum test with Bonferroni correction for multiple testing; Padj < 0.05). We tested only genes that were detected in a minimum of 25% of the cells within the cell type (min.pct = 0.25) and that showed, on average, at least a 0.25-fold difference (log-scale) between the cells of the cell type and all remaining cells (logfc.threshold = 0.25). Marker genes of the high-resolution cell types or states were determined separately for each major cell type class. We additionally compared the EC excitatory neuron subtypes to cell type annotations previously reported previously, which were computed using ACTIONet, and compared microglial markers to previously reported subtypes,.

44G2/M and S phase cell cycle scores were determined using the function CellCycleScoring in Seurat. Histograms showing the distribution of the G2/M- and the S phase module scores in each major cell class were generated using Prism 9 software. The statistical analyses comparing the number of genes detected per cell and the number of unique transcripts (UMIs) detected per cell between cell types was performed using Prism 9 software.

45Single-cell transcriptomic data from the human dLGN were obtained from the Allen Brain Institute (https://portal.brain-map.org/atlases-and-data/rnaseq/comparative-lgn). Single-cell transcriptomic data from multiple cortical areas and the hippocampal formation of the mouse brain were obtained from the Allen Brain Institute (https://portal.brain-map.org/atlases-and-data/rnaseq/mouse-whole-cortex-and-hippocampus-10x). Single-cell transcriptomic data across nine regions in the adult mouse brain were obtained from the McCarroll and Macosko Labs (http://dropviz.org/). Single-cell transcriptomic data from the mouse nervous system were obtained from the Linnarsson laboratory (http://mousebrain.org/adolescent/downloads.html). The external datasets and the human multiregion data presented in this study were integrated using the reciprocal PCA (RPCA) pipeline with the default parameters in Seurat (https://satijalab.org/seurat/articles/integration_rpca.html). The integration of single-cell data was performed separately for astrocytes, excitatory neurons and inhibitory neurons. For the integration of GABAergic neurons, the single-cell transcriptomic data from multiple cortical areas and the hippocampal formation of the mouse brain were downsampled to 50,000 GABAergic neurons. For the integration of excitatory neurons, the human multiregion dataset was downsampled to 5,000 neurons per high-resolution cell type. The mouse cortical data were downsampled to 50,000 excitatory neurons. The frontal cortex, posterior cortex, HC and thalamus data of the DropViz dataset were combined and downsampled to 50,000 neurons. Downsampling of the data was performed using the Seurat function subset with the default parameters. The comparison of microarray data from different human brain regions was performed using the Differential Search tool of the Allen Brain Atlas data portal (https://human.brain-map.org/microarray/search). The thalamus was selected as the target structure and compared to the cerebral cortex as the contrast structure. The differential search results including the fold change values and P values of the top 2,000 genes were downloaded from the data portal.

46We would like to find gene expression modules by calculating gene–gene correlations in single-cell data and using these to detect communities of similarly expressed genes. However, in single-cell data, which often contain an unbalanced composition of cell types, modules computed using this approach will be dominated by the most common cell type markers and pathways. Moreover, correlation values will often be inflated for pairs of sparsely expressed genes. We developed a method which accounts for these pitfalls to call multiresolution gene expression programs from single-cell data using an SVD-based approximation of zero-phase component analysis (ZCA) and gene sparsity-dependent thresholding,.

47The preprocessing transformations alternately called decorrelation, whitening, or sphering, transform a matrix X with a matrix W such that the covariance of XW is the identity matrix. In particular, ZCA is the transformation which maximizes the similarity of the transformed data to the original, which is achieved by setting W = C−1/2, where C is the covariance of X.

48In single-cell data, given a count matrix X with n cells (rows) and g genes (columns), we would like to perform ZCA decorrelation on the samples as a preprocessing step for calling modules. Computing and storing the n by n sample-wise covariance Cn = XXT/g is prohibitively expensive for modern datasets (with n > 1 × 106), even without centring X. Instead, we can analytically approximate the covariance with the SVD of XT = UnSnVnT as Cn ≈ (UnSnVnT)T(UnSnVnT)/g = VnSn2VnT/g and therefore Cn−1/2 = g1/2VnSn−1VnT. The ZCA transformation can then be computed as XZCA = Cn−1/2X = g1/2VnSn−1VnTX before calculating the covariance of XZCA for downstream analysis. While this approximation is already tractable for small single-cell datasets, we may not be able to compute the matrix multiplications or centred SVD for larger datasets. Here, we can use the SVD of X = USVT, which is commonly calculated in single-cell analyses, to approximate Cn−1/2 as g1/2US−1UT and XZCA = g1/2US−1UTX. From this, the non-centred covariance of XZCA is CZCA = XZCATXZCA/n = g × (US−1UTX)T(US−1UTX)/n. By substituting the SVD in for X, this reduces to CZCA = g × VVT/n, which is a simple and efficient approximation for very large single-cell data. As this approximation commonly drops out the largest identity program in the data due to the decorrelation approach, we allow computing CZCA = g × VSpVT/n, for any p, to tune the relative involvement of the larger eigenvalue components of the SVD.

49To control for inflated correlation estimates in highly sparse genes, we bin the estimated correlations (CZCA) for each pair of genes according to their sparsities (fraction of cells expressing the gene, binned on a log10 scale). We calculate the mean and s.d. for each 2D bin and smooth the estimates by fitting bivariate splines to the binned statistics, weighted by the log number of examples in the bin. We use the smoothed estimates to z-score the correlation matrix (CZCA), which we then threshold with a single z-score cut-off to create an adjacency matrix for a gene–gene graph. Graphs are laid out using the Fruchterman–Reingold algorithm and we remove connected components with fewer than four genes. We then use the leidenalg package and the Leiden algorithm with an RBConfiguration vertex partition to cluster the graph into gene modules. To robustly estimate modules for each set of cells in our analyses, we first performed a grid search for the optimal number of SVD components for cells of that type. We then computed the z-scored matrices for each of 10 bootstraps, selecting 90% of batches for each bootstrap and using only genes expressed in over 5% of cells in the full dataset for the cell type. We thresholded the average of the bootstrapped z-score estimates with z = 4.5 to build a graph. To balance the contributions of modules across the compositional spectrum, we calculated and thresholded separate graphs for eigenvalue powers p = 0, 0.25, 0.5, 0.75 and 1 and combined them using multigraph Leiden clustering to call modules with leiden resolution = 3. Although we identify smaller modules, here we only report modules with at least 10 genes. We also ran the modules method on three published datasets, for which we ran the method with the same parameters on each dataset (k = 100, z = 4.5, resolution = 2.5), used genes with >5% sparsity for the COVID and brain datasets and genes with >10% sparsity for the Tabula Sapiens dataset, and report modules with 10 genes or more.

50Module enrichments for cell subtypes and brain regions were performed using the hypergeometric test by calculating whether cells with a module score above 1 s.d. from the mean score were significantly enriched in a specific subtype or region. For plotting scores against other modules, averaged modules scores to either the subtype by sample level (within the same major cell type) or at the sample level alone (across cell types) and calculated Pearson correlations and P values using the cor.test function of R. To create the module–module network across microglial and immune modules, we calculated module–module score Pearson correlations using the logged module scores at the aggregated subtype by sample level, using a one-sided test with Padj < 0.01 as a cut-off. To generate contour plots of module expression on a UMAP, we first smoothed cell-level expression on a 500 × 500 grid with a 2D Gaussian kernel (size = 25 × 25 and σ = 1) and then plot contours for smoothed values (0.1 to 0.8).

51Gene expression programs underlying both cell type identity and cellular activities were determined according to the consensus NMF (cNMF) analysis pipeline established previously using the default parameters. The number of components (K) to use for cNMF was determined on the basis of a diagnostic plot showing the stability of the solution and the Frobenius reconstruction error as a function of K. To reduce runtime and working memory requirements, the data were downsampled using the Seurat function subset with the default parameters. The data were downsampled to 200 cells per major cell type. For the cNMF analysis at the level of high-resolution cell types, the analysis was performed separately on excitatory neurons, inhibitory neurons and astrocytes. For these analyses, the data were downsampled to 2,000 cells per astrocyte subtype and 1,000 cells per excitatory and inhibitory neuron subtype. Statistical significance of the overlap between the top 200 genes of a gene expression program and cell type marker genes was computed using Fisher’s exact tests.

52The gene regulatory network analysis was performed using pySCENIC with the default parameters. To reduce runtime and working memory requirements, the data were downsampled to 2,000 cells per major cell type. For the SCENIC analysis at the level of high-resolution cell types, the analysis was performed separately on excitatory neurons, inhibitory neurons and astrocytes and the data were downsampled to 1,000 cells per high-resolution cell type. To identify the top cell-type-specific regulons, we calculated regulon specificity scores as described by previously and ranked the regulons based on their regulon specificity score. Finally, we calculated the activity of each regulon in each cell using the AddModuleScore function of Seurat. The calculation of regulon module scores for major cell types was performed on a random sample of 50% of the cells (676,537 cells). For the analysis at the level of high-resolution cell types, the regulon module scores were determined based on all the cells of a major cell type class. For the statistical analysis of differences in the activity of regulons between cell types, the average regulon module score per individual and major cell type or high-resolution cell type was computed, respectively. The statistical analyses comparing the regulon module score activity was performed using Prism 9 software.

53GABAergic and glutamatergic module scores across all neuronal cell types were determined on the basis of a set of GABAergic and glutamatergic neuron marker genes, respectively, using the AddModuleScore function of Seurat. The sets of GABAergic and glutamatergic neuron marker genes were determined based on the human multiple cortical areas SMART-seq dataset from the Allen Brain Institute (https://portal.brain-map.org/atlases-and-data/rnaseq/human-multiple-cortical-areas-smart-seq). We identified marker genes using the FindAllMarkers function from Seurat (Wilcoxon rank-sum test with Bonferroni correction for multiple testing; Padj < 0.05). We tested only genes that were detected in a minimum of 25% of the cells within the cell type (min.pct = 0.25) and that showed, on average, at least a 0.25-fold difference (log-scale) between the cells of the cell class of interest and all remaining cells (logfc.threshold = 0.25). To quantify the intermediate character of thalamic excitatory and inhibitory neurons, we first computed the average GABAergic and glutamatergic module score values for each neuronal cell type and for each individual. We then used the resulting data to determine the first principal component (PC1) scores (the coordinates of the individual observations on the first principal component axis) using the princomp function in R. The ridgeline plot showing the distribution of PC1 score for each neuronal cell type was generated using the ggplot2 package in R. To determine the association between the average glutamatergic and the GABAergic module score across neuronal cell types, we performed a simple linear regression analysis using Prism 9 software.

54Marker genes significantly upregulated in extratelencephalic projection neurons (exc. L5 ET) compared with near-projecting excitatory neurons in layers 5 and 6 (exc. L5/6 NP) were determined using the FindAllMarkers function from Seurat (Wilcoxon rank-sum test with Bonferroni correction for multiple testing; Padj < 0.05). We tested only genes that were detected in a minimum of 25% of the cells within the cell type (min.pct = 0.25) and that showed, on average, at least a 0.25-fold difference (log-scale) between exc. L5 ET cells and exc. L5/6 NP cells (logfc.threshold = 0.25). The exc. L5 ET module score was computed based on the identified marker genes using the AddModuleScore function of Seurat. To determine the exc. L5 ET module scores across excitatory and inhibitory neurons, cells were downsampled to 2,000 cells per high-resolution cell type.

55Cell–cell communication events were predicted using the LIgand-receptor ANalysis frAmework (LIANA) in R. Specifically, the ligand–receptor analysis was performed using liana_wrap. The methods included were CellPhoneDB, NATMI and SingleCellSignalR. liana_aggregate with the argument ‘aggregate_how’ set to ‘magnitude’ was run to find consensus ranks of different methods. Only interactions (ligand–receptor pairs) with a robust rank aggregation (RRA) score smaller than 0.05 (aggregate_rank < 0.05) were considered in downstream analyses. The interaction score of ligand–receptor pairs was calculated by applying −log10 transformation to the RRA score (aggregate_rank). To determine the number of interactions and the overlap of interactions between regions, liana_wrap was run on the pool of cells isolated from all individuals, with separate analyses conducted for each brain region. To determine cell–cell communication events that are brain-region specific, liana_wrap was run separately for each individual. We then used a linear mixed-effects model to evaluate the association between the interaction scores of individual ligand–receptor pairs obtained from LIANA and the respective brain region serving as the predictor variable. To account for potential confounding factors and individual variability, we included age, sex and post-mortem interval as covariates in the linear mixed-effects model. These variables were added as fixed effects to the model. Moreover, we included a random effect for the individual to capture the participant-specific variability in the data. Linear mixed-effects models were implemented using the R software packages lme4 and lmerTest. The lmer function from the lme4 package was employed to fit the models. To obtain P values for the fixed effects in these models, we used the lmerTest package, which incorporated Satterthwaite’s degrees of freedom approximation. To account for multiple hypothesis testing, the obtained P values were further adjusted using the Bonferroni method.

56For comparing the relative abundance of major cell types across brain regions, the fraction of cells of a major cell type class was computed relative to all the cells isolated from a region. For the statistical analysis of cell type composition differences between brain regions, we also computed the relative abundance of major cell type classes separately for each study participant. To this end, the fraction of cells of a major cell type class was computed relative to all the cells isolated from a brain region of an individual. At the level of high-resolution cell types or cell states, two distinct measures of relative abundance were computed. First, the relative abundance of each subtype of a major cell class was computed as the proportion of a subtype relative to all cells of the corresponding major cell class isolated from a brain region. Second, for the statistical analysis of differences between brain regions, the fraction of cells of a subtype was computed relative to all the cells isolated from a brain region of an individual. The statistical analyses comparing the relative abundance of major cell types and subtypes between brain regions was performed using Prism 9 software.

57We calculated compositional differences in individuals with AD versus individuals without AD (or AD dementia versus no dementia) by modelling the number of cells of a certain cell type or subtype in a specific sample (individual by region) relative to the total number of cells using a quasi-binomial regression model. We modelled AD status by binary ascertainment variables (cogdx 4–5, NIA-Reagan score 1–2, Braak stage 5–6 versus others, as well as any detected presence of NFTs, neuritic plaque or diffuse plaques in the region) while adjusting for brain region and sex. We used the emmeans package in R to assess the significance of the regression contrasts and used p.adjust with the fdr method to adjust P values. We modelled the effects of fraction of cells on cognitive performance in multiple domains with gaussian linear regression of cognitive performance on last visit versus the log10-transformed fraction of cells in the subtype or major cell type jointly with covariates for age, sex, APOE-ε4 and post-mortem interval, with false-discovery rate P-value correction (p.adjust in R). We compared the fractional abundances of pairs of neuronal subtypes between two regions using Kendall’s τ only in individuals with AD (NIA-Reagan score). Significance was assessed using beta regression (R library betareg) controlling for sex, APOE genotype, post-mortem interval and age of death, and we adjusted P values using p.adjust in R with the fdr method.

58We performed differential expression analyses with three separate methods: MAST, Nebula and Wilcoxon testing,. For all methods, we subset the tested genes to only genes present in over 20% of cells. For MAST and Nebula, we calculated and included in the regression the top 10 components of unwanted variation calculated using RUV on the pseudo-bulk-level data (individuals by regions). For these methods, we also included as covariates the individual’s sex, age of death and post-mortem interval, each cell’s counts per gene and number of captured genes and, where applicable, the high-resolution cell subtypes and the brain region. For Nebula we used a Poisson mixed-model on the counts data with an offset of the log10-transformed total counts per cell. For MAST and Wilcoxon, we normalized each cell to a total library size of 10,000 counts. We ran Wilcoxon tests on both the cell and individual levels. We adjusted P values for multiple testing in all cases by using the p.adjust function in R with the fdr method. For our final set of differential genes in each analysis, we took all genes that were significant (Padj > 0.05) and concordant in both the MAST and Nebula results. We also provide the results for Wilcoxon tests, but did not use these to determine concordant results as they do not correct for many covariates. We computed differential expression against five AD ascertainment variables: continuous measurements of NFT, plaq_n, and plaq_d measured in each region except the thalamus (excluded from these analyses) and binary cognitive impairment (cogdx no dementia = 1 and 2 versus AD dementia = 4 and 5) and NIA-Reagan score classifications (non-AD = 3 and 4 versus AD = 1 and 2). We provide differential expression results for each of the 14 major cell types (with T cells, CNS macrophages, and each vascular subtype separately) for all regions jointly and for each region separately. We also provide results for each of the excitatory subtypes either in its most prevalent region for EC, HC or TH subclasses, or across the neocortex for neocortical subtypes (Supplementary Table 9). We also computed DEGs for the interaction between pathologic diagnosis of AD and sex in each major cell type, both across all regions and in each region separately. For the glial energy metabolism analyses we recomputed all DEGs in glial glycolysis-associated modules separately (keeping all genes, with no cell percentage cut-off). Glycolysis pathway diagram is from the glycolysis and gluconeogenesis pathway from WikiPathways.

59Pathology-biased DEGs were based on neuritic plaque or NFT pathology measurements in each region and were computed in each major cell type across all regions and in each region (except for the TH). Genes were ordered by the residual between NFT effect size and predicted NFT effect size from a regression using plaque effect size and region. Genes were retained if they were consistently up (or down) in 3+ regions for either NFT or plaque but in fewer than 2 regions for the other pathology measurement (shared genes are genes found in 2+ cell types).

60We compared our differential expression results to results from seven different previously published studies,,,–. We compared the published DEGs both to: (1) cross-region DEGs calculated in each major cell type for individual-level AD status (NIA-Reagan score or clinical diagnosis of AD) and for quantitative measurements of AD pathology (neuritic plaques, diffuse plaques and NFTs); and (2) region-specific DEGs calculated in each major cell type and in endothelial cells, computed relative to pathologic diagnosis of AD (NIA-Reagan score, AD, 1–2; non-AD, 3–4). As some studies reported only the significant genes, we compared the log-transformed fold change estimates for our DEGs and reported DEGs by a Pearson correlation test.

61To assess the enrichment of upregulated, downregulated non-differentially expressed genes in each module, we first assigned each tested DEGs to its closest module by correlation to the module’s average expression profile. We then performed a hypergeometric enrichment test for the number of genes in a category (upregulated, downregulated, not differentially expressed) assigned to the module, against the total number of tested genes assigned to the module, the total number in the category and the total number tested and corrected P values using p.adjust (Benjamini–Hochberg). Enrichments of pathology-biased DEGs in modules were performed in the same manner.

62To partition neuronal DEGs into non-vulnerable and vulnerable-associated subclasses, we calculated each genes’ average expressions and differential expression effect sizes at the subtype level and compared these to the relative depletion of the subtypes. For each gene that was differentially expressed in late-AD (stratified by Braak stage, late AD, 5–6 versus non-AD or early-AD, 1–4) in at least 25% of all neuronal subtypes, we calculated the correlation of its average subtype expression in late-AD with each subtype’s compositional stability (log2[OR] in late AD) across excitatory subtypes, separating non-vulnerable-associated genes (correlation > 0.2) from vulnerable-associated genes (correlation < −0.2). We calculated functional enrichments on neuronal DEG partitions using the top 250 genes ordered by effect size in each category. We further separated DEGs with higher effect sizes in vulnerable subtypes from those with similar effect sizes across all neuronal subtypes by calculating the correlation of their differential effect sizes in each subtype with that subtype’s depletion (log2[OR] in late AD). To perform enrichments along the continuum of genes associated with vulnerability to non-vulnerability, we kept only genes with biased effect sizes (effect size correlation < −0.2) and binned them along the axis of expression correlation (window size 0.2 for at 0.05 intervals) and performed functional enrichments for all bins jointly.

63We performed DEG enrichments for each differential expression run using the gprofiler2 package in R, with multi-query for upregulated and downregulated genes, as unordered queries, a P-value cutoff of 0.05, and using GO, REAC, WP, KEGG and CORUM as annotation sources, and retained enriched terms with fewer than 500 genes. Module and module cluster enrichments were performed in the same manner, using the core genes identified for each module and for genes found in more than two modules within a module cluster.

64We identified markers associated with excitatory neuron subtype vulnerability by performing linear regression to predict the log2[OR] of each subtype’s depletion in late AD based on its expression at the subtype-aggregate level (log[X + 1], averaged normalized expression in each subtype by region by individual batch), controlling for age, sex and post-mortem interval and adjusting P values with p.adjust (fdr).

65Processed snRNA-seq data (DLPFC, experiment 2) were obtained from Synapse (syn51123521) and integrated with our own PFC snRNA-seq dataset comprising 427 individuals. To identify vulnerable inhibitory neuron subtypes, we examined the association between the relative abundance of cell types and the measure of NFT density (variable tangles). We used a quasi-binomial regression model to model the number of cells belonging to a specific cell type in a given sample (individual study participant), relative to the total number of cells in that sample. We fitted the regression model using the glm function in R, including age, sex and post-mortem interval as covariates. P values were corrected for multiple testing using the Benjamini–Hochberg procedure as implemented in the R function p.adjust. The results are presented in the form of association scores (signed −log10-transformed Benjamini–Hochberg-adjusted P value, where the sign was determined by the direction (positive or negative) of the association). Inhibitory neuron subtypes demonstrating a significant negative association with tangle density (Benjamini–Hochberg-adjusted P value < 0.05) were classified as vulnerable subtypes, while all other subtypes were categorized as non-vulnerable. Genes exhibiting differential expression between vulnerable and non-vulnerable inhibitory neuron subtypes were identified on the basis of our PFC snRNA-seq dataset. This analysis was restricted to individuals without a pathologic diagnosis of AD. The differential expression analysis comparing vulnerable to non-vulnerable inhibitory neuron subtypes was performed using the R package dreamlet (https://diseaseneurogenomics.github.io/dreamlet/). We used the dreamletCompareClusters function with the argument ‘method’ set to ‘fixed’ for this analysis. Adjusted P values for multiple testing were obtained using the topTable function of dreamlet, with the ‘adjust.method’ set to ‘BH’.

66Intersection of regional expression and pathology-specific DEGs (across all regions) was performed for 149 annotated AD GWAS familial and AD risk loci from recent GWAS,–. We calculated the disease-relevance score of each cell in the dataset against a recent Alzheimer’s GWAS, using scDRS (based on MAGMA),,. For the scDRS results, we counted the fraction of cells with significant scDRS scores (FDR < 0.05) in each cell type, subtype and region. To test for overlap with microglia/immune modules, we compared the set of immune cells with significant expression of each module (z score > 2.5) and with the set of cells with significant scDRS scores (FDR < 0.05) and performed a hypergeometric test for significance of the overlap (Padj < 0.01, Benjamini–Hochberg correction). To identify region-specific GWAS genes, we performed an analysis of variance for the effect of region on average gene expression at the pseudobulk level.

67To quantify CR, we computed a CR score as the difference between the observed cognition and the cognition predicted by a linear regression model, given the level of pathology (Fig. 5a). Using this approach, we computed cognitive resilience (CR) scores based on the measure of global cognitive function and CDR scores based on the measure estimating the rate of change of global cognitive function over time (Fig. 5a). Four distinct CR and CDR scores were derived using four distinct measures of AD pathology, namely global AD pathology and, separately, neuritic plaque burden, NFT burden and tangle density.

68We performed differential expression analyses using the R package muscat to identify genes associated with CR in the PFC. Low-expressed genes were excluded and only genes with more than one count in at least ten cells were considered. To take advantage of robust bulk RNA-seq differential expression frameworks, such as edgeR, in a first step, muscat aggregates measurements for each sample (in each cluster) to obtain pseudobulk data. Using this approach, single-cell measurements were aggregated per study participant and cell type using the sum of raw counts option. Differential expression analysis was run using the edgeR method as implemented in muscat. We included as covariates the individual’s age at death and post-mortem interval. We report adjusted P values for multiple testing in all cases by using the p.adjust function with the Benjamini–Hochberg method as implemented in muscat. The multiple testing correction was performed locally, that is, on each of the cell types separately with the number of tests equal to the number of genes considered. These differential expression analyses were performed on the entire set of 427 individuals except for the group-based differential expression analysis based on our categorical definition of CR. In this case, we focused on comparing two distinct groups determined by their pathologic and clinical diagnoses of AD. First, we identified individuals with a pathologic diagnosis of AD, using the NIA-Reagan pathology criteria. Subsequently, these individuals were further categorized on the basis of their clinical consensus diagnosis of cognitive status at the time of death. Specifically, we compared individuals with no cognitive impairment (NCI, final consensus cognitive diagnosis (cogdx) value of 1) against individuals with a cognitive diagnosis of AD dementia and no other cause of cognitive impairment (cogdx value of 4) among individuals with a pathologic diagnosis of AD.

69To confirm the differential gene expression results based on the CR and CDR scores, we also evaluated the association between gene expression and global cognitive function or the rate of change of global cognitive function adjusting for AD pathology as a covariate. The AD pathology measures considered as a covariate were global AD pathology (gpath), neuritic plaque burden (plaq_n), NFT burden (nft), or tangle density (tangles). Thus, together with the DGE analysis based on CR and CDR scores, we performed a total of 16 different tests assessing the association between gene expression and CR.

70We used the model-based analysis of single-cell transcriptomics (MAST) tool to investigate whether the CR genes identified in PFC astrocytes were also associated with CR in astrocytes from other regions of the human brain. To ensure robust analysis, we initially filtered the genes under investigation, selecting only those with more than one count in at least 10 cells. The analytical model incorporated the condition variable of interest, as well as several covariates known to influence gene expression. These covariates included the cellular detection rate (cngeneson), age at death (age_death), post-mortem interval (pmi), and sex (msex). We also accounted for potential participant-specific variation in the data by incorporating a random effect term for the individual (1|individual). To account for multiple comparisons, the P values were adjusted using the FDR method as implemented in the p.adjust function.

71Differential expression analysis of bulk RNA-seq data from the ROSMAP cohorts was performed using DESeq2 (plotted) and edgeR. Age at death and post-mortem interval were converted into z scores and included as covariates in the regression equation. Both approaches (DESeq2 and edgeR) provided similar results.

72The average expression level of each gene within each major cell type was determined using the ‘AverageExpression’ function from the Seurat R package. The genes considered in the differential expression analysis for each major cell type were categorized into ten subsets based on their average expression level within the corresponding cell type. We next intersected the genes in each of the ten subsets with genes identified as significantly associated with either neuritic plaque burden (in our dataset) or the CPS score (in the MTG SEA-AD dataset). This intersection enabled us to determine the number of significant DEGs in each subset. The process was performed separately for genes positively and negatively associated with AD pathology. Subsequently, we randomly sampled the determined number of significant DEGs from each of the 10 subsets, ensuring that the expression level distribution of the DEGs was preserved in the random samples. This random sampling step was repeated for a total of 1,000 iterations. These steps were performed separately for both our dataset and the SEA-AD dataset. For each of the 1,000 random samples, we determined the overlap of genes between datasets and compared it to the observed overlap between the two datasets. To assess the significance of the observed overlap, we computed z scores, which represent the difference between the observed value of overlap and the mean value of overlap based on the permutation results, divided by the s.d. of the permutation results.

73To further validate our differential expression results, we evaluated the correlation between the effect sizes of gene expression changes observed in our study and those identified through quantitative proteomics. We specifically examined the correlation between the effect sizes of genes associated with neuritic plaque burden in our study and the effect sizes of overlapping differentially expressed proteins in the quantitative proteomics analysis of bulk tissue. The correlation was computed using the cor.test function in R with the argument ‘alternative’ set to ‘two.sided’ and the argument ‘method’ set to ‘pearson’. P values were adjusted for multiple testing using the Benjamini–Hochberg method as implemented in the R function p.adjust.

74We determined genes significantly associated with global AD pathology for each glial cell type, using single-nucleus RNA sequencing data from the PFC of 427 participants in the ROSMAP study. We then calculated module scores for these gene sets in astrocytes, microglia, oligodendrocytes, and OPCs using the Seurat ‘AddModuleScore’ function. The module scores were determined separately for genes positively and negatively associated with global AD pathology. To assess the progression of these scores across the spectrum of global AD pathology burden, we averaged the module scores of all cells of a specified cell type isolated from the brain region of interest of an individual. For visualizing the relationship between the global AD pathology burden and mean module scores, we employed Locally Estimated Scatterplot Smoothing (LOESS) using the ggplot2 package in R, with the ‘geom_smooth’ function and the ‘method’ parameter set to ‘loess’. The correlation of mean module scores between regions was determined using the cor.mtest function of the R package corrplot. P values were adjusted for multiple hypotheses testing using the Benjamini–Hochberg method as implemented in the R function p.adjust.

75Frozen human post-mortem brain samples were embedded in Tissue-Tek OCT compound (VWR; 25608-930), sectioned on a Leica cryostat at a thickness of 20 μm, and mounted onto Fisherbrand Superfrost Plus microscope slides (Thermo Fisher Scientific; 12-550-15). Slides were fixed in 4% paraformaldehyde at 4 °C for 30 min, and dehydrated in ethanol. The RNAscope 2.5 HD Chromogenic Duplex Detection Kit and RNAscope Multiplex fluorescent V2 Kit (ACDBio) were then used according to the manufacturer’s instructions. Tissue samples were hybridized using the following chromogenic RNAscope probes: GAD2, FOXP2, MEIS2, AQP4, GRM3, ADCY8, PFKP, PNPLA6, GPCPD1 and CHDH (ACDBio). For in situ hybridization of Reelin, tissue samples were hybridized using the following fluorescent RNAscope probes: vGlut and Reelin. Cell nuclei were stained with 50% haematoxylin (for chromogenic experiments) or with Hoechst (for fluorescent experiments). For fluorescence RNAscope analysis, sections were incubated in TrueBlack (Biotium; 23007) for 10 s before Hoechst staining to quench auto-fluorescence. Images were acquired using the Zeiss LSM 900 confocal microscope, with a 63× oil objective. Two images were acquired per tissue sample. For both chromogenic and fluorescence RNAscope experiments, puncta were manually counted by researchers blinded to the experimental group of each image.

76All experiments were performed according to the Guide for the Care and Use of Laboratory Animals and were approved by the National Institute of Health and the Committee on Animal Care at Massachusetts Institute of Technology. Sample sizes were determined on the basis of previous work from our laboratory, without power analysis calculation or randomization. In the experiment comparing App-KI (C57BL/6-App<tm3(NL-G-F)Tcs>, RBRC06344) and WT mice, the App-KI group consisted of 7 mice (5 male and 2 female) and the WT group included 6 female mice. In the experiment comparing Tau P301S (The Jackson Laboratory, 008169) to WT mice, both groups consisted entirely of male mice, with 6 mice in the Tau(P301S) group and 5 mice in the WT group. Mice were transcardially perfused with ice-cold phosphate-buffered saline, followed by 4% paraformaldehyde for fixation. Brains were dissected out and post-fixed in 4% paraformaldehyde overnight at 4 °C. Brains were sectioned horizontally on the Leica vibratome at a thickness of 40 μm. Slices containing the EC were selected under a dissecting microscope to ensure consistent anatomical structure across all cohorts. Brain sections were incubated in antigen retrieval solution (pH 6, 100 mM sodium citrate buffer, prewarmed to 80 °C) for 20 min, and then cooled to room temperature. The sections were then washed twice with phosphate-buffered saline, and blocked in buffer (0.3% Triton X-100, 2% bovine serum albumin, 10% normal donkey serum in phosphate-buffered saline) for 1 h at room temperature.

77The sections were incubated in primary antibodies (anti-Reelin, 1:200; anti-NeuN, 1:200, anti-phospho tau, 1:200; and anti-amyloid-β, 1:500) overnight at 4 °C. After primary antibody incubation, the sections were washed three times with PBS, twice with blocking buffer and incubated in secondary antibody (1:1,000) for 2 h at room temperature. The sections were then washed three times with PBS, incubated with Hoechst (1:1,000) for 10 min and washed once more with PBS.

78Confocal tile scans of the EC were acquired on the Zeiss LSM 900 using a 20× objective, with consistent laser setting across all cohorts. Layer II–III EC was identified based on previous criteria. Orthogonal projections of the confocal tile scans were exported to Fiji for signal quantification. In Fiji, layer II–III of the EC was set as a region of interest, and a macro was used to count Reelin-positive cells in the region of interest and quantify the mean fluorescence intensity for each cell. The signal intensity of the Reelin channel was subsequently normalized to the NeuN signal. Researchers were blinded to animal genotype.

79Processed snRNA-seq data generated by the Seattle Alzheimer’s Disease Brain Cell Atlas (SEA-AD) consortium (SEAAD_MTG_RNAseq_final-nuclei.2023-05-05.h5ad) were obtained from the Seattle Alzheimer’s Disease Brain Cell Atlas (SEA-AD) (https://sea-ad-single-cell-profiling.s3.amazonaws.com/index.html#MTG/RNAseq/). The SEA-AD DLPFC data (SEAAD_DLPFC_RNAseq_final-nuclei.2023-07-19.h5ad) were downloaded from https://sea-ad-single-cell-profiling.s3.amazonaws.com/index.html#DLPFC/RNAseq/. Additional processed snRNA-seq datasets (specifically the h5ad files Neurons.h5ad and Nonneurons.h5ad) were obtained from the Linnarsson laboratory (https://console.cloud.google.com/storage/browser/linnarsson-lab-human;tab=objects?authuser=0&prefix=&forceOnObjectsSortingFiltering=false).

80Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

81Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41586-024-07606-7.