Single cell spatial transcriptomics reveals distinct patterns of dysregulation in non-neuronal and neuronal cells induced by the Trem2R47H Alzheimer’s risk gene mutation

INTRODUCTION The R47H missense mutation of the TREM2 gene is a strong risk factor for development of Alzheimer’s Disease. We investigate cell-type-specific spatial transcriptomic changes induced by the Trem2R47H mutation to determine the impacts of this mutation on transcriptional dysregulation. METHODS We profiled 15 mouse brain sections consisting of wild-type, Trem2R47H, 5xFAD and Trem2R47H; 5xFAD genotypes using MERFISH spatial transcriptomics. Single-cell spatial transcriptomics and neuropathology data were analyzed using our custom pipeline to identify plaque and Trem2R47H induced transcriptomic dysregulation. RESULTS The Trem2R47H mutation induced consistent upregulation of Bdnf and Ntrk2 across many cortical excitatory neuron types, independent of amyloid pathology. Spatial investigation of genotype enriched subclusters identified spatially localized neuronal subpopulations reduced in 5xFAD and Trem2R47H; 5xFAD mice. CONCLUSION Spatial transcriptomics analysis identifies glial and neuronal transcriptomic alterations induced independently by 5xFAD and Trem2R47H mutations, impacting inflammatory responses in microglia and astrocytes, and activity and BDNF signaling in neurons.


BACKGROUND
Genome-wide association studies (GWAS) have identi ed multiple genetic variants associated with Alzheimer's disease (AD) [1].One key discovery identi ed in the TREM2 (Triggering receptor expressed on myeloid cell 2) gene was the R47H missense variant, a strong risk factor for development of Late-Onset Alzheimer's Disease (LOAD) [2,3].TREM2 is an immunomodulatory cell surface receptor expressed primarily in microglia in the brain [4,5], and is activated by a variety of ligands including amyloid-beta (Aβ), APOE, and phospholipids [6].The R47 residue of TREM2 is located within a poly-basic region of the extracellular Ig-like domain, and may modify interactions of TREM2 with its associated ligands [7,8] In AD, microglia exhibit an in ammatory response to Aβ plaques both in human AD brains and in animal disease models [9,10].Evidence increasingly implicates regulation of microglia activation in several ADrelated processes including plaque formation and growth [11], plaque compaction [11,12], protection against dystrophic neurites [13], regulation of development and spread of Tau pathology [14], destruction of perineuronal nets [15,16], and synaptic and neuronal loss [15,[17][18][19][20], though the role of microglia in suppressing or aggravating impacts of AD is currently unclear, and may vary with disease progression [21].
Recent efforts have produced a mouse model of the Trem2 R47H mutation in which bulk RNA-seq analysis has identi ed a unique Trem2 R47H induced interferon signature believed to be associated with microglia in response to Aβ pathology [22].Activation of microglia signi cantly impacts neuronal function and can be neurotoxic [23], but resolving the cell type speci c impacts of this mutation on neuronal populations requires single-cell analysis.Additionally, proximity to Aβ plaques directly impacts glial activation and neuronal transcriptomes [24], requiring spatial transcriptomics to analyze the combined in uence of these effects.
To analyze the impacts of the Trem2 R47H mutation in the context of plaque pathology, we utilized the 5xFAD mouse model, noted for exhibiting strong Aβ pathology at relatively early ages [25].The 5XFAD mouse model has been comprehensively evaluated for preclinical testing applications [26,27].It has been shown that different brain regions (i.e.cortex and hippocampus) in the 5xFAD model have both common and unique gene expression responses to the pathology, and that these changes recapitulate the human AD brain with increased age [27].Hemizygous 5xFAD/homozygous Trem2 R47H (Trem2 R47H ; 5xFAD) mice enable transcriptomic analysis of concerted Trem2 R47H and 5xFAD induced patterns of transcriptomic dysregulation across the brain.
In this study, we utilized MERFISH (Multiplexed Error-Robust Fluorescence In Situ Hybridization) on wildtype (WT), 5XFAD, Trem2 R47H , and Trem2 R47H ; 5xFAD mice to probe spatial gene expression in single cells [28,29].Previous studies using MERFISH have analyzed spatial transcriptomics of neurodegeneration in aging [30], and in microglia activation [31], proving the e cacy of MERFISH in investigating transcriptomic dysregulation in the brain.However, studies analyzing spatial impacts of Alzheimer's disease in mouse models have been limited either in cell type speci c impact analysis [32], spatial resolution [33], or the size and number of imaged regions [24].In contrast, we present here a single-cell resolution, spatial transcriptomic atlas of 5xFAD and Trem2 R47H induced alterations across whole mouse coronal sections.Our ndings reveal spatially localized cell type speci c plaque and Trem2 R47H induced transcriptome dysregulations in both glia and neurons in multiple cortical and subcortical brain regions.

Animals
The Trem2 R47H mice used in this study were derived from the same Trem2 R47H NSS (Normal Splice Site) mouse colony with C57/B6 background as previously reported (Jax stock #034036 [22]).All animals were bred and raised by the Transgenic Mouse Facility at UCI under a regular light/dark (12h/12h) cycle with ad libitum access to food and water.All animal care and related experimental procedures were conducted following the highest ethical standards and were approved by the UC Irvine Institutional Animal Care and Use Committee.
Preparation of mouse brain sections and MERFISH Mice (wild-type, Trem2 R47H , 5xFAD, and 5xFAD; Trem2 R47H ) were euthanized at 12 months of age via carbon dioxide inhalation followed by transcardiac perfusion with chilled phosphate-buffered solution (PBS, pH7.2).Brains were quickly collected with hemispheres bisected along the midline and separately embedded with another hemisphere from the designated genotype pair (WT with 5xFAD and Trem2 R47H with 5xFAD;Trem2 R47H ) in a square tissue mode with Tissue-Tek® OCT mounting medium.
Each pair of brains was ash frozen in dry ice-chilled isopentane and stored at -80°C until cutting.
To prepare cryosections for MERFISH, two hemisphere OCT blocks containing 4 samples were combined and sectioned at -20°C on a Leica CM1850 cryostat.A 10-mm-thick coronal slice containing both hippocampus and subiculum regions was collected onto a specially coated 4cm-diameter coverslip (Vizgen Item# 10500001), xed in 4% paraformaldehyde in PBS in a 6cm petri dish, and stored in 70% ethanol at 4°C until MERFISH probe hybridization after a brief rinse with PBS.
MERFISH was performed according to Vizgen's protocol.Brie y, merslides with mouse brain sections were rinsed with Vizgen Sample Preparation Buffer (SPB) after the removal of 70% ethanol, incubated with Vizgen's Formamide Wash Buffer (FWB, 30 min at 37°C), hybridized with a customized mouse gene panel containing speci c binary-coded probes for selected 300 mouse genes (VZG171, Supp Table 1) in a para lm-sealed plate (~ 40h at 37°C), and washed with FWB twice (30 min at 47°C).After the aspiration of FWB, the brain sections were be embedded with Vizgen gel mix, incubated in the clearing mix (5mL with 50 mL protease K, overnight at 37°C), stained with Vizgen DAPI/poly(T) mix from the Vizgen 300gene imaging kit (10 min at RT) after rinse with SPB, and washed with FWB (10 min at RT). Subsequently, the merslide was thoroughly rinsed with SPB before being carefully assembled into a gasket chamber.
Once assembled, the merslide was then uploaded into the MERSCOPE for imaging.The MERFISH imaging was done on the MERSCOPE with a Vizgen 300-gene imaging kit after the addition of the imaging buffer activator and RNase inhibitor (100 mL) as stated in the manual.The imaging process was controlled by the MERSCOPE program.Once MERFISH imaging process was completed, the output les were transferred for in-depth analysis with MERFISH Visualizer and our customized bioinformatic pipeline.
Imaging occurred in 5 total batches: 1) one 5xFAD animal (F) and one WT animal (F), two brain slices each; 2) one Trem2 R47H ; 5xFAD animal (F) and one Trem2 R47H animal (M), two brain slices each; 3-5) one slice from each genotype, no duplicated animals, all male.One WT sample failed imaging QC in the nal batch.Samples from the same animal are aggregated for subsequent analysis.

Processing of MERFISH datasets
Following automated transcript decoding and error correction via the MERSCOPE software pipeline, individual cells were segmented using the machine learning model cellpose [34], which was custom trained on DAPI stained slices collected previously, captured at the same resolution.Next, transcripts were assigned to individual cells.Cells exhibiting volume more than 1800 µm 3 , or less than 50 transcripts per cell were removed.
Data was then processed using our Scanpy [35] based custom pipeline, namely, library size normalization, log transformation, regression of sequencing depth per cell as a confounding variable, standard scaling, PCA transformation, batch integration using harmony [36], and UMAP dimensionality reduction.Next, marker genes were computed (scanpy's sc.pl.rank_genes_groups function) and matched to cell type reference atlases including the Allen institute cortex and hippocampus dataset [37], and the mousebrain.orgsingle-cell reference atlas [38].Cell type annotations were re ned and corrected using spatial coordinates, particularly ensuring cortical-layer-speci c neuron types were correctly organized.As part of the quality control process, cell types were subclustered, and cell subclusters exhibiting markers for other (typically glial) cell types were excluded from analysis as contaminated.
Region annotation of major spatial domains was performed using a semisupervised approach based on superposition and manual annotation of cell types on the appropriate Allen mouse brain coronal atlas slice [39], guided by the spatially localized cell types for ne region selection, particularly in hippocampus.
Subclustering was performed by subsetting to the desired cell type and running the same pipeline on the individual cell types.Subclusters containing less than 5% of the total cells were excluded from analysis.Markers for identi ed subtypes were identi ed using sc.pl.rank_genes_groups.Genotype bias was computed by rst normalizing the number of cells in each subcluster by the total number of cells contributed from that genotype, and then normalizing by genotype, to obtain cell type proportions in each subcluster.The following thresholds were used to identify genotype bias: if a single genotype proportion for a given subtype exceeded 33% (25% being uniform distribution), the subtype was considered upregulated in that genotype.Additionally, if the combined proportion of two genotypes exceeded 60% (50% being uniform distribution), the subtype was considered upregulated in that pair of genotypes.This latter was restricted so to identify only Trem2 R47H , and 5xFAD speci c upregulation (e.g.any subtypes coupregulated in WT and Trem2 R47H ; 5xFAD were not considered for further analysis).

Pseudobulk Differential Expression Analysis
Cells were divided by cluster, genotype and batch.Genes present in fewer than 15% of cells were not analyzed for differential expression, due to limited accuracy of differential expression in low frequency genes [40].Data for each cell type was aggregated by genotype and batch to construct pseudobulk replicates [41].Samples with fewer than 50 cells of a given cell type were removed.
Due to unbalanced genotype proportions in individual batches, pairwise differential expression was performed separately for each genotype pair studied.Subsetting pseudobulk replicates to those associated with the compared genotypes, we utilized a linear mixed effects model (lme4 and multcomp packages [42,43]), with batch as the random effect.Gene signi cance was identi ed using an absolute log fold change of 0.35, and an adjusted p -value of 0.05.The Benjamini-Hochberg method [44] was used for p-value correction.Fold changes were computed based on the inferred values by the LME model.
Log fold changes were computed using a base of 2, unless otherwise stated.

Continuous plaque proximity differential expression
Continuous plaque proximal differential expression was computed at the single cell level for all cell types.Omitting the conversion to pseudobulk differential expression DESeq2 [45] was utilized to analyze differential expression contiguous to plaques, by including distance to plaque as a continuous covariate, and computing differential expression as a function of distance to plaque.Genotype was not used as factor in the model due to low sample numbers, though only 5xFAD and Trem2 R47H ; 5xFAD samples were utilized in this analysis.Only adjusted p-value (< 0.05) was used for identi cation of differentially expressed genes in this analysis.

Differential expression between regions
To perform differential expression of glial cell types between different spatial regions in the same genotype, we utilized a pseudobulk approach.After subsetting to individual cell types and removing genes expressed in < 15% of cells, we excluded sample-region combinations in which fewer than 50 cells are identi ed.We note that this completely excludes OPC analysis in the dentate gyrus due to lack of cells.At this point, we created pseudobulk samples for each sample-region combination (typically 10 pseudobulk samples for each slice, based on the 10 annotated regions).We then computed differential expression within the same genotype, comparing expression in each region to the average expression in all other regions, to identify spatially variable genes.We utilized a linear mixed effects model, with batch as the random effect, and gene signi cance was identi ed using an absolute log fold change of 0.35, and a p -value of 0.05.

Cellpose segmentation of plaques
Initial visualization of DAPI staining in merslides appeared to show plaque staining in addition to cell nuclei.We con rmed this by costaining additional prepared 5xFAD merslides with DAPI and thio avin S, to con rm that DAPI stains Aβ plaques in addition to cell nuclei.To segment Aβ plaques in the DAPI image, we utilize the cellpose GUI to identify plaques based on DAPI brightness, plaque size (frequently larger than cells), and the presence of brils and irregular cellular shape.We con rmed the accuracy of the model predictions based on false positive rates and F1 scores.We also demonstrate that this model does not segment cells, by showing that transcript density is signi cantly lower in plaques than in cells, and by checking individual elds of view for plaque segmentation, comparing annotated images to model predictions.A total of 100 images were used for training, and 25 for validation.

Filtering of differentially expressed genes
Due to irregular cellular shapes, cellular processes distal to the soma, microglia phagocytosis, and possible segmentation errors, overall upregulation in expression within one cell type, may overlap into another.This pattern was particularly noticeable in glia, which frequently exhibited expression overlap in marker genes, as well as neuronal speci c genes (such as Slc17a7).As an example, Gfap was ubiquitously expressed, and differentially expressed in 5xFAD samples compared with WT in many neuronal cell types.Yet its expression is limited almost exclusively to astrocytes in previous mouse brain cell atlases [37,38].
To lter out these false positive differentially expressed genes, we heavily annotated all major glia subtypes for gene expression of all 300 genes in the panel, using both the mousebrain.org[38] and Allen institute datasets [37] (Supp Table 2).Initial annotation utilized the Allen institute dataset, requiring a trimmed mean expression greater than zero.This was supplemented using the mousebrain.orgdataset, which is not restricted to the hippocampus and contains activated microglia and astrocyte cell types.Differential expression of glial cells was then subset to genes known to be expressed in these cell types as annotated.Differential expression in neurons was ltered for expression of glial cell type markers, and disease associated microglia and astrocyte markers.While we recognize this process may introduce bias into the differential expression results, raw results exhibited signi cant cell type induced biases complicating analysis of results.Raw differential expression results for glia are shown in supplemental gures, and raw differential expression results for all cell types are available as supplemental tables.
Additionally, hippocampal and inhibitory differential expression is impacted by hippocampus size, which varies across samples.Subclustering identi ed spatially localized subgroups for each cell type, primarily in the ventral hippocampus.Marker genes for these clusters were identi ed, and differentially expressed genes overlapping these subsets were removed from the analysis.

Computation of cell density within major regions
To compute the area of identi ed regions in individual slices, we utilize the python package alphashape [46], a method for automatically constructing concave bounding envelopes of point clouds.For each region, we subset to all cells contained in that region, and utilize alphashape with an alpha of 0.015, and compute the area of the resulting polygon.
For computing density of glia in regions proximal and distal to plaques, all cells in each individual slice within 100 µm (proximal) and between 100 and 500 µm (distal) were used to compute the combined area of all regions distal and proximal to plaques.This was then used as the normalizing factor to obtain cell densities proximal and distal to plaques.When investigating density of individual cell types, we used the region identi ed using all cells within 50 µm of a single cell of the given subtype, as some cell types were insu ciently dense for area inference via alpha shape.

Statistical methods
Excluding differential expression (described above), statistical tests are described in the text.We utilized two-sided tests unless otherwise stated.

Impact of 5xFAD and Trem2 R47H mutations across major brain cell types
In this study we investigate regional, plaque proximal, and genotype speci c gene expression changes induced by the Trem2 R47H mutation.As this mutation is not su cient to induce amyloid plaque pathology in mice, we utilize a hemizygous 5xFAD/homozygous Trem2 R47H mouse model which induces Aβ pathology in concert with the Trem2 R47H mutation, compared with matched 5xFAD (Aβ pathology only), Trem2 R47H , and WT controls (Fig. 1A).By comparing these four genotypes, we uncover the transcriptomic alterations induced speci cally by 5xFAD transgenes (independent of Trem2 R47H mutation), speci cally by Trem2 R47H (independent of 5xFAD), and those induced by a combination of 5xFAD and Trem2 R47H mutations.
We performed spatial transcriptomic analysis using MERFISH on 19 coronal half sections from 15 total animals with WT, 5xFAD, Trem2 R47H and Trem2 R47H ; 5xFAD genotypes at 12 months of age.After quality control, this dataset results in 432,794 cells.Using a 300 gene panel, we identi ed 37 major cell types, and transcriptomically and spatially mapped 5xFAD and Trem2 R47H transcriptomic alterations at the single-cell level.We also identi ed Aβ plaque locations in the same samples and assessed their relationship to spatial gene expression (Fig. 1B).
After spatial transcript decoding (Fig. 1C), cells were processed using our single-cell pipeline (Supp Fig. 1A-B), and clusters were identi ed based on reference to known cell type markers (Supp Fig. 1C), in conjunction with spatial location (Fig. 1D).Color coding genotype information on the UMAP shows strong 5xFAD induced cell type composition changes, particularly in astrocytes and microglia (Fig. 1E).Hierarchical clustering identi ed initial splits between non-neuronal and neuronal cells, followed by excitatory vs inhibitory, and spatial (subcortical, hippocampal, cortical) based splits in excitatory cell types (Fig. 1F).
Visualization of neuronal cell types (Fig. 2A, Supp Fig. 2) shows strong spatial localization, commensurate with previous region-based studies and atlases.Hippocampal excitatory cells de ne the primary structures of the hippocampal formation (DG, CA1, CA3), while cortical excitatory neurons divide into distinct layers across most of the cortex.We identify and visualize cell type speci c markers for these distinct neuron types (Fig. 2B) to verify spatial delity with raw decoded transcripts.
Next, we segmented major brain regions, subdividing the cortex into three subregions: the neocortex (somatosensory, visual, parietal, retrosplenial, and auditory cortices), the limbic cortex (perirhinal, ectorhinal, entorhinal, and piriform cortices), and the cortical amygdala, and identify major structures in hippocampal and subcortical regions.This resulted in 10 identi ed major brain regions (Fig. 2C).
Overall, MERFISH spatial transcriptomics enables detection of high-level cell type clusters, visually identi able and quanti able transcriptomic differences in microglia and regional annotation and assignment of individual cells to speci c coarse-grained spatial regions.
Glial and neuronal transcriptomes are affected by nearby plaques Spatial transcriptomics can reveal local effects of pathology, such as Aβ plaques, on the regulation of gene expression in nearby cells.By co-staining coronal brain slices with both DAPI and thio avin S (a canonical stain for Aβ plaques), we observed that DAPI brightly labels Aβ plaques in addition to nuclei [48] (Fig. 3A).We therefore applied DAPI staining to MERFISH prepared coronal slices and a machine learning approach to automatically detect and segment plaques in each of the MERFISH samples.
DAPI stained plaques are visually distinguishable from nuclei by their large size, greater brightness, and brous morphology and lack of circular cell soma shape (Fig. 3B-C).These features enable manual annotation of plaques in individual elds of view.We trained a modi ed cellpose model [34] to detect plaques, but not cells (Fig. 3B).In testing the model on holding out annotated plaque data, we identi ed only 1 false positive across 25 FOVs, and 28 total plaques, with an F1 score of 0.89.
We analyzed each 5xFAD and Trem2 R47H ; 5xFAD section using this model (Supp Fig. 3A-B) and veri ed that 1) the model does not detect cells (Fig. 3B), 2) the predicted plaques are morphologically distinct from cells (Fig. 3C), and 3) the predicted plaques have signi cantly lower transcript density when compared to cells (Fig. 3D, Wilcoxon rank-sum test, p < .01).In total we identi ed 5,616 plaques distributed across multiple brain regions.
Next, we assessed whether plaques appear proximal to neurons.The distance from a plaque to the closest neuron was signi cantly larger than the distance from a neuron to its closest neuronal neighbor (5xFAD: minimal plaque to neuron distance 56.4 +/-10.7 µm, minimal neuron to neuron distance 21.6 +/-1.24µm, p = 0.012; Trem2 R47H ; 5xFAD: minimal plaque to neuron distance 48.7 +/-2.53 µm, minimal neuron to neuron distance 21.0 +/-0.733µm, p = 1.47 x 10 − 4 , mean +/-s.e., plaques in corpus callosum excluded from analysis due to lack of nearby neurons, t-test).We examined the typical distance of each neuronal cell type to the nearest plaque.This analysis showed that subiculum excitatory neurons, layer 5 and layer 6 excitatory neurons have the lowest median distance to plaques among identi ed cell types (Supp Fig. 3C).However, none of the top 5 neuron types (ranked by median distance to plaque, excluding subiculum excitatory and SST-Chodl cells due to low cell numbers), exhibited signi cant density variation between plaque proximal (< 100 µm) and distal (100-500 µm) regions.This implies that neuronal plaque proximity is driven primarily by plaque density in the associated regions.Additionally, plaques on average form in regions nearly twice as far from the nearest neuron as the typical distance between neurons, but the average neuronal density does not appear to be decreased in plaque proximal vs plaque distal regions, implying a variation in plaque to neuron distance at the microscale (< 100 µm), but not at larger scales (< 500 µm).
We next tested whether cells within 100 µm of the nearest Aβ plaque have altered patterns of gene expression (Fig. 3H).Due to the relatively low cell abundance proximal to plaques, we aggregated 5xFAD and Trem2 R47H ; 5xFAD samples, and separated individual cells by cluster.We analyzed plaque proximity based differential expression with two techniques.First, we identi ed cells within 100 µm of a plaque center.Using DESeq2 and treating cells as independent samples, we identi ed genes whose expression correlated with proximity to the nearest plaque (Fig. 3H, Supp Table 3).Continuous effects were identi ed primarily in microglia and astrocytes, with microglia showing an upregulation of typical DAM associated genes (e.g.Csf1, Apoe, Cst7), and a downregulation of P2ry12, a homeostatic microglia associated gene [50].Similarly, C4b, Clu, and Gfap, markers of a previously known disease associated astrocyte (DAA) phenotype were also upregulated near plaques [51].
To validate these ndings and to account for variability across biological replicates, we additionally performed a pseudobulk analysis of differential expression between plaque-proximal (within 100 µm of the closest plaque) and plaque-distal (100-500 µm to closest plaque).We applied a linear mixed effects model to pseudobulk expression for each cell type in each sample, accounting for batch as a random effect (Fig. 3I, Supp Table 4).Additionally, we ltered genes based on their known expression in each cell type from previous single-cell atlases [37,38], to avoid spurious identi cation of differentially expressed genes due to technical (errors in segmentation) or biological (phagocytosis, overlapping cellular processes) effects [52].
Pseudobulk analysis was generally consistent with the DESeq2 results and identi ed both glial and neuronal changes (Fig. 3I, top row).Microglia and astrocytes exhibited typical disease associated pro les in cells located proximal (< 100 µm) to plaque centers.However, Nnat expression in astrocytes and oligodendrocytes, and Mmp14 expression in astrocytes decreased near plaques.This result contrasts with previous ndings in humans and other mouse models showing upregulation of Mmp14 in reactive astrocytes in Alzheimer's disease (AD) [53].
The pseudobulk analysis also revealed notable changes in gene expression affecting neurons proximal to plaques (Fig. 3I, bottom row).L6b neurons showed lower Ngf expression near plaques, a gene therapy target in AD [54].Nr2f2, upregulated near plaques, is known to be dysregulated by AD associated single nucleotide polymorphisms in the APOE enhancer [55].L2 intratelencephalic (IT) neurons near plaques showed downregulation of Dkk3 (a WNT signaling modulator whose presence reduces Aβ pathology in mouse models [56]) and of the potassium ion channel subunit Kcnd2 [57] near plaques.L5 NP cells show Grm1 upregulation and Chrna7 downregulation near plaques.Parvalbumin-expressing inhibitory cells shows Grin2a, Zbtb20, and Plagl1 downregulation near plaques.Excitatory neurons in the cortical amygdala exhibited downregulation of Ntf3 (associated with nervous system maintenance[58]), Nptx1 (associated with synapse remodeling, but typically downregulated in previous studies of cortical neurons near plaques [59]), and Camk2g (implicated in synaptic plasticity[60]).Because there were few plaques in subcortical regions, we did not test plaque-associated differential expression for neuronal cell types in this region.
We identi ed 19 differentially expressed genes in microglia and 8 in astrocytes across all 4 pairwise comparisons.By contrast, we found 1-2 differentially expressed genes in oligodendrocyte (OGC) and oligodendrocyte precursor cells (OPC) cell populations (Fig. 4A), and none in the other non-neuronal cell types.Microglia and astrocytes primarily exhibited 5xFAD dependent changes (similar differential expression results for both 5xFAD vs WT, and Trem2 R47H ; 5xFAD vs Trem2 R47H ), replicating the DAM/DAA gene upregulation and homeostatic gene downregulation identi ed in the plaque proximity analysis.
Interestingly, neither Itgax nor Cd74 were identi ed as differentially expressed in plaque proximity analysis of microglia, whereas they are upregulated 9.58 and 15.7-fold in 5xFAD compared with WT animals, and 19.7 and 26.0-fold in Trem2 R47H ; 5xFAD compared with Trem2 R47H animals.The Trem2 gene itself showed a small reduction in expression dependent on the Trem2 R47H mutation, contrary to previous studies [22].This may be due to reduced binding of gene probes overlapping the mutated region.
The consistent variation in both Trem2 R47H vs WT and Trem2 R47H ; 5xFAD vs 5xFAD comparisons indicates this may be a plaque independent effect and corroborates the overall lower plaque burden in Trem2 R47H ; 5xFAD samples.
To explore the effects of AD risk genes in speci c glial subtypes, we subclustered the microglia and astrocyte subpopulations.We identi ed several small clusters of microglia that appear to express neuronal or other glial markers, and we con rmed that these cells are located near cells expressing these markers.We removed these cells from this portion of the analysis.After removal, subclustering identi es 7 microglia clusters (Fig. 4B1).Pseudotime analysis identi ed a single linear trajectory across all microglial cell types (Fig. 4B2).We next examined the genotype proportions of these clusters.After normalizing by the number of cells per sample, we averaged across samples of the same genotype, and computed cluster proportions.This identi es a clear 5xFAD dependent bias, with two clusters (labeled homeostatic) exhibiting > 80% proportion coming from non 5xFAD (i.e.WT and Trem2 R47H ) samples.The remaining ve clusters corresponded to disease associated microglia (DAM) enriched in 5xFAD and Trem2 R47H ; 5xFAD mice (Fig. 4B3).
We next aggregated homeostatic and DAM subgroupings and identi ed regional spatial biases (Fig. 4B4).DAMs were enriched in hippocampal area CA1.They were also enriched in thalamus, and midbrain, despite the relative lack of plaque density in these regions compared to the CA1 and CC (Fig. 3E).Finally, we identi ed markers for the individual microglia subpopulations, and plot normalized expression (Fig. 4B5).
We focused on the analysis of the genes differentially associated with late-stage DAMs as several genes exclusive to late-stage DAM (DAM2) were included (Itgax, Cst7, Csf1, Ccl6), as well as genes present across both stages (Apoe).Except for Ccl6, all of these genes are differentially expressed in 5xFAD and Trem2 R47H ; 5xFAD, with primary expression of DAM2 genes in C3-5 (later pseudotime).Apoe is evenly distributed across C2-5, re ecting its overexpression across the DAM developmental timeline (Fig. 4B5) [61].
Subclustering the astrocyte subpopulations, we aggregated clusters not exhibiting genotype speci c bias (see methods for thresholds) into a single cluster (C1), retaining the genotype biased clusters (Fig. 4C1).
Pseudotime trajectory analysis (Fig. 4C2) did not yield a distinctive pattern, however after analysis of genotype bias (identifying C1 as unbiased, C2/C3 as DAA, and C4/C5 as upregulated in WT/Trem2 R47H samples, Fig. 4C3), we note that C5 and C4 exhibited distinct spatial distributions, with C4 appearing exclusively in cortex and hippocampus, and C5 appearing in subcortical regions (Fig. 4C4).The DAA exhibited a similar regional speci city, with C2 restricted to cortex and hippocampus.Cluster markers are identi ed and plotted (Fig. 4C5).
We next examined the spatial distribution of DAM and DAA cells by region.Disease associated microglia were enriched in the CC, subiculum and subcortical regions (Fig. 4D1).Computing the proportion of microglia identi ed as DAM by region (Fig. 4D2) showed similar proportions of DAMs between 5xFAD and Trem2 R47H ; 5xFAD samples by region, except in the DG, thalamus, midbrain, and hypothalamus.This corresponds with the plaque density bias in 5xFAD samples (Fig. 3G).In the cortex, we saw a signi cant increase in DAMs in the lower cortical layers (L5/L6) compared with the upper cortical layers (L2/L3) (p < .0001,linear mixed effects model, Fig. 4D3).
Disease associated astrocytes were concentrated in the CC and surrounding areas (Fig. 4E1).The only regions exceeding 40% DAA proportion are the CC and CA1 (Fig. 4E2).Virtually no disease associated astrocytes are present in upper cortical layers, but this population was signi cantly upregulated in deeper cortical layers (Fig. 4E3).We did not nd signi cant genotype speci c effects in other glial cells.
To further analyze the spatial variation of gene expression, we performed direct pseudobulk differential expression analysis of microglia, astrocytes, oligodendrocytes and oligo-precursors across the 10 identi ed major brain regions (Supp Fig. 4, Supp Table 6).We compared each region with the average across the remaining 9 regions.We also computed regional cell density and cell proportion for each of these cell types.
Analysis of microglia (Supp Fig. 4A) shows that ~ 60% of spatially variable genes are also differentially expressed across genotypes.For example, the canonical late-stage DAM markers Cst7 and Itgax were signi cantly upregulated in the corpus callosum.A small number of other genes (Ctss, C1qa, Zbtb20, Ly9, Tmem119) had spatially variable patterns of expression that were consistent across WT and Trem2 R47H mice and dysregulated in 5xFAD and Trem2 R47H ; 5xFAD mice.Microglia cell populations also show drastic increases in both cell proportion and density across all brain regions.Astrocytes (Supp Fig. 4B) exhibited large numbers of spatially variable genes with consistent patterns of expression across all genotypes (e.g.Erbb4, Nnat, Grin3a, Mmp14, Id4, Pax6, etc).We also found spatial variation in several disease associated genes (Aqp4, Gfap).These spatial variations were primarily observed between cortical and subcortical (thalamus, midbrain, hypothalamus) regions.However, astrocytes exhibited little genotype speci c cell proportion or density variations between regions.
Oligodendrocytes exhibited two separate gene groupings (Supp Fig. 4C).One set (Snca, Dlg4, Nnat, Robo1, S100b, Ptgds) showed spatially variable expression across multiple regions, particularly between cortical/hippocampal and subcortical regions.The other set of genes (Adam10, Psen1, Olig1, etc) is primarily upregulated in CC and downregulated in amygdala, with very little variation in other regions.This pattern is not 5xFAD or Trem2 R47H dependent and was observed even in WT oligodendrocytes cells.This second pattern is not replicated in oligodendrocyte precursor cells, though a signi cant cortex vs subcortical divide is present in OPCs (Supp Fig. 4D).Neither cell type exhibits signi cant genotype dependent cell proportion changes within regions.
Overall, our data show that spatial variation in microglia and astrocyte gene expression is more affected by 5xFAD than by Trem2 R47H .Both disease-associated microglia and astrocytes exhibit speci c spatial distributions.DAMs were distributed across the coronal section, but concentrated in the CC and subcortical regions, and DAA were biased almost exclusively to the CC and surrounding regions.Regional transcriptional variations were primarily impacted in 5xFAD for microglia and astrocytes, and both 5xFAD and Trem2 R47H mutations were independent of regional variations in oligodendrocytes and oligodendrocyte precursors.

Neurons exhibit complex transcriptomic impacts of 5xFAD and Trem2 R47H mutations
We next performed differential expression analysis for each of the four comparisons (5xFAD vs WT, Trem2 R47H ; 5xFAD vs Trem2 R47H , Trem2 R47H vs WT, and Trem2 R47H ; 5xFAD vs 5xFAD), as well as subclustering analysis, for each of the neuronal cell types (Supp Table 5).Analysis of cortical neurons identi es differentially expressed genes for all these comparisons in each cell type (Fig. 5A-C), as well as genotype biased subclusters for most neuron cell types (Fig. 5D-E, Supp Fig. 5).We rst considered genes consistently identi ed as differentially expressed across multiple cortical neuronal cell types.
Cdh12, associated with calcium ion binding [62], was differentially expressed in 5 out of 9 cell types.Both L3 IT and L5 IT neurons exhibit upregulation of Cdh12 in Trem2 R47H ; 5xFAD over 5xFAD genotypes.Ntrk2, which encodes TrkB, a high a nity receptor for BDNF [63], is upregulated in Trem2 R47H ; 5xFAD vs 5xFAD and Trem2 R47H vs WT comparisons in L2 IT, L3 IT, L6 IT and L6 CT excitatory cell types.On the other hand, Bdnf itself, expected to decrease in the 5xFAD context, was identi ed as signi cantly decreased only in L2 IT and L6b neurons.Fos, a molecular marker of neuron activity [64], was consistently identi ed as differentially expressed across 6 of the 9 cell types, for at least one comparison.In each case this gene was downregulated, implying downregulation of Fos induced by both 5xFAD and Trem2 R47H mutations.
All other genes were differentially expressed in at most 3 cortical excitatory cell types.Of these, the most interesting are Wfs1, recently implicated in Tau clearance in AD [65], which is downregulated in Trem2 R47H ; 5xFAD compared with 5xFAD, and Grm2, a glutamate receptor downregulated in Trem2 R47H ; 5xFAD compared with 5xFAD in L3 IT and L6 IT neurons.
Subclustering cortical neurons identi ed genotype speci c subpopulations in 7 cortical excitatory cell types (Fig. 5C,D; Supp Fig. 5).In all but one cell type (L3 IT), this represents a WT (or WT/Trem2 R47H ) enriched, and thus 5xFAD/Trem2 R47H ; 5xFAD reduced subpopulation.For IT cell populations, these genotype speci c subtypes did not exhibit signi cant spatial localization.However, 5xFAD/Trem2 R47H ; 5xFAD reduced subtypes in L2 IT are spatially localized to the retrosplenial and visual cortices, and 5xFAD/Trem2 R47H ; 5xFAD reduced subtypes in L5 NP are spatially localized in the retrosplenial cortex near the subiculum, a plaque dense environment (Fig. 5E, Supp Fig. 5).This latter subtype upregulated Sulf2 (antibody staining has shown this is reduced in AD [66]) and Cplx1 (regulates synaptic transmission by preventing neurotransmitter release prior to action potential [67]) as top differentially expressed genes.
Hippocampal CA1 excitatory neurons (Fig. 6B) exhibited upregulation of the expression of Ntrk2 and Mapk1, associated with the MAPK signaling pathway, in Trem2 R47H ; 5xFAD compared with 5xFAD animals.CA3 excitatory neurons showed several differentially expressed genes (Rph3a, Itga7, Hs3st1), but almost exclusively in the 5xFAD vs WT comparison, though Ntsr1 was upregulated in both Trem2 R47H ; 5xFAD vs 5xFAD, and Trem2 R47H vs WT comparisons.The dentate gyrus showed upregulation of Dkk3 and downregulation of Adgra1 in both 5xFAD dependent comparisons.
Subclustering subcortical and hippocampal neurons (Supp Fig. 6) reveals greater genotype proportion heterogeneity than in cortical excitatory neurons.In contrast with cortical excitatory neurons, hippocampal and thalamic excitatory neurons clustered into large numbers of variable genotype proportion subclusters.Thalamic excitatory neurons subcluster into three genotype enriched sets, including a 5xFAD/Trem2 R47H ; 5xFAD enriched subtype a Trem2 R47H /Trem2 R47H ; 5xFAD enriched subtype and a WT/Trem2 R47H enriched subtype.CA1 excitatory neurons identi ed 6 genotype enriched subclusters, though two of them are spatially localized in the ventral CA1, which is not included in some samples.These include two 5xFAD/Trem2 R47H ; 5xFAD enriched subtypes and two WT/Trem2 R47H enriched subtypes.CA3 excitatory neurons subclustered into 5 genotype biased subclusters, including several localized to the ventral hippocampus (Fig. 6D, ventral hippocampus clusters not shown).One subcluster, labeled C2, present primarily in 5xFAD and Trem2 R47H ; 5xFAD samples, is spatially positioned in the intersection of the CA3 and dentate gyrus.This cluster upregulated Rph3a and Dkk3 (Supp Fig. 6C).
Overall, neuronal populations exhibit transcriptional alterations associated with both 5xFAD and Trem2 R47H mutations.In cortical excitatory neurons, these changes are frequently replicated across cell types.Thalamic excitatory neurons, uniquely among subcortical populations, exhibit signi cant 5xFAD and Trem2 R47H induced transcriptomic alterations.In the hippocampus, the CA1 shows the most transcriptional alteration among genes measured in this study.Neuronal subclusters show both genotype enrichment, and spatial localization, implying regional population variations induced by 5xFAD and Trem2 R47H mutations.

Discussion
Trem2 R47H is strongly associated with the development of Late-Onset Alzheimer's Disease.Here, we investigated the spatial transcriptomic impacts of this critical mutation in the context of the 5xFAD mouse model of amyloidosis, as well as in a WT mouse background.Across 19 coronal slices we pro led over 400,000 cells and examined transcriptome dysregulation in neuronal and glial cell types.This analysis provides a broad perspective, enabling analysis of regional and cell type-speci c transcriptome dysregulation at the single-cell level.We improved on previous spatial transcriptomic analyses of Alzheimer's disease mouse models, which were generally limited either in imaged area or in spatial resolution.
Previous characterization of Trem2 R47H ; 5xFAD brains found an initial effect whereby the presence of the variant impeded the microglial response to plaques and exacerbated surrounding dystrophic neurites.However, this suppression of microglial response to plaques subsided at later disease stages, resulting in expected numbers of plaque-associated microglia, and no signi cant changes in DAM gene expression via bulk-tissue RNA-seq between Trem2 R47H ; 5xFAD and 5xFAD hippocampi [22].This is consistent with our ndings here in regions exhibiting early plaque development (hippocampus, limbic cortex).Despite a seemingly appropriate microglial response to plaques at these later stages, this previous study observed unique emergence of an interferon signature by TREM2 R47H , coupled with increased plasma neuro lament light chain, a marker of neuronal damage [22].In agreement with this, we found that the presence of TREM2 R47H has signi cant impacts on gene expression within neurons.
As indicated, both microglia and astrocytes showed much higher impacts from the 5xFAD transgenes than from the Trem2 R47H mutation, resulting in a larger number of differentially expressed genes, and higher magnitude expression changes in both 5xFAD dependent comparisons than in the Trem2 R47H dependent comparisons.Of those few genes that were differentially expressed in a Trem2 R47H dependent manner, several are associated with a homeostatic microglial state (P2ry12, Tmem119), in Trem2 R47H ; 5xFAD compared with 5xFAD mice.TREM2 is required to transition microglia from a homeostatic to a DAM state, and the lack of downregulation of these homeostatic markers is consistent with a partial loss of TREM2 ability to mediate this transition with the R47H mutation.
Both disease-associated microglia and astrocytes exhibited distinct spatial distributions, with signi cant concentrations of both populations at the corpus callosum (CC), though DAM were distributed more evenly across the brain, while DAA were restricted almost exclusively to the CC.While proportions of disease associated astrocytes were consistent between 5xFAD and Trem2 R47H ; 5xFAD mouse models, 5xFAD proportions of disease associated microglia were higher than Trem2 R47H ; 5xFAD in the dentate gyrus, midbrain, thalamus and hypothalamus, correlating with plaque density.
Analysis of plaque proximal transcriptomic dysregulation identi ed the prominent disease associated microglia and astrocyte signatures, with slight distance dependent variations (e.g.P2ry12 was identi ed as downregulated near plaques for microglia within 100 µm, but not between cells within 100 µm, and those 100-500 µm from the nearest plaque, with the reverse true for Tmem119).Cd74 was not identi ed as upregulated near plaques but was strongly upregulated in 5xFAD and Trem2 R47H ; 5xFAD mice (compared with WT and Trem2 R47H respectively).Upregulation of this gene is thought to precede nal differentiation into DAM states [75], and exhibits homogeneous upregulation induced by 5xFAD independent of plaque proximity.
Examining genotype induced differential expression across cortical excitatory cell types, we identi ed Fos, Wfs1, and Grm2 as downregulated by the Trem2 R47H mutation and Ntrk2, Bdnf, and Cdh12 as upregulated.Downregulation of Fos and Grm2 is indicative of a decrease in activity and synaptic transmission, and is consistent with the LTP impairments observed in these mice at this age [22].Wfs1 is a marker for a unique neuronal population in the entorhinal cortex that modulates spatial memory and is implicated in late stage induced hypoactivity in the hippocampal formation, and is also linked to Tau pathology [65,76].Nr2f2 and Bdnf are both critical to the BDNF signaling pathway, which activates ERK and Akt pathways to maintain neuronal survival and synaptic plasticity [77].
Thalamic excitatory neurons exhibited both 5xFAD and Trem2 R47H impacts, and many differentially expressed genes between WT and 5xFAD mice, including downregulation of Bdnf, have previously been linked to Alzheimer's.CA1 excitatory neurons showed the largest number of differentially expressed genes among non-cortical neurons, among genes in this panel.However, relatively few differentially expressed genes were identi ed between 5xFAD and WT mice.The only consistent differentially expressed genes induced by the Trem2 R47H mutation were Ntrk2 (upregulated in Trem2 R47H and Trem2 R47H ; 5xFAD) and Fos, similar to results in cortical excitatory neurons.Inhibitory neurons exhibited consistent differential expression of Epha10 (upregulated in 5xFAD compared with WT and downregulated in Trem2 R47H ; 5xFAD compared with 5xFAD).This gene is part of the ephrin family and is critically involved in memory formation, with knockout impacts of related receptors (Epha3/4) resulting in reduction in context dependent memory [68].
Most neuron types also exhibit genotype-speci c enriched and reduced subclusters, primarily in 5xFAD and Trem2 R47H ; 5xFAD mice.Several subclusters exhibit spatial localization, such as the L5 NP subcluster spatially localized near the subiculum, and the L2 IT subcluster spatially localized in the visual and retrosplenial cortices.
Together, our spatial transcriptomic analysis of the effect of the Trem2 R47H mutation on transcriptional dysregulation across both cortical and subcortical brain regions identi ed plaque and genotype dependent transcriptional alterations, cell type speci c transcriptome alterations, and the identi cation of genotype speci c cell sub types spatially localized across the brain., 5xFAD and Trem2 R47H ; 5xFAD mice.B: Integration of cell by gene matrix with RNA spatial location enables spatial analysis of transcriptomic variation on a regional and genotype basis.C: 300 gene overlay on a single coronal section, at increasing resolutions.D: UMAP displaying 37 annotated cell types after integration across all samples.E: UMAP of cell genotypes.Note the distinct subpopulations speci c to 5xFAD and Trem2 R47H ; 5xFAD genotypes, particularly in microglia and astrocyte cell populations.F: Hierarchical organization of cell clusters, combined with raw cell type proportions per genotype.

Figures
Figures