Integrative in situ mapping of single-cell transcriptional states and tissue histopathology in an Alzheimer’s disease model

Amyloid-β plaques and neurofibrillary tau tangles are the neuropathologic hallmarks of Alzheimer’s disease (AD), but the spatiotemporal cellular responses and molecular mechanisms underlying AD pathophysiology remain poorly understood. Here we introduce STARmap PLUS to simultaneously map single-cell transcriptional states and disease marker proteins in brain tissues of AD mouse models at subcellular resolution (200 nm). This high-resolution spatial transcriptomics map revealed a core-shell structure where disease-associated microglia (DAM) closely contact amyloid-β plaques, whereas disease-associated astrocytes (DAA) and oligodendrocyte precursor cells (OPC) are enriched in the outer shells surrounding the plaque- DAM complex. Hyperphosphorylated tau emerged mainly in excitatory neurons in the CA1 region accompanied by the infiltration of oligodendrocyte subtypes into the axon bundles of hippocampal alveus. The integrative STARmap PLUS method bridges single-cell gene expression profiles with tissue histopathology at subcellular resolution, providing an unprecedented roadmap to pinpoint the molecular and cellular mechanisms of AD pathology and neurodegeneration.


Introduction
Alzheimer's disease (AD) is a progressive neurodegenerative disorder and the most common cause of dementia in the elderly 1 . Widespread deposition of extracellular amyloid-β (Aβ) plaques and intracellular neurofibrillary tangles (hyperphosphorylated tau deposits), especially in the neocortex and hippocampus, are the neuropathologic hallmarks of AD [1][2][3][4] . In addition, AD pathology features reactive changes of microglia and astrocytes and white matter abnormalities [5][6][7] . A key question in AD research is how the histopathological hallmarks are correlated with molecular disturbance that drives neurodegeneration across different cell populations.
Conventional experimental approaches are disadvantageous for uncovering the molecular and cellular complexity of AD: Bulk-tissue analyses mask the heterogeneity of cell populations in the 3 brain; standard imaging methods visualize only a few genes/proteins and access limited cell types and brain regions. The recent application of single-cell RNA sequencing (scRNA-seq) to AD brain tissue has been transformative, revealing substantial and heterogeneous changes of gene expression in major brain cell types in patients and mouse models [8][9][10][11][12][13] . A sub-population of microglia with a distinctive transcriptomic state (termed disease-associated microglia or DAM) was identified by bulk and scRNA-seq studies of mouse AD models and other neurodegenerative disorders 10,[14][15][16][17][18][19] . Besides DAM microglia, disease-associated astrocytes (DAA) with characteristic gene signatures emerge in AD, suggesting a major transcriptional response in AD by multiple cell types 9,20-25 .
However, there are major limitations of scRNA-seq methods: they cannot preserve the spatial pattern of single cells (or their relationships to localized tissue pathology); and the isolation of single cells or single nuclei from the brain tissue can introduce significant bias in cellular representation and artifactual changes in gene expression 26,27 . Therefore, a fundamentally different technology platform capable of integrating spatially resolved single-cell transcriptomics with histology and immunostaining in intact tissue is needed to fully understand the scope and heterogeneity of diverse cellular responses to amyloid plaque, tau aggregation, cell death, and synapse loss, and to investigate the spatial relationships between the above localized pathologies and cellular responses. Successful development and implementation of this technology would be highly informative for AD research.
Spatially resolved transcriptomic technologies are capable of mapping single-cell transcriptomic profiles within tissue architecture [28][29][30][31][32] , but many existing ones are incompatible with protein detection in the same tissue sections or limited by the spatial resolution and/or gene coverage.
For example, a recent study uncovered amyloid plaque-induced genes (PIG) using Spatial Transcriptomics with fluorescent staining of adjacent brain sections to correlate the positions of 4 plaques with local gene expression 33 . However, the resolution was limited to 100 µm and only a small set of genes were verified at cellular resolution.
We previously developed an image-based in situ RNA sequencing method called STARmap (spatially-resolved transcript amplicon readout mapping) for single-cell transcriptional state profiling in three-dimensional (3D) intact brain tissues 34,35 . Here we introduce STARmap with protein localization and unlimited sequencing (STARmap PLUS), enabling simultaneous highresolution spatial transcriptomics concomitantly with specific protein localization in the same tissue section. We employed STARmap PLUS to draw a comprehensive transcriptomic atlas of AD at 200 nm resolution across all brain cell types during the development of amyloid plaque and tau pathology (Fig. 1a). In this study, we applied STARmap PLUS to TauPS2APP triple transgenic mice, an established mouse AD model that exhibits both amyloid plaque and tau pathologies.
By mapping a targeted list of 2,766 genes, we created a spatial cell atlas of 8-and 13-month old

5
The original STARmap method is incompatible with histological staining (immuno-or chemical staining) and limited to detecting 1,024 genes 32 . To overcome such limitations, in STARmap PLUS, we developed the experimental protocol to incorporate antibody staining (in this study, AT8 antibody, detecting hyperphosphorylated tau) and dye staining (X-34, detecting Aβ plaque) into the library preparation and in situ sequencing steps (Fig. 1b, Extended Data Fig. 1a-b). In brief, STARmap PLUS entails the following steps: 1) mRNAs within fixed brain sections are detected by a pair of SNAIL (specific amplification of nucleic acids via intramolecular ligation) probes (Fig.1b), and enzymatically amplified as cDNA amplicons; 2) specific proteins of interest are labeled with primary antibody; 3) the cDNA amplicons, primary antibodies, and endogenous proteins (e.g., plaques and tau) are chemically modified and copolymerized in a hydrogel matrix; 4) each cDNA amplicon contains a gene-unique identifier (barcode) that is read-out through in situ sequencing with error-reduction by dynamic annealing and ligation (SEDAL); 5) fluorescent secondary antibody and small-molecule dye X-34 are then applied to visualize specific proteins and their localization. Besides protein labeling capability, we expanded the gene-coding barcode in the DNA probes from 5 to 10 nucleotides (resulting in 10 6 coding capacity for STARmap PLUS) that is sufficient to encode more than 20,000 genes (Fig. 1b).
We applied STARmap PLUS to investigate how AD-related pathology, including amyloid deposition and hyperphosphorylated tau, influences brain cell states at the transcriptomic level at single-cell resolution in intact brain tissue in which spatial relationships between protein pathology, cell location and mRNA changes are maintained and readily measured. A curated list of 2,766 genes were extracted from previous bulk and single-cell RNA-seq studies and diverse AD-related databases (Supplementary Table 1) 8,9,12,13,[36][37][38][39][40][41] . We performed 8 rounds of in situ sequencing to map RNAs and 1 round of post-sequencing imaging to locate Aβ plaques and hyperphosphorylated tau (p-Tau) in coronal sections of the brains from TauPS2APP mice and control mice ( Fig. 1c and Extended Data Fig. 1c).
6 Transgenic TauPS2APP mice develop progressive amyloid plaques and tau pathology starting at 4.5 months of age, growing exponentially from 6 to 8 months of age, and rising steadily from 9 months of age. Prominent neurodegeneration, measured by disintegrative staining or volumetric MRI, becomes apparent by 9 months 34,35 . Thus, brain sections were collected and analyzed at 8 months (when tau and Aβ pathology have set and are expanding) and 13 months (a more advanced disease stage with severe pathology and elevated neuroinflammatory activity) 35 . We focused on the cortex and hippocampus, two of the most susceptible regions in AD. STARmap PLUS revealed that Aβ plaques were prominent in both cortex and hippocampus whereas AT8 p-Tau immunoreactivity was strongest in the CA1 region of the hippocampus (Fig. 1c, Extended Data Fig. 1c), which is consistent with previous reports 34,35 . A total of 19,932 cells from two TauPS2APP mice and 17,135 cells from two control mice (non-transgenic littermates) were imaged at 200 nm resolution (Fig. 1a). After projecting 3D RNA reads to two-dimensional planes, segmenting cells, and filtering single-cell transcriptional profiles by quality control (see Methods), the remaining 33,106 cells pooled from all 4 samples were subjected to downstream analysis (Extended Data Fig. 1d).

Hierarchical cell classification and spatial analysis
To identify cell types from the STARmap PLUS data, we adopted a hierarchical clustering strategy, where top-level clustering served to classify cells into common cell types shared by all samples, and sub-level clustering served to further identify disease-associated subtypes. During top-level clustering, the Leiden algorithm was applied to the Uniform Manifold Approximation and Projection (UMAP) representation of all transcriptomic profiles 42,43 . Through Leiden clustering, we identified 13 major clusters and annotated the cell types according to their spatial distribution and previously reported gene markers 38,39,41 (Fig. 2a and Extended Data Fig. 2a). For example, excitatory neurons were annotated by their high expression levels of genes related to ion channels and synaptic signaling such as Vsnl1, Snap25, and Dnm1. Inhibitory neurons were separated by 7 their enrichment of gamma-aminobutyric acid (GABA) transporters Slc6a1. Other non-neuronal cell type specific markers such as Aldoc for astrocyte, Bsg for endothelial cell, Ctss for microglia, Plp1 for oligodendrocyte were used to annotate corresponding clusters (Extended Data Fig. 2a).
The UMAP plots of TauPS2APP samples showed differential distribution of cells within the astrocyte, microglia, oligodendrocyte and dentate gyrus (DG) clusters in comparison with controls, suggesting possible disease-associated cell subtypes (Fig. 2a). Thus, we further investigated the transcriptomic heterogeneity within each major cell type, which further identified 27 sub-level clusters based on their transcriptomic signatures (Fig. 2b).
Because spatial information of RNAs was preserved at subcellular resolution, we were able to generate a high-resolution spatial cell atlas of the cortex and hippocampus in conjunction with precise localization of AD-related histopathology ( Fig. 2c and Extended Data Fig. 2b,c). All spatial cell atlases of TauPS2APP and control non-transgenic mice showed similar anatomic structure with major cell-type distribution in cortical and hippocampal regions, confirming the robustness, reproducibility and reliability of top-level clustering results. In TauPS2APP mice, the percentage of plaque area in the whole tissue (cortex and hippocampus) rose from 0.17% (73 Aβ plaques identified with an average size of 89 ± 44.9 µm 2 ) in 8-month mouse brains to 0.54% (151 plaques identified with an average size of 157 ± 110.7 µm 2 ) in 13-month mouse brains. Meanwhile, no plaques were visible in non-transgenic controls (Extended Data Fig. 1c). p-Tau signal was observed mostly in the hippocampal region of the TauPS2APP mice at 8 months and became more pronounced at 13 months in both hippocampal and cortical regions ( Fig. 2c and Extended Data Fig. 2c).
How are cell distributions and cell states affected by the nearby presence of amyloid plaque and tauopathy? First, given the different cell-type compositions in different parts of the brain, we analyzed the cell distributions and states in the cortex and hippocampus separately. Next, to quantify cell-type distribution in the spatial relationship to Aβ plaque, we calculated the density of 8 different cell types at the first five 10 µm intervals surrounding plaques (Fig. 2d, Methods). Among the 13 major cell types, microglia, astrocytes, oligodendrocytes, and oligodendrocyte precursor cells (OPC) showed relative enrichment in the vicinity of plaques in comparison with the overall density of that cell type in the cortex and hippocampus in TauPS2APP mice ( Fig. 2e and Extended Data Fig. 2d). The most striking change was detected in microglia: they were the most prevalent cell type in the 10 µm ring around the plaque, where they were markedly concentrated compared with the normal density in that brain region and often seemed to be in direct contact with the plaques. Astrocytes, oligodendrocytes and OPC showed modest but significant enrichment at the 10-30 µm distance compared to their overall average density but were depleted in the immediate neighborhood of plaques (<10 µm), perhaps crowded out by microglia. Both excitatory and inhibitory neurons were also depleted within 10 µm around Aβ plaques in comparison with the overall density of each neuron type.
In contrast to the localized macroaggregates of Aβ plaques, p-Tau (AT8) immunoreactivity was dispersed throughout neuronal cell bodies and axon bundles. To investigate the spatial correlation of p-Tau intensity with different cell types, we divided the spatial atlas into blocks of 20 µm-spaced grids ( Fig. 2f) and analyzed the covariation of cell-type composition versus p-Tau density.
Oligodendrocytes were found enriched in the blocks with high p-Tau intensity regardless of the presence or absence of Aβ plaques (Fig. 2g).
The top-level cell clustering and spatial analyses revealed that microglia, astrocytes, OPCs, oligodendrocytes, and neuronal cells showed the biggest changes in their transcriptional profiles, spatial distributions, or both. These cell types were thus selected for in-depth sub-clustering analyses to pinpoint disease-associated subtypes, cell states, and gene pathways in the following sections (Figs. 2b and 3-6).

Disease-associated microglia directly contact Aβ plaques from early disease stage
We first investigated the heterogeneity within the microglia population by sub-clustering analysis of transcriptomic profiles. Three subpopulations were identified and categorized as Micro1 (45%), Micro2 (24%), and Micro3 (31%) (Fig. 3a). Micro1 and Micro2 subtypes were present in both TauPS2APP and control samples, while the Micro3 subpopulation cells were almost absent in controls and expanded greatly from 8 months to 13 months in the TauPS2APP brains (Fig. 3a). It is noteworthy that clustering analysis effectively identifies distinct cell subtypes but does not capture well multi-step cell-state transitions. To uncover cell-state transitions during disease progression and determine the relationship among different subtypes, we deployed Monocle pseudotime analysis 44 , a widely used computational tool for reconstructing cell differentiation trajectory, as a complement to subtype analysis in the following sections. We aimed to reconstruct the presumptive path along which microglia alter their state by pseudotime trajectory analysis of the single cell transcriptomic profiles (Fig. 3c). By this computational approach, the microglia population showed a linear pseudotime trajectory that aligned well with the real disease progression timeline: microglia in control mice were enriched at the starting point of the trajectory while those in TauPS2APP mice shifted along the trajectory from 8 months to 13 months (Fig. 3c). To get a more comprehensive understanding of microglia response to AD pathology at the molecular level, we analyzed differential expressed genes (DEGs) of microglia along disease progression and in relation to plaques (Fig. 3i,j and Extended Data Fig. 3e,f, Supplementary Table   3). Most of the up-regulated microglial DEGs identified in the TauPS2APP mouse brains overlap with Micro3 (DAM) gene markers and with up-regulated "spatial DEGs" near plaques; these genes are involved in biological processes such as regulation of immune response (i.e. Apoe, C1qa, Cd74), regulation of cell activation (i.e, Cd9, Cd83, Cst7), cellular protein catabolic process (i.e. Ctss, Ctsb, Sumo2, Cul1) (Extended Data Fig. 3e,f). Notably, relating to a recent study of human AD patient samples 45 , our study revealed that a few upregulated genes (i.e. Ccl3, Trem2,Csf1r) in the DAM microglia of the TauPS2APP model may be mediated by the activation of ERK1/ERK2 cascade (Extended Data Fig. 3f).

Disease-associated astrocytes emerge near the plaque-DAM complex at a later stage
Astrocyte was another non-neuronal cell type that showed a significant difference in TauPS2APP versus control. Sub-clustering analysis of the astrocytes identified three transcriptionally distinct subpopulations Astro1 (37%), Astro2 (44%), and Astro3 (19%) (Fig. 4a). Astro3 cells were rarely found in control mice but the Astro3 subpopulation greatly expanded from 8 months to 13 months in TauPS2APP mice (Fig. 4a). Transcriptomic changes in Astro3 astrocytes closely resembled those of disease-associated astrocytes (DAA) previously described in AD models, characterized by upregulation of genes such as Gfap, Vim, Apoe (Fig. 4b) 9 . Based on this transcriptomic similarity and their association with TauPS2APP mice we annotated the Astro3 subtype as DAA.
Astro1 and Astro2 are present in both control and TauPS2APP mice; their gene markers correspond with previously reported low-and intermediate-Gfap cell populations 9 . In comparison with controls, the overall cell density of Astro1 declined in the TauPS2APP mice whereas the density of Astro2 increased in the cortex.
Despite the apparent linear gradient of many 'marker' genes across Astro1-3 subtypes, we noted a bifurcation path based on pseudotime trajectory analysis of the astrocyte population (Fig. 4c).
By visualizing the subtype annotation along the pseudotime trajectory, we identified that the longest path (lower path) matches with the gradient from Astro1 (the starting point) and  Table 3). The gene expression profiles of top DEGs identified on the pseudotime embedding showed that the expression pattern of Vim best correlates with the DAA population, whereas the molecular gradients of Gfap resembled the disease-associated trajectory from Astro 2 to Astro 3 ( Fig. 4i). Spatial analysis of gene expression around plaques in the 13month TauPS2APP mice further support the spatial enrichment of DAA near plaques (Fig. 4j).

OPCs are enriched in the intermediate vicinity of Aβ plaques and oligodendrocyte subtypes infiltrate the hippocampal alveus colocalizing with p-Tau.
Sub-clustering analysis identified four subtypes within the oligodendrocyte lineage ( Fig. 5a Table 2). At 13 months, the density of Oligo1 increased by 58%-80% at the 10-30 µm distance from the plaques in TauPS2APP mice in comparison with the overall density ( Fig. 5h); the average densities of Oligo2 and Oligo3 also increased by 1.5-3 fold at the 10-30 µm 14 distance from plaques but they were only sparsely present near (10-30µm) 13-14% of all plaques. Spatial analysis of the pseudotime values of oligodendrocytes lineages (mostly Oligo1) showed a global upregulation from 8 months to 13 months regardless of their distance to plaques (Extended Data Fig. 5c-e), indicating that the shift of cell states in oligodendrocytes is not spatially confined near plaques.
Given that the density of oligodendrocytes is positively correlated with the density of p-Tau ( Fig.   2g), we sought to pinpoint which oligodendrocyte subtype is spatially associated with tauopathy.
Using the aforementioned grid-based spatial correlation analysis, we found that in the TauPS2APP mice, the cell density of Oligo1 and Oligo3 increased by 1.5-3 fold in regions with higher p-Tau signals regardless of Aβ pathology ( Fig. 5i and Extended Data Fig. 5f). In contrast, the cell density of OPC and Oligo2 increased by 30% and 140% respectively only in the presence of plaques, suggesting that these oligodendrocyte lineage populations respond to amyloid rather than Tau pathology. We observed that p-Tau signals were concentrated in the alveus of the hippocampus where axon bundles of hippocampal neurons course (alveus, Extended Data Fig.   5g). Compared with 13-month control samples, the Oligo1 and Oligo3 densities increased markedly (2-4-fold) in the alveus region of 13-month TauPS2APP mice (Fig. 5j). In total, the spatial analysis of cells of the oligodendrocyte lineage in relation to p-Tau revealed a strong association between tauopathy and the accumulation of Oligo1 and Oligo3 subtypes in the hippocampus alveus.
Through analyzing the DEGs in oligodendrocytes along disease-associated pseudotime trajectory and around plaques, as well as DEGs from the TauPS2APP mice versus control mice, we identified and verified a group of genes (such as Klk6, C4b, Cd9, Serpina3n) that were strongly upregulated in 13-month old TauPS2APP mice compared to the same age control mice (Supplementary Table 3). Interestingly, oligodendrocytes near plaques exhibited increased expression of these genes, suggesting subpopulation(s) of oligodendrocytes (with elevated expression of C4b, Serpina3n, Cd9) might be responding to plaques or interacting with other cells or cellular compartments (microglia, astrocytes and/or dystrophic neurites) that are affected by amyloid pathology (Fig. 5l). In contrast, the fold change of DEGs between 8-month old TauPS2APP and control mice was less significant (Supplementary Table 3 Table 3). Dapk1 is involved in neuronal death regulation whereas Tnik has been implicated in dentate gyrus neurogenesis 46 , so these results may indicate altered neurogenesis in DG in TauPS2APP mice, which is related to a recent report showing that the adult hippocampal neurogenesis activity in DG sharply declines in human patients of AD 47 .
To investigate the neuronal transcriptomic alterations induced by tauopathy, we first quantified the ratio of p-Tau positive pixels to the total pixel area of neuronal cell bodies and defined p-Tau positive neurons as those in which this ratio is greater than 7% (see Methods). In TauPS2APP mice, there were 4 times more p-Tau positive neurons at 13 months than 8 months. At 8 months the majority of p-Tau positive neurons were CTX-Ex2 excitatory neurons, whereas at 13 months the majority of p-Tau positive neurons were the CA1 excitatory neurons. Inhibitory neurons account for <20% p-Tau positive neurons: most of them were Pvalb neurons at 8 months whereas at 13 months most of them came from the Sst population (Fig. 6i). We also quantified the p-Tau signal around the plaques and found that the p-Tau signal was enriched within 10 µm distance near plaques (Fig. 6h). Considering neuronal cell bodies were relatively depleted within the 10 µm range, the observed p-Tau signals likely correspond to dystrophic (injured) neurites 35 .
We finally examined the DEGs of all neuronal types in response to Aβ plaques and p-Tau. GO term analysis of general neuronal DEGs (TauPS2APP versus control) revealed upregulation of protein kinase activity, downregulation of protein transport and section, and downregulation of mRNA processing (Extended Data Fig. 6g). The DEGs of the dentate gyrus region were related to downregulation of brain development (Extended Data Fig. 6h). SynGO analysis 36 of DEGs in p-Tau positive CA1 neurons at 13 months are clustered in synaptic vesicle cycle, synaptic vesicle endocytosis and postsynapse organization (Fig. 6j). Neuron-type resolved DEG analysis (Supplementary Table 3) revealed that Ccnb2 (G2/mitotic-specific cyclin-B2) was identified among top upregulated DEGs from all neuronal types in TauPS2APP models, as well as genes involved in cell death (i.e. Ddit3, Dapk), suggesting a reactivation of the cell cycle and supporting the notion that neuronal reentry into the cell cycle is a mediator of neurodegeneration and cell death 48,49 .

Integrative analysis of disease-associated cells and genes in AD pathology
The analyses above focused on dissecting disease-associated subtypes and DEGs within major brain cell types. In order to synthesize a comprehensive picture of AD gene pathways from multiple cell types, we performed cell-type-resolved Gene Set Enrichment Analysis (GSEA) 50,51 using DEGs from four major cell types (microglia, astrocytes, oligodendrocytes, and neurons) comparing TauPS2APP versus control. The results (Fig. 7a,b) showed that most of the significantly enriched terms of upregulated genes in non-neuronal cell DEGs were in the biological processes of inflammatory response, glial cell differentiation, gliosis, cytokine production, and regulation of neuronal death, while those of neurons involve synaptic signaling and protein kinase activity. Neurons in the TauPS2APP mice also showed significant downregulation of RNA processing/splicing and protein transport (Fig. 7a,b).
Besides the DEGs defined by comparison of TauPS2APP vs. control samples, we also identified plaque-associated spatial DEGs (SDEGs) by clustering genes based on their expression profiles as a function of distance from amyloid plaques (Fig. 7c,d). We identified six clusters of SDEGs:  Table 4). We identified 57 PIGs at 8 months and 81 PIGs in the TauPS2APP sample at 13 months (Fig. 7g). 37 out of 57 PIGs in the 8-month TauPS2APP sample were also found in the 13-month samples (Fig. 7g). The PIGs enriched within 10 µm distance from plaques were mainly DAM marker genes, such as Trem2, Cst7, Ctsb, Apoe and Cd9 (Fig. 7e,f). The DAA marker Vim was upregulated in the regions 10 to 20 µm away from plaques in AD 13-month samples (Fig. 7f), which is consistent with our previous finding that DAA enriched in the 10-20 µm ring (Fig. 4). Although SDEGs in Cluster 4-6 (Supplementary Table 4) are down-regulated at the 0-30 µm distance from plaques and not defined as PIGs, they are also valuable to reveal potential disease mechanisms. For example, an interesting set of SDEGs was induced at 8 months with peak expression at 30-40 µm from plaque, including Elavl4, an RNA binding protein that is implicated in early dysregulation of splicing in frontotemporal dementia 52 . Comparing our PIGs with previously reported PIGs identified from 18month AppNL-G-F mice 33 , 22 out of 81 and 19 out of 57 PIGs in our AD 8-month and 13-month samples, respectively, overlapped with previously reported PIGs (Fig. 7g).
To validate the findings from the study, we repeated the STARmap PLUS experiment on another set of TauPS2APP mice, this time using a focused subset of genes (64 key cell-type marker genes and disease-associated genes). This focused 64 gene analysis yielded cell clustering results, DEG data and spatial map information that was highly consistent with the major conclusions derived from the 2766-gene datasets (Supplementary Tables 1,2, Extended Data Fig. 7) including the spatial enrichment of various glial cells surrounding plaques (Extended Data Fig. 7d), the infiltration of oligodendrocytes in hippocampus alveus (Extended Data Fig. 7e), as well as celltype and spatially resolved DEGs (Extended Data Fig. 7f-h).

Discussion
Here we developed STARmap PLUS for in situ detection of RNAs and proteins in the same tissue section at subcellular resolution. Compared to STARmap, we have substantially improved the number of in situ sequenced genes and enabled a simultaneous profiling of RNAs and proteins in intact hydrogel-tissue scaffolds. This method development enables a multimodal analysis of spatially expressed RNAs and proteins, bringing an opportunity to study biological systems in a more comprehensive manner. This method is readily applicable to a variety of samples for study of physiological and pathological mechanisms in healthy and diseased tissues by integrative mapping of single-cell transcriptomic states, tissue morphology and disease markers.

Analysis framework
We developed STARmap PLUS to enable high-resolution spatial transcriptomics in combination with microscopic protein localization. In a mouse model of AD exhibiting amyloidosis, tauopathy and neurodegeneration we applied this high-resolution and multimodal in situ mapping with four computational analysis strategies to identify disease-associated cell populations and gene programs: (1) hierarchical cell clustering analysis to pinpoint disease-associated cell subtypes, (2) pseudotime trajectory analysis to reconstruct the transitions of cell states during disease progression, (3) spatial analysis to recognize changes in cell types and cell states in physical proximity of Aβ plaques and tauopathy; (4) differential gene expression analysis in spatiotemporal relationship with amyloid and p-Tau and disease stage to identify disease-associated gene pathways. This analysis framework successfully traced how AD hallmark pathologies correlate with gene pathways that drive inflammation, gliosis, and neurodegeneration across different cell types. 20

Core-shell structures surrounding plaques and cascades of cell-state transitions
Applying our methods to the TauPS2APP mouse AD model at two different ages, we have constructed a cell-type and cell-state resolved spatiotemporal map of the TauPS2APP mouse model (schematized in Fig. 7i). Specifically, in response to the early emergence of Aβ plaques (at 8 months or earlier), microglia are the primary responders, closely aggregating around plaques (<10 µm distance away from plaques, Fig. 7i). Integrative spatial and pseudotime analysis of microglia subtypes suggested a cell state transition from Micro1 and Micro2 to Micro3 subtype (DAM) associated with microglia cells accumulating around Aβ plaques ( Fig. 3 and Extended Data   Fig. 3). The transition from Astro2 to Astro3 (DAA) in the layer (10-20 µm) next to the DAM microglia was seen at a later stage (13 months), which is consistent with an induction of DAA by DAM. Indeed, a previous study showed that reactive astrocytes with the high expression of DAA marker genes Gfap and Vim, were induced by activated microglia 21 . Like DAA, Oligo2/3 accumulation near plaque (10-30 µm) occurred later (13 months) (Fig. 7h). Interestingly, we discovered that OPCs were enriched in the regions 10-30µm away from plaques at both 8-month and 13-month stages, indicating potential in situ proliferation and differentiation of OPCs to oligodendrocytes. In contrast to glia, neuron density around plaque declined from 8 months to 13 months (Fig. 7h,i).
Based on the collective data, we propose a core-shell structure of glial cells surrounding Aβ plaques where the DAM emerge early in disease near plaques as the core, and the shell is a gliosis zone enriched for DAA, OPC, and oligodendrocytes that develops at a later disease stage, perhaps dependent on formation of the inner ring of reactive microglia. Because STARmap PLUS only captures snapshots from different disease stages, the dynamic cell-type and state transitions inferred from pseudotime trajectory reconstruction and spatial patterns need future verification by live-cell imaging or in vivo cell-fate tracing approaches.

Oligodendrocytes and tauopathy
In the TauPS2APP model, hyperphosphorylated tau (as detected with AT8 antibody) was mainly found in CA1 excitatory neuronal bodies and its upper layers of axon tracts (alveus), where it is strongly associated with the infiltration of Oligo1 and Oligo3 subtypes regardless of the presence or absence of plaques (schematized in Fig. 7i). It is unclear whether the infiltration of Oligo1 and Oligo3 is a reactive mechanism to support tauopathy-injured axons, repair damaged myelin, or whether Oligo1 and Oligo3 instead exacerbate tauopathy. It is possible that Oligo1 and Oligo3 might feature distinct responses and exert different impacts, for example, Oligo 1 being homeostatic and protective while Oligo 3 is detrimental in response to tau-driven changes. The reciprocal functional interactions between oligodendrocytes and neurons (beyond oligodendrocytes providing electrical insulation and metabolic support for axons) are increasingly appreciated and studied 53 . Single-cell RNAseq studies of human AD brain tissues have noted major transcriptomic changes in oligodendrocytes, but without spatial information in relation to plaque and p-Tau 12 . Enabled by high-resolution spatial transcriptomic analysis of oligodendrocytes in the TauPS2APP mice, we identified oligodendrocytes subtypes that are associated with tauopathy. Further pathway analysis implicated oligodendrocyte DEGs in the regulation of neuron death, cytokine production, synaptic plasticity, and anterograde synaptic signaling (Extended Data Fig. 5h), pointing to a potential link between inflammation and axonal tauopathy mediated by oligodendrocytes.

Implications for neurodegeneration mechanisms in AD
The overall density of neuronal cells in cortex declines near Aβ plaques (Fig. 6c,d), implying a loss of neurons due directly to Aβ toxicity or indirectly via the effects of amyloid on microglia and other glia surrounding the plaque. This notion is supported by the accumulation of p-Tau in dystrophic neurites near to plaque. We also cannot exclude that neurons are merely displaced 22 from plaque-adjacent regions by reactive glia; however, we note that the TauPS2APP mice show macroscopic brain volume loss by 13 months of age 35 . Consistent with neuronal loss around plaques, we note that neurons in DG, where a large amount of Aβ plaques appeared by 13 months, showed altered transcriptional profile in intracellular protein transport, RNA processing and cell cycle pathways related to cell stress and neurodegeneration as well as molecular pathways related to synapse structure and organization. Meanwhile, all neuronal types examined in cortex and hippocampus of TauPS2APP share some common DEGs (i.e. Ccnb2, Ddit3), perhaps suggesting a general mechanism of cell cycle reactivation in neurodegeneration. The core-shell glial structure around plaques implies that microglia may cause neurodegeneration in part through their crosstalk with astrocytes 21   Iodide (PI) staining was performed following the manufacturer's instruction for the purpose of cell segmentation. Another round of imaging was performed to detect spatial protein signals.

STARmap PLUS Image Processing
All of the image processing steps were implemented using MATLAB R2019b and related opensource packages in Python 3.6 and applied according to 32 .

Image Preprocessing
A multi-dimensional histogram matching was performed on each tile with MATLAB function `imhistmatchn`. It used the image from the first color channel in the first sequencing round as a reference to uniform the illuminance and contrast level of all other images in the current imaging position. Additionally, a customized tophat filtering was applied to the sequencing images to further enhance the signal and suppress the background noise.

Image Registration
Image registration was applied according to 32 . Global image registration was accomplished using a three-dimensional fast Fourier transform (FFT) to compute the cross-correlation between two image volumes at all translational offsets. The position of the maximal correlation coefficient was identified and used to translate image volumes to compensate for the offset. Then a non-rigid registration was applied with MATLAB function `imregdemons` to further align images in different sequencing rounds.

Spot Calling
After registration, individual dots were identified separately in each color channel on the first round of sequencing. For this experiment, amplicon dots were identified by finding local maxima in 3D with MATLAB function `imregionalmax`. Dots with intensity at their centroids less than the threshold were removed. Because the dot was approximately 6 pixels in diameter on the xy plane, the dominant color for that dot across all four channels on each round was determined by a 5x5x3 voxel volume surrounding the dot centroid. The integrated intensity of the voxel volume in each channel was used for color determination. In this case, each dot in each round had a L2 normalized vector with four elements. The color of each dot was determined by the corresponding channel with the highest value in the vector. Dots with multiple maximum values in the vector were discarded.

Barcode Filtering
Dots were first filtered based on quality scores (average of -log(color vector value in dominant channel) across all sequencing rounds). The quality score quantified the extent to which each dot on each sequencing round came from one color rather than a mixture of colors. The barcode codebook was converted into color space based on the expected color sequence following the 2base encoding of the barcode DNA sequence. Dot color sequences that passed the quality 47 threshold and matched sequences in the codebook were kept and identified with the specific gene that that barcode represented; all other dots were rejected. The high-quality dots and associated gene identities in the codebook were then saved for downstream analysis.

2D Cell Segmentation
Nuclei were automatically identified by applying a pre-trained 2D machine learning model from the StarDist package 54 to a maximum intensity projection of the stitched DAPI channel following the final round of sequencing. Then the cell locations were extracted from the segmented DAPI image and used as markers for cell body segmentation. Cell bodies were represented by an overlay of stitched Nissl staining and merged amplicon images. Finally, a marker-based watershed transform was applied to segment the thresholded cell bodies. Points overlapping each segmented cell region in 2D were then assigned to that cell, to compute a per-cell gene expression matrix.

Protein image preprocessing
p-Tau images were processed using a customized Fiji macro 55 . A rolling ball background subtraction with a radius equal to 5 was applied to each image. Then each image was processed with a Gaussian filter with a sigma value equal to 2 and followed by a maximum entropy binarization.

Cell type classification
A two-level clustering strategy was applied to identify both major and sub-level cell types in the dataset. Processing steps in this section were implemented using Scanpy v1.4.6 56 and other customized scripts in Python 3.6 and applied according to Wang et al (2018) 32 . After filtration, normalization, and scaling, principal component analysis (PCA) was applied to reduce the dimensionality of the cellular expression matrix. Based on the explained variance ratio, the top 48 PCs were used to compute the neighborhood graph of observations. Then the Leiden algorithm was used to identify well-connected cells as clusters in a low dimensional representation of the transcriptomics profile. The cells were displayed using the Uniform Manifold Approximation and Projection (UMAP) and color-coded according to their cell types. The cells for each interested toplevel cluster were then extracted and sub clustered using the same workflow again.

Pseudotime and Trajectory Analysis
An R package Monocle3 44 (https://cole-trapnell-lab.github.io/monocle3/) was utilized for pseudotime calculation. Monocle3 firstly learned a principal graph (via learn_graph() function) in the three-dimensional UMAP space using a dimensionality-reduced representation of original data. In detail, Monocle3 ran a k-means clustering algorithm and selected "landmark cells" by first mapping cells to their nearest k-means points, then choosing the cell with highest local density as landmark cell. A principal graph was learned within these landmark cells rather than the whole dataset to reduce running time using simplePPT algorithm 57 !"After generating the principal graph, the graph was refined by initially performing a depth-first visitation of the nodes in the graph. For nodes with degree (the number of edges connected to this node) ≤ 2, no operation was performed.
For nodes with degree > 2, Monocle3 computed the diameter path for each subtree of the node not visited, and removed the subtree if its path length is less than a dedicated length specified by the user.
To compute pseudotime from the principal graph, firstly all cells were mapped to their nearest principal points based on Euclidean distance in the UMAP embedding. Next, for each edge in the graph, all cells mapped to the endpoints were orthogonally projected to the nearest points on the graph edge. Based on the projected position of cells, suppose the endpoints are a and b and here we have cells named ci , cj , with projection p(ci ) , p(cj ) . We can order the cell projections on the edge and suppose the order is a < p(ci ) < p(cj ) < b . Then edge (a,ci ) and (b,cj ) will be added to the final graph. After specifying the root node(s) by user, the pseudotime for each cell can be assigned by calculating the geodesic distance from root node to the cell node. In this case, raw expression matrix of each cell type: Microglia, Astrocyte, Oligodendrocyte/Oligodendrocyte precursor cell (OPC) and Dentate gyrus (DG), was used separately as the input of Monocle3. To visualize the results, we used the Monocle3-generated UMAP embedding with palettes representing normalized pseudotime (rescaled to 0 ~ 1), cell types, and sample identities.

Plaque Segmentation
The plaques were segmented from the thresholded image of plaque channel by using the `bwlabel` function in EBImage package 58 . Then, the size and the center of each plaque were calculated by using `computeFeatures.moment` and `computeFeatures.shape` functions respectively. Finally, plaques with an area greater than 400 pixels (~ 40 µm 2 ) were kept for the downstream analysis.

Cell Composition around Plaques
As filtered plaques were acquired in the last step, we dilated the plaque images 5 times with steps of 10 microns. Next, we counted the number of cells for every cell type that fall into different intervals. The ranges were set from 0-10µm (Ring 1) to 40-50µm (Ring 5). The statistics were normalized by calculating the percentage and density of each cell type in a ring. The graphical explanation of this analysis is shown in Fig. 2d. For the overall statistics, we calculated the percentage of each cell type in the whole sample. To test the significance of cell enrichment in both density and proportion aspects in specific intervals, two kinds of statistical tests: one sample t-test and chi-square test (Fisher's Exact test) were utilized. One sample t-test was used to compare the density values for each cell type in one interval around every plaque with the average density. According to the difference of mean value of density sequence and the average density, 50 an alternative hypothesis was set correspondingly to test whether the density was significantly higher or lower. For the chi-square test, raw cell counts in each interval were used for the test.

Differential Expression Analysis
Before performing DE analysis, the dataset was normalized according to the following steps: 1.
Divide the gene counts in a sample by the median of total counts per cell for that sample and multiply by the scale factor, which was defined as the mean value of median of total counts per cell for all samples; 2. Perform log2 transformation by adding a pseudo-count of one. DE genes were identified by performing Wilcoxon-Rank Sum test between two groups of cells using the `FindMarkers` function in Seurat. For 'disease vs. control' comparison of specific cell types, the two groups of cells were extracted from TauPS2APP and control samples and compared. As for the comparison between cells close to plaque versus cells far away from plaque, the cell was defined as 'near plaque' if its distance to the nearest plaque was less than 25µm and all other cells with distance to plaque greater than 25µm were defined as 'away from plaque'. In the comparison of 'CA1 Tau+ vs Tau-', Tau+ CA1 cells were defined according to the fraction of Tau signal area to the cell body's area. The threshold was set to 0.07.
In order to filter out lowly expressed genes, genes that were expressed by less than 10% of cells in either group of the comparison were excluded. We also applied the following threshold values on the generated gene list to filter out non-significant genes: absolute value of log fold change > logFC > 0 were colored in red while others were colored in blue. For those significant genes (pvalue < 0.05) but failed to pass the LogFC threshold were green tinted. All other non-significant genes were colored in gray. Note, some genes with extremely high -log(P-value) or logFC were capped.

Gene Set Enrichment Analysis
Gene Set Enrichment Analysis (GSEA v4.1.0) 50,51 was used to perform pathway analysis for each cell type with a LogFC ranked gene list whose genes were detected in a minimum fraction of 5% cells in the targeted cell type. Gene Ontology Biological Processes (GO:BP) database was selected as the gene sets database to perform GSEA. To limit the size of gene sets subjected to enrichment analysis, the minimum and maximum size of gene sets were set to 20 and 1000.
Number of permutations was set to 1000. Results were exported in tab-delimited format to be used as input for further functional enrichment analysis. The SynGO enrichment tool 36 was used to further characterize the synapse functions enriched in DEGs from excitatory neurons, inhibitory neurons, and CA1 cells with Tau pathology. Brain expressed genes were used as the background gene list.

Grid Analysis
To analyze the scattered p-Tau signal, we firstly performed rasterization on p-Tau images with resolution at 20µm (See heatmap Extended Data Fig. 5f). For each block, we counted the number of cells for each major cell type and the p-Tau signal intensity. Based on the distribution of p-Tau intensity and the existence of plaque, blocks were classified into 4 groups: No p-Tau, Low p-Tau (intensity of p-Tau in 1 to 50% range of population with non-zero p-Tau signal), High p-Tau w/ plaque (51-100% range with plaque) and High p-Tau w/o plaque (51-100% range without plaque).
The significance of cell type enrichment was tested using one-sample t-test comparing the density values of each group with average density.
The averaged expression matrix was calculated separately for disease samples in 8 months and 13 months. Next, K-means clustering was applied on the matrix respectively and 6 clusters were obtained in each sample. We tested whether the averaged expression level of the cells within the high expression interval was significantly higher than that of other cells outside the interval using Wilcoxon signed-rank test (With p-value < 0.05, Supplementary Table 4). Plaque induced genes (PIGs) were defined as significantly enriched genes within the 0-30 µm distance surrounding plaques (Cluster 1-3). 'ComplexHeatmap' package 59 is utilized for heat map visualization.

Quantification and Statistical analysis
The statistical tests and number of independent replicates per experiment are indicated in the One-sample t-test was used to test the significance of cell type enrichment by comparing the density values of each cell type in each interval with average density.

Code availability
The code for the STARmap PLUS image analysis is available on Zenodo (https://zenodo.org/record/5842625#.YeAvYljML0o)