Human microglial transitions at the Aβ–tau inflection point associate with divergent pathways to dementia and resilience

Postmortem human brain tissues

We selected a total of 56 individuals from the NBB (Netherlands Institute for Neuroscience) or from the Dutch 100-plus Study⁴⁵, all with an RNA integrity number (RIN) above 6.8. All donors provided a written informed consent for brain autopsy and research use of material and clinical data following the guidelines from the ethics committee of the VU University Medical Center and the KU Leuven ethics committee (study number S63259). Autopsy procedures were approved by the medical ethics committee of the VU University Medical Center. The study included cognitive assessments (Clinical Dementia Rating (CDR), Reisberg scale or MMSE), neuropathological evaluations, APOE genotyping and bulk tissue transcriptomic profiling⁹. Detailed clinicopathological information is provided in Supplementary Table 1.

We selected eight controls (OCT_HC) with no/minimal AD pathology and 16 individuals with mild-to-severe Aβ and tangle deposition (OCT + DEM and OCT − DEM) (Supplementary Table 1). OCT_HC and OCT − DEM individuals were cognitively normal, whereas OCT + DEM had cognitive impairment. Donors with non-AD psychiatric or neurological disease were excluded. A second cohort (CEN) from the 100-plus Study comprised centenarians with mild-to-severe Aβ but minimal tangle pathology and a cognitive spectrum from normal to impaired, as assessed by MMSE. Cohorts were balanced for sex and age (Supplementary Table 1). Frozen superior frontal gyrus cylinders (8-mm diameter, covering all six cortical layers and white matter) were collected via biopsy punch (Ted Pella, 15111-80), shipped on dry ice to the De Strooper laboratory in KU Leuven and stored at −70 °C. Cohorts were balanced for sex; however, analyses were not sex stratified due to limited subgroup sizes, as the aim was to identify shared pathology-associated cellular trajectories.

Pathological features and scoring

Aβ and pTau pathology were assessed using NBB-provided neuropathological scores (Thal amyloid phase, Braak stages for Lewy bodies and neurofibrillary tangles, CERAD neuritic plaque scores and ABC scores)⁶⁰ (Extended Data Fig. 1a and Supplementary Table 1). Cored plaques were defined as cortical amyloid plaques with a compact central core and diffuse Aβ halo, and diffuse plaques were sharply delineated extracellular Aβ aggregates without an amyloid core⁵⁹. All other plaque types with ill-defined borders or linked to specific anatomical areas such as subpial or white matter Aβ aggregates⁵⁹ were considered as unclassified plaques in the context of this study. Aβ plaque density was significantly higher in OCT + DEM, OCT − DEM and CEN compared to OCT_HC (Fig. 1d), with no significant differences among the three pathology groups, indicating similar frontal lobe plaque loads. By contrast, OCT + DEM cases had significantly more neuritic plaques than OCT − DEM, OCT_HC and CEN (Extended Data Fig. 1a). Because Thal and Braak scores reflect brain-wide pathological distribution rather than individual Visium spot loads, Aβ and pTau abundance and morphology were determined per Visium spot via immunohistochemistry on the same tissue sections used for spatial transcriptomics (Fig. 1e).

Tau morphology across experimental groups

pTau pathology was assessed using an AT8 antibody against pTau (Ser202/Thr205; Pierce, 1/100)^61,62 and revealed large qualitative differences across experimental groups. In OCT + DEM cases, pTau regions of interest (ROIs) exhibited neuronal threads, neurofibrillary tangles and dystrophic neurites strongly associated with Aβ plaques. Conversely, pTau in OCT_HC brains was diffusely distributed and predominantly localized to the soma. OCT − DEM and CEN brains showed pTau morphologies resembling those in OCT_HC, suggesting less advanced pathology. Age-related tau astrogliopathy (ARTAG), often observed in the white matter of cognitively healthy seniors, was also detected. Differentially expressed genes associated with pTau and neuritic plaques showed significant correlations within the same experimental group (for example, r = 0.5 for OCT + DEM) but differed substantially between groups (for example, r = 0.1 between OCT + DEM and OCT − DEM). This suggests that variations in pTau types contribute to differences in gene expression profiles across experimental groups.

Tissue collection

For Visium Spatial Gene Expression, frozen brain blocks were cryosectioned at 10 µm (Leica, CM3050S cryostat; chamber −16 °C to −18 °C, specimen holder −12 °C to −14 °C). One section per donor was mounted onto a Visium Spatial Gene Expression Slide (10x Genomics, PN-1000185). To minimize batch effects, sections from different donors were distributed across slides such that each slide combined at least one OCT + DEM, one OCT − DEM and one OCT_HC or two CEN + DEM and two CEN − DEM. Additionally, 7–10 adjacent sequential sections (10 µm) were collected onto SuperFrost Plus glass slides (Thermo Fisher Scientific). All slides were stored at −70 °C. For snRNA-seq, 15 sections (50 µm each) were collected per donor from the same frozen blocks and stored at −70 °C. The tissue distance between Visium and snRNA-seq sections ranged from 0.1 mm to 1 mm across donors.

Nuclei isolation and single-nuclei RNA library preparation

Single nuclei were extracted from frozen tissue sections. In brief, OCT was removed by incubating sections in 1 ml of cold salt-Tris solution (146 mM NaCl, 10 mM Tris (pH 7.5), 1 mM CaCl₂, 21 mM MgCl₂, 1 mM 2-mercaptoethanol, 1× cOmplete protease inhibitor, 0.2 U µl⁻¹ RNasin Plus) for 1 minute on ice, followed by centrifugation (500g, 3 minutes, 4 °C). The pellet was resuspended in 500 µl of ice-cold homogenization buffer (salt-Tris base supplemented with 25 mM KCl, 0.03% Tween 20, 0.01% BSA, 250 mM sucrose, 0.5 U µl⁻¹ RNasin Plus) and homogenized in a KIMBLE Dounce grinder (Sigma-Aldrich) with 10 strokes of pestle A and 5–10 strokes of pestle B. The homogenate was filtered through a 70-µm strainer (Greiner Bio-One), incubated 5 minutes on ice and pelleted (500g, 5 minutes, 4 °C). The pellet was resuspended in 2.65 ml of homogenization buffer without Tween 20. All buffers and equipment were prechilled.

Nuclei were isolated via OptiPrep density gradient⁶³. Homogenate was mixed with equal volume gradient medium (50% OptiPrep, 1 mM CaCl₂, 5 mM MgCl₂, 10 mM Tris (pH 7.5), 75 mM sucrose, 1 mM 2-mercaptoethanol, 0.5× cOmplete protease inhibitor, 0.5 U µl⁻¹ RNasin Plus), layered onto 4 ml of 29% gradient medium and centrifuged (10,160g, 30 minutes, 4 °C; Beckman Coulter, SW41 Ti rotor). The nuclei pellet was resuspended in PBS with 1% BSA and 1 U µl⁻¹ RNasin Plus and filtered through a 40-µm Flowmi strainer (Sigma-Aldrich). Final volume was 100 µl per sample.

Nuclei count and viability were assessed using AO/PI staining on a LUNA-FL counter (Westburg). We targeted 5,000–6,000 nuclei per sample, pooled per three samples and loaded onto a Chromium Next GEM Single-Cell 3′ Chip (10x Genomics). Library preparation followed the v.3.1 protocol (CG000204, Rev D). After quality control on an Agilent Bioanalyzer 2100, libraries were sequenced by BGI Tech Solutions using DNBSEQ paired-end sequencing.

DAB staining

Adjacent cryosections (10 µm) were used to quantify local Aβ plaque and pTau load. Sections were post-fixed in 4% PFA/PBS for 10 minutes. For Aβ staining, antigen retrieval was performed by boiling in citrate buffer (0.01 M (pH 6.0), 700 W microwave, 10 minutes) followed by 70% formic acid (10 minutes). Sections were blocked with 5% fetal calf serum (30 minutes, room temperature). Primary antibodies (Thermo Fisher Scientific, anti-pTau AT8, 1:2,000; Signet, anti-amyloid 4G8, 1:10,000) were incubated overnight at 4 °C. Secondary antibodies (Dako, horse anti-mouse HRP, 1:400, or Dako, EnVision Detection System, K8023) were applied for 1 hour at room temperature, followed by ABC (Vector Laboratories) and DAB development (Dako). Sections were scanned at ×20 on an Axio Scan.Z1 (Zeiss). Two gray matter ROIs were selected per section. For quantification, a threshold of 3× background optical density was applied. The corrected integrated optical density (cIOD) was calculated as the optical density of positive area multiplied by its surface area, divided by total ROI area.

Immunofluorescence, image acquisition and Visium Spatial Gene Expression

To detect Aβ, pTau and neuronal nuclei on the tissue sections of the OCT cohort used for Visium, we modified the immunofluorescence staining protocol provided by 10x Genomics (CG000312, revision version A). In brief, after methanol fixation at −20 °C for 30-minute and 5-minute incubation with blocking buffer (3 × SSC supplemented with 2% w/v BSA, 0.1% Triton X, 1 U ml⁻¹ RNase inhibitor, TruStain FcX (1:35) and 20 mM ribonucleoside vanadyl complex) at room temperature, tissue sections were incubated consecutively with four different antibody solutions, prepared in blocking buffer. Tissues were first incubated with 10 µg ml⁻¹ anti-Aβ 17-24 (4G8) antibody (BioLegend, 800701) and then with 10 µg ml⁻¹ donkey anti-mouse IgG (H + L) Alexa Fluor 555 antibody (Thermo Fisher Scientific, A-31570). The third incubation step was a mixture of 30 µg ml⁻¹ Alexa Fluor 647-conjugated anti-RBFOX3/NeuN antibody (1B7; Novus Biologicals, NBP1-92693AF647) and 3.33 µg ml⁻¹ biotin-conjugated anti-pTau antibody (AT8) (Thermo Fisher Scientific, MN1020B), followed by a final incubation with a mixture of 3.4 µg ml⁻¹ DAPI and 10 µg ml⁻¹ Alexa Fluor 488-conjugated streptavidin. Each incubation lasted for 30 minutes at room temperature. After each incubation, five sequential wash steps with wash buffer (3 × SSC supplemented with 2% w/v BSA, 0.1% Triton X, 1 U ml⁻¹ RNase inhibitor and ribonucleoside vanadyl complex) were applied as suggested by the company. To quench lipofuscin autofluorescence, the sections were treated for 30 seconds with 1× TrueBlack solution (Biotium, 23007), diluted in 70% ethanol, after which the slide was rinsed by dipping 15 times in 3 × SSC buffer. Finally, 200 µl of 80% glycerol with RiboLock RNase inhibitor 2 U µl⁻¹ (Thermo Fisher Scientific, EO0382) was added. A coverslip was applied before imaging. To reduce the staining procedure time for the CEN cohort, we formed a pre-complex by incubating 3.33 µg ml⁻¹ biotin-conjugated anti-pTau antibody (AT8) (Thermo Fisher Scientific, MN1020B) and 6.67 µg ml⁻¹ Alexa Fluor 488-conjugated streptavidin for 30 minutes at room temperature prior to adding it to the tissue sections, in combination with 30 µg ml⁻¹ Alexa Fluor 647-conjugated anti-RBFOX3/NeuN antibody (1B7) (Novus Biologicals, NBP1-92693AF647) and 3.4 µg ml⁻¹ DAPI.

Next, 16-bit fluorescent images of DAPI, AT8, 4G8 and NeuN were acquired on an Axio Scan.Z1 slide scanner (Hamamatsu Orca Flash 4.0 V3 camera, Plan APO ×20/0.8 numerical aperture objective, ZEN Blue v.3.1). Channels were excited at 353 nm, 493 nm, 553 nm and 653 nm and collected with BP335-470-nm, BP500-550-nm, BP570-640-nm and BP665-715-nm emission filters, respectively. Images covered the entire Visium fiducial frame (8 × 8 mm) for downstream registration. Coverslips were removed by immersing slides at 45° into 3 × SSC buffer, followed by an additional wash. The entire procedure from methanol fixation through coverslip removal was kept under 5 hours. RiboLock RNase inhibitor (2 U µl⁻¹; Thermo Fisher Scientific, EO0382) and/or ribonucleoside vanadyl complex (20 mM; New England Biolabs, S1402S) were included throughout to minimize RNA degradation.

After coverslip removal, tissues were immediately processed using the Visium Spatial Gene Expression Kit (10x Genomics, PN-1000186) and the Library Construction Kit (10x Genomics, PN-1000190), following the demonstrated protocol (CG000239, Rev D for OCT or Rev C for CEN). Tissues were permeabilized for 18 minutes (OCT) or 5 minutes (CEN) to release poly-adenylated mRNA onto Visium spots, followed by spatially barcoded reverse transcription. Second-strand synthesis was performed on-slide, after which cDNA was transferred to tubes for polymerase chain reaction (PCR) amplification and quality control on an Agilent Bioanalyzer 2100 High Sensitivity chip. Libraries were constructed through enzymatic fragmentation, size selection, end repair, A-tailing, adaptor ligation and sample index PCR (annealing at 67 °C for OCT/Rev D or 54 °C for CEN/Rev C, with recommended cycle numbers). Final libraries were quality controlled on a Bioanalyzer High Sensitivity chip and sequenced via paired-end DNBSEQ technology (BGI Tech Solutions).

Xenium in situ gene expression, immunofluorescence and image acquisition

Xenium in situ gene expression was performed at LISCO (Leuven, Belgium) on five OCT + DEM and six OCT − DEM samples. Tissue sections on Xenium slides were fixed and permeabilized using Xenium Sample Prep reagents (10x Genomics, PN-1000460) following the demonstrated protocol (CG000581). In situ gene expression was performed using a standalone Xenium Custom Gene Expression Panel following protocols for probe hybridization, ligation and amplification (CG000582) and Xenium Analyzer operation (CG000584). In brief, fixed and permeabilized sections were incubated overnight with the custom probe panel, followed by post-hybridization washes, ligation to generate circular DNA probes and rolling circle amplification. An autofluorescence quencher and DAPI nuclear stain were applied before loading Xenium cassettes into the analyzer for imaging and decoding (consumables PN-1000487).

After the Xenium run, sections were stored in 750 µl of 1× PBS at 4 °C (maximum 1 week) before downstream immunofluorescence for Aβ, pTau and NeuN. The autofluorescence quencher was removed with 10 mM sodium hydrosulfite (10 minutes, room temperature) and 3× water rinses. Sections were blocked (1× PBS, 5% normal donkey serum, 0.5% Triton X; 90 minutes, room temperature) and sequentially incubated with (1) 5 µg ml⁻¹ anti-Aβ 4G8 (BioLegend, overnight, 4 °C); (2) 5 µg ml⁻¹ donkey anti-mouse Alexa Fluor 555 (90 minutes, room temperature); (3) 15 µg ml⁻¹ Alexa Fluor 647-conjugated anti-NeuN mixed with 2 µg ml⁻¹ biotin-conjugated AT8 (2 hours, room temperature); and (4) 1 µg ml⁻¹ DAPI mixed with 10 µg ml⁻¹ Alexa Fluor 488 streptavidin (1 hour, room temperature). Sections were washed three times for 10 minutes with wash buffer (1× PBS, 0.5% Triton X) between incubations. Lipofuscin autofluorescence was quenched with 1× TrueBlack in 70% ethanol (30 seconds), rinsed in 3 × SSC and mounted with 20% SlowFade Diamond (Thermo Fisher Scientific, S36967) in 80% glycerol.

Then, 12-bit fluorescent images were acquired on a Nikon Ni-E Eclipse (Digital Sight 10 camera, Plan APO ×20/0.8 numerical aperture, NIS-Elements v.3.1). DAPI, AT8, 4G8 and NeuN were excited at 400 nm, 470 nm, 555 nm and 635 nm (CoolLED, pE800) and collected with BP335-470-nm, BP500-550-nm, BP570-640-nm and BP665-715-nm emission filters, respectively.

Image analysis

For Visium, 16-bit fluorescence images were aligned with the sequenced Visium array using the fiducial frame. The 555-channel contrast was adjusted to visualize fiducial markers and tissue boundaries; images were rotated to correct orientation; and the channel was exported as an 8-bit grayscale single-page TIFF at full resolution. These TIFFs were uploaded to Loupe Browser 5.0 (10x Genomics) for manual fiducial alignment, after which Visium spot coordinates were automatically calculated based on the aligned frame and slide ID. Tissue boundaries were indicated, and an aligned JSON file containing tissue-overlapping spot IDs and pixel coordinates was generated and passed to Space Ranger to compute per-spot transcriptomic profiles.

For Xenium, fluorescent pathology images (Aβ, pTau, NeuN and DAPI) were aligned to the Xenium DAPI image in Fiji⁶⁴ using the ‘register virtual stack slices’ plugin (affine model). The higher-resolution Xenium DAPI was downscaled to match the pixel size of the pathology images before serving as alignment reference.

Pathology annotations were performed in QuPath (v.0.2.3)⁶⁵ on aligned fluorescent images. A tissue mask was manually annotated to remove artifacts and subdivided into white matter and cortical layers based on NeuN and DAPI morphology. Vessels and pia were annotated separately to account for non-specific antibody staining. Aβ plaques and pTau pathology were annotated using pixel intensity thresholds and size filters. Plaques were further classified as core plaques, neuritic plaques or diffuse plaques based on neuropathological and morphological criteria: neuritic plaques were characterized by dystrophic neurites surrounding plaques, core plaques by circled-dot morphology (Fig. 1b) and diffuse plaques by lower Aβ intensity and heterogeneity. Diffuse plaques were distinguished from neuritic plaques/core plaques using the 15th percentile of ranked pixel intensity standard deviations. Annotations overlapping vessels and pia were removed.

For Visium, spot boundaries were computed, and fluorescent channel intensity measures (mean, median and standard deviation for Aβ, pTau, NeuN and DAPI) were recorded per spot. The resulting metadata matrix detailing ROI overlap per spot was used for downstream transcriptomic analysis.

Visium spot selection and analysis

Unbiased clustering of Visium spots by transcriptomic profiles revealed patterns consistent with cortical layer identity, individual variability and pathological heterogeneity (Fig. 1f). To relate local transcriptomic variation to protein aggregation, we quantified the proportion of each spot overlapping with Aβ and pTau pathology. As pTau accumulation was predominantly observed in cortical layers 3–5, particularly in OCT + DEM cases (Fig. 1c), subsequent analyses focused on these layers.

Centenarian cohort analysis

The centenarian cohort was not subdivided by dementia status due to insufficient separation in clinical data. CEN brains showed pTau and plaque morphologies with no significant differences in Aβ plaque density compared to OCT ± DEM cases.

Human FFPE brain samples

FFPE brain samples from three human autopsy cases were analyzed (clinical and pathological details in Supplementary Table 1). Individuals were selected based on absence of encephalitis/meningitis and presence of symptomatic AD, defined by clinical dementia and intermediate-to-high AD neuropathological changes. Autopsies were performed at UZ Leuven (Belgium) in compliance with national regulations, with written informed consent and ethical approval from the UZ Leuven Ethics Committee. The right hemisphere was dissected for gross neuropathological examination and stored at −80 °C. The left hemisphere was fixed in 4% formaldehyde for 2–4 weeks and then dissected and macroscopically evaluated. Frontal cortex samples were collected for further analyses. CDR scores⁶⁶ were retrospectively assigned from medical records.

Left hemisphere tissue blocks were paraffin embedded and sectioned at 5 µm onto Flex IHC adhesive slides (Dako), dried at 55 °C and stained with hematoxylin and eosin for diagnostic evaluation. AD pathology was assessed via immunohistochemistry using established criteria: Aβ plaque deposition (Aβ phase)⁶⁷, tau neurofibrillary tangle distribution (Braak neurofibrillary tangle stage)⁶⁸ and neuritic plaque density (CERAD score)⁶⁹. The NIA−AA score⁶⁰ was determined from these three indices.

Immunofluorescence for HLA and SPP1

Immunofluorescence was performed on frontal cortex FFPE AD and frozen centenarian brain sections. FFPE sections were deparaffinized and subjected to heat-induced epitope retrieval (pH 6.1). Frozen sections were heated briefly (1 minute, 37 °C) and fixed in 4% PFA (20 minutes). Amyloid plaques were stained with X34 (20 minutes, room temperature). Frozen sections were blocked and permeabilized (5% BSA, 0.2% Triton X-100 in PBS). Primary antibodies were applied overnight at room temperature, followed by fluorophore-conjugated secondary antibodies or streptavidin. IBA1 served as microglial marker for FFPE tissue and CD45 for frozen tissue. Autofluorescence was quenched with TrueBlack (30 seconds; Biotium). Slides were mounted in Glycergel (Agilent Technologies). Antibody details are listed in Supplementary Table 10. Images were acquired on a Nikon A1R confocal system coupled to a Nikon Eclipse Ti inverted microscope using NIS-Elements. Image processing and figure assembly were performed in ImageJ (National Institutes of Health) and Inkscape (https://inkscape.org/).

snRNA-seq processing

Raw FASTQ files were processed using Cell Ranger v.6.0.1 (10x Genomics) with default settings, aligning to the human GRCh38 genome (prebuilt reference 2020-A) for filtering, barcode counting and unique molecular identifier (UMI) quantification.

Because snRNA libraries contained pooled individuals, feature−barcode matrices were demultiplexed using Souporcell (v.2.0)⁷⁰ based on genetic profiles from the Illumina Global Screening Array (GSA). Cell Ranger BAM files were used to calculate allele frequency tables per library. Souporcell identified variants, computed allele frequencies and clustered cells by individual. Clusters were matched to reference genotype profiles via Pearson correlation, with assignments made at greater than 90% similarity. Sample IDs were added as metadata to the feature−barcode matrix.

Count matrices were imported into Scanpy (v.1.8.2)⁷¹. Cells with fewer than 200 genes, genes in fewer than five cells and cells with more than 20% mitochondrial reads were filtered. Counts were normalized (total count normalization, scale factor 10,000) and log transformed.

Doublets were removed using Scrublet⁷² (threshold 0.39). Remaining cells were clustered using scVI (scvi-tools v.1.1.6)⁷³. Highly variable genes (HVGs) were selected accounting for sequencing batch (Scanpy highly_variable_genes: min_mean = 0.0125, max_mean = 3, min_disp = 0.5, batch_key = batch). This yielded 112,698 single nuclei with a mean of 3,399 genes per nucleus (IQR: 1,563–4,663) for the OCT cohort. An overview is provided in Supplementary Data 1.

Visium spatial transcriptomics processing

Raw Visium FASTQ files were processed using Space Ranger 1.2 (10x Genomics) with the GRCh38 reference (2020-A). The filtered count matrix was combined with spot positions and scaling factors for integration with Visium image analysis data. Inclusion criteria were as follows: ≥200 counts per gene, ≥10 genes per spot and ≥500 reads per spot. Staining artifacts, low-quality tissue regions, vessels and samples with poor RNA quality were removed. The curated dataset comprised 60,129 spots from 24 individuals (mean 1,831 genes, IQR: 1,168–2,352). The same approach yielded 61,739 spots from 20 centenarian individuals.

After removing spots with imaging artifacts or outside the tissue mask, the top 2,000 HVGs were selected. Samples were integrated using scVI (v.1.1.6)⁷³, batch correcting for RIN (continuous covariate) and gyrus frontalis superior (GFS) location and Visium slide (categorical covariates). The scVI latent representation was used to construct a k-nearest neighbor graph for UMAP⁷⁴ visualization (min_dist = 0.2) and Leiden clustering (resolution = 0.2).

Cell deconvolution

To map cell types and states from snRNA-seq onto Visium data, we used cell2location (v.0.1)⁴⁸. The snRNA-seq reference was trained using the RegressionModel (batch_size = 2,500, train_size = 1, lr = 0.002, max_epochs = 200), with sequencing lane ID as batch covariate, GFS as categorical covariate and RIN_NBB as continuous covariate. Reference cell state signatures were estimated by simulating 1,000 posterior samples with the same batch size. For spatial deconvolution, cell2location was run with N_cells_per_location = 8, detection_alpha = 200, max_epochs = 30,000, batch_size = None and train_size = 1, using Visium slide ID as batch covariate alongside GFS and RIN_NBB. Posterior samples were drawn using the same configuration as the reference model. In downstream analyses, the 5% quantile of the posterior distribution (q05_cell_abundance) was used as a conservative estimate of cell type abundance per Visium spot.

TD clustering and visualization

Clustering was performed on the cell2location cell type abundance matrix. Abundances were normalized to z-scores based on total cell abundance per spot. Mid-layer (layers 3–5) Visium spots from 24 OCT individuals were used. To mitigate layer-specific variation while preserving pathology-associated biology, Harmony integration⁴⁹ was applied on layer (Scanpy harmony_integrate, 20 principal components, max_iter_harmony = 20). The same procedure was applied to 20 CEN individuals. Leiden clustering was performed on cell type abundance profiles using a k-nearest neighbor graph (20 neighbors) in Harmony-integrated principal component analysis space (random seed 42, resolution 0.4), yielding six TDs for OCT (OCT-TD0 to OCT-TD5) and seven TDs for CEN (CEN-TD0 to CEN-TD6). Results were visualized using UMAP.

To address overplotting, hexbin plots (Matplotlib⁷⁵) were used, with hexbins containing fewer than five spots grayed out. Three coloring schemes were applied: (1) percentage score (fire colormap) for categorical variables, showing the proportion of spots associated with a given category; (2) average score for the aggregated mean of continuous variables per hexbin; and (3) mean difference score (blue/red colormap) reflecting the difference in z-normalized means of continuous variables or gene set module scores between OCT + DEM and OCT − DEM spots per hexbin.

Differential gene expression analysis

Differential gene expression analysis was performed on Visium data at different pathology and phenotype levels using edgeR⁷⁶. Pathology load was calculated as the log1p-transformed percentage of each spot covered by the specific pathology ROI. Spots from layers 3–5 (100% overlap) were extracted. A quasi-likelihood negative binomial generalized log-linear model was applied for genewise statistical tests with robust prior quasi-likelihood dispersion estimation.

Differential cell abundance analysis

Differential abundance analysis was performed using edgeR on the same layer 3–5 Visium spots used for differential gene expression analysis, leveraging negative binomial generalized linear model methods to model the overdispersed cell abundance matrix. Cell2location-predicted abundances were normalized by total cells per sample. The percentage overlap of pathology per spot, derived from image analysis, was included as continuous covariates in the design formula and tested using glmQLFTest for each phenotype separately. The analysis yields log₂ fold changes and FDR-corrected P values (Benjamini−Hochberg), enabling assessment of cell abundance changes in response to plaque and pTau pathology across phenotypes.

TD prediction in the CEN cohort

We trained an MLP from scikit-learn (http://jmlr.org/papers/v12/pedregosa11a.html) trained on data from the OCT cohort to predict the assigned OCT TD based on the underlying, normalized, cell type abundance matrix. We subsequently applied this predictor to the CEN cell2location abundance matrix to predict in what OCT TD each CEN Visium spot would fall.

Gene co-expression analysis of Visium Spatial Gene Expression

Signed co-expression networks were constructed using WGCNA (R, v.1.72.5)⁷⁷ on Visium spatial transcriptomic data from gray matter regions of the OCT cohort (17,271 genes, 41,943 spots) and the CEN cohort (18,090 genes, 39,358 spots) independently. Raw count matrices were normalized and variance stabilized using SCTransform (Seurat 4.0). A soft power of 14 was selected via pickSoftThreshold. Modules were identified using cutreeDynamic (deepSplit = 2) with a minimum of 30 genes per module, yielding 16 modules for the OCT cohort and 20 for the CEN cohort.

Xenium data processing

Xenium in situ sequencing was performed on 11 samples (six OCT − DEM and five OCT + DEM) using a panel of 329 genes, including cell type/subtype markers selected via scMAGS (v.1.5)⁷⁸ and genes of interest identified from Visium analysis. scMAGS selects markers highly expressed in a specific cell type with low expression in others, making them well suited for spatial transcriptomics. Raw Xenium data were preprocessed using Xenium Ranger (v.1.7). Cell segmentation was performed on DAPI nuclear staining across z-stack planes, generating non-overlapping two-dimensional nuclear masks projected onto the x−y plane. Masks were expanded by 5 µm or until intersecting another cell boundary, providing approximate cell segmentation. Transcripts within boundaries were assigned to cells, producing a per-cell count matrix.

For cell type identification, quality filters were applied: nuclear area 9.5–130 µm², transcript counts 13–900, cell area 11.92–444 µm² and ≥0.8 transcripts per µm². Doublets were removed using Scrublet (v.0.2.3, threshold 0.82), retaining 424,056 cells. Clustering was performed using Scanpy (v.1.9.1) and scVI (v.0.19.0). HVGs were identified based on mean expression (min_mean = 0.125, max_mean = 30) and dispersion (min_disp = 0.5) across cassette batches, retaining genes variable in more than two batches. An scVI model was trained with cassette name as batch covariate, GFS as categorical covariate and RIN_NBB as continuous covariate. The scVI latent space was used for k-nearest neighbor graph construction, UMAP visualization (min_dist = 0.5) and Leiden clustering (resolutions 0.4, 0.6 and 1.0). Identified populations included oligodendrocyte precursor cells (11,794), vascular cells (45,816), astrocytes (53,569), excitatory neurons (64,010), inhibitory neurons (23,401), microglia (25,926) and oligodendrocytes (199,540).

Xenium differential expression analysis

For differential expression analysis relative to pathology in the Xenium dataset, we quantified each cell’s spatial proximity to pathological features. Using pathology masks defined during image analysis, the distance from each cell to the nearest edge of plaque or pTau pathology was calculated as the average Euclidean distance and then log1p-transformed and included as continuous covariates in gene-level differential expression analysis using edgeR. The A axis differential expression was calculated by selecting cells >25 µm from pTau pathology and testing against Aβ distance. For the T axis, cells <25 µm from Aβ were selected, and differential expression was performed against pTau distance.

scANVI label transfer

To further annotate our single-cell data, we performed label transfer from the SEA-AD study³⁴ using scANVI (scvi-tools v.1.1.6)⁵³. For each cell type separately, both datasets were preprocessed and embedded with scVI using raw counts and donor ID as batch key. The scVI model was trained for up to 200 epochs and then used to initialize scANVI with the SEA-AD ‘Supertype’ as label key. scANVI was trained for an additional 200 epochs to predict subtype labels for the octogenarian cells. To assess overlap between scANVI-predicted subtypes and Leiden clusters, we performed a χ² test of independence on a cross-tabulation of cluster and label assignments. Standardized residuals were calculated by normalizing observed versus expected values and dividing by the square root of expected values.

MAGMA analysis of GWAS results

We employed MAGMA⁵¹ v.1.10 to aggregate the GWAS summary data from Kunkle et al.²⁴ to individual genes, selecting genes with an aggregate P value better than 5 × 10⁻⁸.

Statistical analysis

All statistical analyses were performed in Python (SciPy v.1.13.1 and statsmodels v.0.14.5).

Pairwise group comparisons

To compare the average pathology overlapping with tissue domains per sample between phenotypic groups (Fig. 1d), pairwise, two-sided Mann−Whitney U-tests were used.

Enrichment within TDs

Enrichment of a group or statistic within each TD cluster (Figs. 2b,c, 3e,f,i,j and 5e,f and Supplementary Fig. 2) was assessed comparing the distribution of values per hexbin inside a given TD against the pooled distribution from all other TDs using a one-sided Mann−Whitney U-test (alternative = greater, testing for enrichment in a TD).

Delta hexbin plots

Figures 1d, 3g,k and 5g use the same comparison but apply a two-sided Mann−Whitney U-test.

Pathology differences across cortical layers

To examine pathology distribution across cortical layers (Fig. 1c), Aβ and pTau levels were expressed as log-transformed mean percentages per layer and donor. Under the null hypothesis of no group effect, phenotype labels were permuted within each layer while preserving donor-level distributions (5,000 permutations). For each permutation, a between-group test statistic was recomputed and compared with the observed statistic to obtain empirical P values, Bonferroni corrected across layers (α = 0.05).

Cell type enrichment within TD clusters

Cellular enrichment within TDs (Fig. 2d,e) was assessed using cell2location abundance estimates. Expected values under independence were calculated from group means (per TD and per cell type). Enrichment was quantified as log₂(observed / expected). One-sided P values for overenrichment were derived from standardized residuals and Bonferroni corrected for multiple testing.

Multiple testing correction and significance thresholds

For all multiple pairwise or multicluster comparisons, Bonferroni correction was applied to control the familywise error rate at 5%. Unless otherwise noted, tests were two-sided with corrected P < 0.05 considered significant.

Glial Gene Ontology enrichment

Gene Ontology Biological Process enrichment was performed on glial subtypes using GSEA. Ranked gene lists were constructed from pseudobulk differential expression between subtypes and analyzed with gseGO in clusterProfiler (R, v.4.6.0). Gene sets of 10–300 genes were included, and Gene Ontology terms with Bonferroni-adjusted P < 0.05 were considered significant.

To compare enrichment across subtypes, NESs and adjusted P values were merged into a matrix. Dot plots were generated with ggplot2, encoding NES by color and significance by point size. A subset of biologically informative, subtype-distinguishing Gene Ontology terms was selected for final visualizations.

Reporting summary

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