A framework for quantitative multiplexed imaging and feature selection in human brain
Archival, formalin fixed paraffin embedded (FFPE) brain regions are predominantly used for human neurodegenerative disease research. Fluorescence imaging on adult brain FFPE is confounded by endogenous tissue autofluorescence that can overwhelm signals originating from antibody bound fluorophores (Additional file 1: Fig. S1A) [8, 9, 19]. With this in mind, we created a framework for high resolution, quantitative, multiplexed imaging of archival brain regions (Fig. 1a) using MIBI-TOF. Because antibodies are detected using elemental mass tags, MIBI-TOF images possess no equivalent endogenous signal from the biological matrix (i.e., no equivalent ‘autofluorescence’). Brain tissue sections were stained with a cocktail of primary antibodies (36-plex panel), where each antibody is labeled with a unique elemental mass tag (Fig. 1a, Additional file 1: Fig. S1B, Table S1). Stained brain sections were imaged using an ion beam, which liberates these mass reporters as secondary ions (Fig. 1a). The spatial distribution of elemental reporters is converted into an N-dimensional image where each channel of this image corresponds to one of the primary antibodies (Fig. 1a). Quantitative pixel-level information was extracted to produce global expression summaries and fed into downstream pipelines for data-driven analysis of cell phenotypes and protein aggregate classification (see Methods: Global expression pattern) [10]. Ultimately, these data were used to construct tabular summaries of molecular features describing overall tissue architecture, as well as the spatial distribution of cells and proteopathies (Fig. 1a).
To capture salient neuropathologic changes, the antibody staining panel targeted proteins for delineating major neuronal cell lineages, proteopathies, and vascular structures (Additional file 1: Fig. S1B, Table S1). Antibody specificities were validated against multiple brain regions by standard single-plex IHC (Additional file 1: Fig. S1C), benchmarking optimal titers and emphasizing the concordance of chromogenic IHC with MIBI-TOF imaging, like seen on our previous validation study with multiple other tissue organs [16]. MIBI spectral images may not completely recapitulate the visual features of IHC, due to the differences in chromogenic light capturing vs quantitative ion count-based construction of MIBI images. Particularly, in this study we only scanned a single depth of 100–300 nm tissue [16, 20]. Thus, cell type markers which meander in CNS tissue, such as Iba1, may have been partially captured from the 5 µm thick FFPE tissue. However, given we are capturing a large area of tissue in the x and y domains these markers still give valuable information about how different cell types interact locally. We used the resulting 36-plex panel to reveal the cytoarchitecture of an entire coronal section of human hippocampus (Fig. 1b–e). Hippocampus consists of anatomically distinct regions with structural boundaries that are evident from pseudo-colored overlays. Figure 1b–e illustrates a 31 mm2 hippocampal region from a person with AD dementia (ADD, Additional file 1: Table S2B). The granule cell layer of dentate gyrus (DG) is highlighted by MAP2 and HH3 (Fig. 1b, DG white arrow), while the boundary between white and gray matter can be delineated based on expression of MBP and CD56 (Fig. 1b). In line with previous work, glutamatergic terminals within the hilus, Cornu Ammonis (CA), and molecular layers (Lmol) are marked by the presynaptic vesicular proteins VGLUT1 and VGLUT2 (Fig. 1b) [21,22,23]. Calbindin (CB), calretinin (CR), and parvalbumin (PV) appear as contiguous parcels in the hippocampus (Fig. 1d). As seen in rat hippocampus [23], our human data showed that CB borders the DG and hilus, CR structures the Lmol, and PV is within CA4-CA1 [24]. As expected, NFTs and NTs are highlighted by PHF1-TAU immunoreactivity in the CA1 subfield, while Aβ plaques depicted by Aβ42-immunoreactivity predominantly in hilus and Lmol of the ADD hippocampus (Fig. 1e). Spectral images of all simultaneously acquired markers are illustrated in Additional file 1: Fig. S4I. Taken together, this comprehensive survey of brain phenotypes highlights the capability of our approach to visualize multiplexed signatures that reveal salient structural hallmarks and proteopathies in human hippocampus.
Multiplex imaging protein signatures organize anatomical structures
Given the visual differences in expression across hippocampus subfields, we next sought to validate the strength of the 36-plex panel in a quantifiable manner across a broader set of brain regions (Fig. 2). We used the 36-plex panel to analyze a tissue microarray (TMA) containing six different brain regions from cognitively normal individuals (Fig. 2a; Additional file 1: Table S2A). A total of 24 FOVs (4 per brain region) was acquired at low-resolution scan resolving ~ 1 µm2/pixel. A subset of markers (25-plex), namely phenotypic and structural targets, was selected for unsupervised clustering of the FOVs and markers. Firstly, except for medulla oblongata (MO), unsupervised clustering organized the different FOVs into similar anatomical brain regions (Fig. 2b). Secondly, mean pixel intensity z-score heatmap shows that many markers exhibited noticeable gradients that are consistent with canonical brain region compositions. For instance, calcium-binding proteins (PV, CB, and CR), that play important roles in memory processes, show laminar distributions and densities within and across the brain regions (Fig. 2b, c). These patterns look like published protein and gene expression profiles in mouse, rat and human, [25,26,27]. Notably within the cerebellum cluster, we observed an enrichment of CB with inhibitory VGAT, and CR with excitatory VGLUT2 markers possible highlighting presence of Purkinje cells in concordance previous reports on human and non-human primate cerebellum (Fig. 2b, c) [28, 29]. Likewise, TH and GAD65/67 are enriched in FOVs for SN and LC marking presence of dopaminergic neurons (Fig. 2b, d). The lower z-scores across most markers of the CA1 brain region can be explained by rarefaction of neuropil in the TMA core. Furthermore, FOVs from tissue sample of a cognitive impaired individual (CA1-ADD, CA1- AD dementia) were grouped in an unsupervised manner from normative (CA1), with these 25-plex phenotypic and structural targets (Additional file 1: Figs. S2B, S2C). Suggesting, cytoarchitectural alterations (Additional file 1: Fig. S2B), apart from gain of pan AD disease markers like PHF1Tau and Aβ, in the CA1-ADD tissue sampled (Additional file 1: Fig. S2C). Altogether, much like multiplexed analysis of the immune system and other tissues [10, 13, 15, 30,31,32,33], combinatorial protein expression patterns provided a snapshot of functional organization or proteopathy within brain regions. This low-level analysis of multiplexed proteomic data could serve as a guide for ‘fingerprinting' human brain and be used to model progression in neurodegeneration.
Discretized cellular and pathological features identify lineage and disease pathology-specific subclusters
We next discretized cells and proteopathies as biological units: neurons, astrocytes, microglia, vasculature, tau neurofibrillary tangles-neuropil threads (NFT-NTs), and Aβ plaques. Segmentation of planar brain tissue sections has inherent problems due to shape, texture, the disjointed features from neuronal and glial processes in the tissue section, but distant from their cell bodies, and from nonconforming shapes defining protein aggregates. As shown in Fig. 3a (dashed lines), we observed these various facets of planar imaging by both standard IHC and MIBI-TOF. To capture these attributes, two segmentation methods were used to partition neuronal perikaryons, microglia, astrocytes, and their processes, endothelial and their vascular-boundaries, Aβ plaques, and NFT-NTs (Fig. 3b). An adapted version of DeepCell was used for nuclear segmentation (Fig. 3b) [34,35,36]. For features not associated with nuclei, a pixel intensity thresholding-based method was used to partition cell body microglial and astrocyte projections, larger vascular structures, Aβ plaques, and NFT-NTs (Fig. 3b, see Methods: Object segmentation).
The same 36-plex stain of ADD hippocampus from Fig. 1b was chosen for the discretization analysis to capture a greater breadth of morphological and pathological landscapes. With nuclear segmentation, we identified 15,270 nuclei-associated cells across 85 connected imaging FOVs which were classified using manual gating (Additional file 1: Table S4, see Methods: Manual gating ADD hippocampus). With the pixel-based approach, we identified 19,332 features not associated with nuclei that included microglia, astrocytes, endothelial cells, Aβ plaques, and NFT-NTs (Additional file 1: Table S4). Neuron masks identified by nuclear segmentation were integrated with masks of the object segmented data (microglia, astrocytes, vasculature, Aβ plaques, and NFTs-NTs) to obtain a comprehensive repertoire of neuronal and non-neuronal cell types as well as disease features. We then generated a UMAP [37] plot organized by lineage and proteopathy markers to analyze the relationship between these features (see Methods: 1024 × 1024 ADD Images). Extracted cellular and proteopathy features mapped to four cell-type clusters (neurons, microglia, astrocytes, and vasculature) and two pathologic component clusters (Aβ plaques and NFTs-NTs) (Fig. 3d, inset).
NFT-NT masks (cyan) mostly interrelated with the neuron (dark red) cloud (Fig. 3d, Additional file 1: Fig. S3K). This is consistent with known accrual of NFT-NTs in intraneuronal compartments [6, 38]. Aβ plaques shared the UMAP space mainly with astrocyte (dark green, Fig. 3d, Additional file 1: Fig. S3L) and neuronal populations (Fig. 3d, Additional file 1: S3M) as well as scattered with microglia (dark yellow) and vascular groups (Fig. 3d). In accordance with previous reports, both astrocytes and neurons highly express Aβ, and GFAP(+) reactive astrocytes accumulate Aβ in the process of clearance in AD [39,40,41].
By mapping the masks of these gated objects back onto the coordinates of the original segmented images (Fig. 3e), we created a cell phenotype map (CPM) that illustrated a more even distribution of hippocampal microglia relative to astrocytes (Fig. 3e, f1, Additional file 1: Fig. S3A). Astrocytes were enriched in hippocampal white matter and vascular regions (Fig. 3e, f2, Additional file 1: Fig. S3B) as they are integral parts of CNS white matter and blood–brain-barrier (BBB) architecture [42, 43]. Vascular CPM charts large and micro vessel boundaries (Fig. 3e, f3, Additional file 1: Fig. S3C). In the neuron CPM, masks are enriched in the granule cell layer of DG band, and are scattered throughout the image (Fig. 3e, f4, Additional file 1: Fig. S3D). These data highlight the utility of the combined segmentation approaches, enabling analyses analogous to those achieved in other tissue where spatial organization, expression, and molecular identities of single cells and proteopathy can now be determined and organized for human CNS.
Top-down spatial organization of cellular and proteopathy composition distinguishes hippocampal subregions and disease status
We combined our discretization strategies with prior knowledge of human hippocampal anatomy (i.e., top-down) to identify neuropathologic changes across stages of cognitive decline (defined by Thal phase, Braak stage, CERAD and ADNC scores, Additional file 1: Table S2B). Coronal hippocampal sections from three individuals were imaged, focusing on the DG to CA4 through CA1 subregions, and expanding outward (Fig. 4). Samples from cognitively normal (CN), cognitively impaired with no dementia (CIND), and ADD subjects were used to capture a wide spectrum of pathologic changes (Fig. 4a). Both MCI and CIND are attempts to define a state in-between normal cognitive function and dementia. Mild Cognitive Impairment (MCI) is a clinical diagnosis [44]. Cognitive Impairment with no Dementia (CIND) is determined by neuropsychological test results [45, 46]. Single-plex spectral images of CN (Additional file 1: Fig. S4A–C), CIND (Additional file 1: Fig. S4D–F) and ADD (Additional file 1: Fig. S4G–I) hippocampus demonstrated larger, qualitative structures that can be captured by tiling numerous imaged fields together. For instance, VGLUT1 and VGLUT2 positivity demarking the CA and DG borders, respectively.
Given the known differences and progression of disease within AD hippocampal subregions [4], we first segregated all cells and protein aggregates into DG, CA4, CA3, CA2, and CA1 associations (Fig. 4b). We determined these boundaries for the subregions based on prior knowledge of morphological characteristics [47, 48] and expression of neuronal and pan-CNS markers including CALBINDIN, CALRETININ, MAP2, CD56, SYNAPTOPHYSIN (SYP), VGLUT1, and VGLUT2 (Fig. 4a, Additional file 1: Fig. S4) that have previously been shown to highlight boundaries [21, 24, 49]. Spectral images of CALBINDIN, CALRETININ, and MAP2, specifically showed anatomical delineation of CA4–CA2 (CALBINDIN), DG and Alveus (CALRETININ), joined by neuronal cytoarchitecture from CA4 to the subiculum (MAP2) (Fig. 4a, bottom panel). We then assigned cells and protein aggregates to their respective structures (Fig. 4b, Additional file 1: Table S5A). To assign cell identities, we applied a combination of FlowSOM [50] and manual meta-clustering strategies to parse out the objects into six distinct categories: neurons, microglia, endothelial cells, and non-immune glia which included both astrocytes and oligodendrocytes (Fig. 4b, lower panel). Astrocytes and oligodendrocytes were grouped together based on MBP, MAG, and GFAP expression, as their expression patterns overlapped enough at this lower resolution to make it difficult to assign individual labels for these subtypes.
Comparing the proportion of each cell type within the CA and DG subregions, we observed a higher proportion of protein aggregates in the CA1 subfield relative to other regions within the same sample, regardless of cognitive status (Fig. 4c, Additional file 1: Table S5B). In addition, the total number of protein aggregates increased in the tissue of cognitively impaired individuals (Additional file 1: Fig. S4J). We then considered whether proteopathies tended to lay in proximity to cell somas bound by our cellular segmentation. For each cell subtype, we counted the number of cells in direct proximity (i.e., overlapping pixels) to each proteopathy subtype and divided this value by the number of cells with no proteopathy overlap (Fig. 4d). Most cells did not exhibit direct proximity to Aβ plaque objects, particularly in the CN and CIND conditions. In the ADD individual, more cells lay in proximity to Aβ plaques, particularly for neurons in the CA1 and CA2 subfields, similar to previous reports [51]. For PHF1-TAU labeled NFT-NTs, a similar trend was found. In particular, microglia surrounding NFT-NTs doubled the number of microglia independent of NFT-NTs in the ADD CA1 subfield (Fig. 4d, green circle).
Building upon this observation, we contrasted those microglia that were NFT-NT associated (positive) or not (negative) (Fig. 4e). Microglia markers associated with reactivity [51], APOE, Iba1, CD33, and CD45, were consistently higher in NFT-NTs(+) relative to NFT-NTs(−) CA1 microglia, across all samples, with the highest expression in ADD (Fig. 4e). Microglia interacting with Aβ plaques did not show as strong of a difference in these reactive markers compared to microglia independent of Aβ plaques (Fig. 4f). Projecting an equal subset of microglia from all three samples onto an 2D-UMAP embedding, shows subsets of these reactive microglia associated with PHF1-TAU, particularly high Iba1(+) and APOE(+) in CA1 ADD cells (Fig. 4g, green circle); and that while these cells exist across all individuals, they are most prominent in ADD (Fig. 4g, h, Additional file 1: Fig. S4K–L). This differential state of expression further reflects the possibility that immune cell reactivity may play a central role in AD severity, particularly an association with PHF1-TAU formation in CA1 subregion [51].
Bottom-up, data driven neighborhood analysis identifies common regions of neuropathology across individuals and severity
Given the level of neuropathological organization identified using previously defined hippocampal subregions, we investigated what other levels of spatial order could be revealed with a more data-driven approach to our multiplexed images. In addition to including the DG and CA regions, we also included cells from additional areas of the hippocampus captured including the alveus and subiculum. To this end, we employed a bottom-up workflow to isolate common signatures of hippocampal spatial identities, independent of previously known neuroanatomy or cognitive status. Using CytoMAP [18], we grouped local neighborhoods of similar protein aggregates and cell types across the three samples. After neighborhoods were calculated, a minimum number of common regions into which they could be grouped into was determined (see Methods: Object co-proximity analysis). Identified cells and protein aggregates were then assigned to these regions across each sample. A voronoi expansion from cell and object bearing areas was then used to cluster the surrounding neuropil (i.e., cell projections and extracellular matrix) into one of the regional groups (Fig. 5a). The final five regions determined across all samples were annotated as neuronal dominant (N), glial dominant (G), mixed proteopathy (M), Aβ plaque dominant (AP), or PHF1-TAU NFT-NTs dominant (TT) (Fig. 5b). The three proteopathy associated regions: M, AP, and TT were increased in representation with increasing cognitive impairment, while the relatively protein aggregate-free neuronal (N) and glial (G) cell regions decreased with increasing cognitive impairment, making up less than a quarter of the total ADD tissue (Fig. 5c, Additional file 1: Table 5SC). Unsurprisingly, even the cognitively normal hippocampus contained AP and TT regions (Fig. 5b, left purple and blue), consistent with the high prevalence of pre-clinical AD in older people. The mixed proteopathy M region represents a mixture of high proportions of neurons, Aβ plaques, and tau NFT-NTs (Fig. 5b, orange, Additional file 1: Table 5SC). Interestingly, this mixed region grew predominantly in the CA4-CA2 subregion, known for representing a resilient area of hippocampal degeneration [4]. CA1, a degeneration susceptible area, instead contained much more of the tangle dominant TT region as disease status worsened. Therefore, the M region could potentially indicate a transitional state of proteopathy development that could be related to resilience. Moreover, due to the data-driven nature of their identification we saw similar compositions of cells and proteopathy within a given region and across samples, regardless of cognitive diagnosis (Fig. 5d).
To integrate these regions and pathologic changes into a single model, we applied UMAP [37] dimensionality reduction to all cells and protein aggregates, comparing how they overlapped in disease severity or data-driven regions (Fig. 5e). While like cells or protein aggregates clustered tightly, there was some divergence with de novo region, suggesting unique combinations of protein expression with both anatomical location and spatially associated pathologic features. For example, a subset of neurons, non-immune glia, and microglia clustered with the TT region, separate from those in the N or G regions (Fig. 5e, black circle). Given the influence of de novo identified regions on cell and object embedding (Fig. 5e), we also quantified how cell-protein aggregate proximity itself compared between the top-down anatomical (Fig. 4d) and data-driven regions (Fig. 5f). Regardless of sample staging, cells within the N- or G-regions showed very little proximity to pathologic objects, as expected. In the AP-region, nonimmune glia and neurons showed high proximity to neurons and nonimmune glia (Fig. 5f). The same trend held for PHF1-TAU NFT-NTs in proximity to cells in the TT-region. As pathologic changes increased, microglia once again stood out with the highest ratio of cells with NFT-NT proximity. These physical associations increased as severity worsened whether the nuclear associated cell body or protein aggregates were used as the basis for the proximity calculation (Additional file 1: Fig. S5A). In summary, while NFT-NTs development started in neuronal peripheries before accumulating in soma [38, 52, 53] our data suggests soma involvement continues as disease progresses.
With the connection between proteopathy formation and synaptic density, we further investigated the mean pixel expression of these markers across all identified data-driven regions in the hippocampus and how coincidence changed. Focusing on synaptic density, we calculated synaptic positive pixels by either excitatory (SYP(+), PSD95(+), and either VGLUT1(+) or VGLUT2(+)) or inhibitory (SYP(+), VGAT(+), GAD(+)) synaptic pixels. We then looked at the coincidence with PHF1-TAU and observed that while the highest percentage of PHF1-TAU formation in synapses occurred in ADD, all samples that contained the TT-region showed synaptic protein/PHF1-TAU localization (Fig. 5g), similar to what we have recently reported on synaptosomes [54]. Similar results were seen for Αβ plaque proximity to synaptic pixels in the AP region (Fig. 5g). This approach, as summarized in Fig. 5a, shows the ability to identify common pathological changes and neighborhoods across tissues representing a diversity of neuropathologic changes. Moreover, completely data-driven, bottom-up approaches, like the ones employed here, captured, and organized neuropathological features beyond the bounds of conventional anatomical subregions potentially increasing the sensitivity to detect and classify dysfunction.
Clustering identifies mitochondrial MFN2 protein expression in persistent neurons associated with proteopathy-laden tissue regions
Leveraging the 14 proteins simultaneously measured to capture neuron identity and putative function in situ, we clustered these extracted cells using FlowSOM to identify 11 meta clusters. The expression distributions across these groups (Fig. 6a) and as seen by UMAP projection (Fig. 6b, c) and physical mapping (Additional file 1: Fig. S6A) showed that neurons fell broadly into two categories: (1) proteopathy associated (clusters 1, 2, 6, 8, 10), and cluster (2) non-proteopathy associated (clusters 3, 4, 7, 9, 11). While this analysis highlighted the previous trend of increased protein expression in proteopathy-associated cells (i.e., synaptic markers, Fig. 6a, Additional file 1: Fig. S5B) it revealed some functional neuronal subsets based on mitochondrial proteins and ubiquitination, e.g., MITOFUSION 2 (MFN2) (Fig. 6a, clusters 1, 2, 3, 4, 10) and POLYUBIQUITIN K48 (Fig. 6a, clusters 1, 2, 8).
As with the broader analysis of all objects (Fig. 5), neuron subclusters associated with proteopathy grouped together in a UMAP of neurons alone (Fig. 6c, left). As MFN2 was variably expressed across multiple non-proteopathy associated (clusters 3, 4) and proteopathy-associated neuronal clusters (clusters 1, 2, 10), we further investigated its expression and localization in the hippocampus. Interestingly, recent single nucleus analysis of transposon accessible chromatin (ATAC-seq) indicates that MFN2 is epigenetically restricted to neuron-specific lineages in human brain (Additional file 1: Fig. S6B) [55], re-enforcing its potentially unique role in neuronal function here. Neuronal MFN2 expression was positive primarily in non-proteopathy associated neurons, with a clear overlap in proteopathy associated neurons. (Fig. 6d, right, black arrow). These MFN2 positive neurons were more prevalent in proteopathy-associated regions (MD, AP, TT) relative to the non-proteopathy associated regions (ND, GD) (Fig. 6d), suggesting that MFN2 expression may be associated with neuronal survival in the presence of stressors from proteopathy.
To understand MFN2-related persistence in proteopathy-laden regions, we compared the single-cell relationship between proteopathy expression and MFN2. Proteopathy- associated regions have higher proportions of MFN2(+) neurons relative to non-proteopathy associated regions in the CN hippocampus (Fig. 6f, Additional file 1: Table S6A). Additionally, in the ADD sample the percentage of MFN2(+) neurons in the non-proteopathy associated regions has jumped up considerably relative to the CN sample, further suggesting that this MFN2(+) expression may indicate an association with protection from the spread of proteopathy. Filtering by de novo region shows high abundance of MFN2(+), PHF1-TAU(+) neurons in the TT region (~ 20%) compared to the N region (~ 1%) (Fig. 6g, bottom, Fig. 6S6C; Additional file 1: Table S6B–C). These PHF1-TAU(+), MFN2(+) neurons can be seen throughout the TT region of the hippocampus, even appearing adjacent to PHF1-TAU single positive, as well as MFN2 single positive expressing neurons, regardless of disease status (Fig. 6h). Comparatively, the N region contained many more MFN2(+) single-positive neurons. This persistence of MFN2(+) neurons in areas of high tauopathy potentially reflects a neuronal protective response to injury as has been studied in animal and in vitro models of neuronal degradation [56]. MFN2’s relationship with Aβ42 did not show this same L-shaped distribution, especially in the AP region (Fig. 6g, top), suggesting that the presumed MFN2 response to injury may be specific to stressors that accompany pathologic Tau accumulation. We validated these expression patterns with IHC staining of PHF1-TAU, Aβ42, MFN2 on each of the tissue samples (Additional file 1: Fig. S6D–F). Altogether, our study demonstrates both anatomical and data-driven approaches can be used to identify deep spatial summaries of neuropathologic changes while identifying and subclassifying cells and features that compose them.