The landscape of myeloid and astrocyte phenotypes in acute multiple sclerosis lesions

Abstract

Activated myeloid cells and astrocytes are the predominant cell types in active multiple sclerosis (MS) lesions. Both cell types can adopt diverse functional states that play critical roles in lesion formation and resolution. In order to identify phenotypic subsets of myeloid cells and astrocytes, we profiled two active MS lesions with thirteen glial activation markers using imaging mass cytometry (IMC), a method for multiplexed labeling of histological sections. In the acutely demyelinating lesion, we found multiple distinct myeloid and astrocyte phenotypes that populated separate lesion zones. In the post-demyelinating lesion, phenotypes were less distinct and more uniformly distributed. In both lesions cell-to-cell interactions were not random, but occurred between specific glial subpopulations and lymphocytes. Finally, we demonstrated that myeloid, but not astrocyte phenotypes were activated along a lesion rim-to-center gradient, and that marker expression in glial cells at the lesion rim was driven more by cell-extrinsic factors than in cells at the center. This proof-of-concept study demonstrates that highly multiplexed tissue imaging, combined with the appropriate computational tools, is a powerful approach to study heterogeneity, spatial distribution and cellular interactions in the context of MS lesions. Identifying glial phenotypes and their interactions at different lesion stages may provide novel therapeutic targets for inhibiting acute demyelination and low-grade, chronic inflammation.

Introduction

Multiple sclerosis (MS) is a common neurological disease, characterized by the formation of inflammatory demyelinating lesions in the central nervous system (CNS) [26]. Inflammation is driven by infiltrating lymphocytes and monocytes, in concert with resident activated microglia and astrocytes. Macrophages and reactive astrocytes are the most abundant cell types in acute lesions [18, 30]. These cells are highly plastic and can adopt pro-inflammatory, anti-inflammatory, neurotoxic, neuroprotective, and tissue-regenerating functions [4, 6, 7, 20, 23, 29, 30, 42]. Previous studies have identified macrophage phenotypes in MS lesions based on the expression of single classical (M1) and alternative (M2) activation markers; however, those studies have produced limited and sometimes conflicting results [6, 42]. It is now increasingly clear that the M1/M2 polarization paradigm, which originated as an in vitro concept, is of limited value for distinguishing myeloid phenotypes in inflamed tissue [35]. Recent studies, including one of our own, have used single-cell or single-nucleus RNA sequencing on CNS tissue to comprehensively assess the complex phenotypes of human glial cells in healthy and diseased brains [14, 22, 27]. Similarly, in these studies, myeloid cell/microglial phenotypes did not separate into categories in which M1 and M2 markers are of organizational value.

Several novel histological techniques now make it possible to perform high parameter imaging of tissue sections and to evaluate complex cellular phenotypes in situ [5, 8, 9, 11, 41]. We have used imaging mass cytometry (IMC), a technique that like mass cytometry (CyTOF) relies on metal isotope-labeled antibodies, and combines immunohistochemistry with high-resolution laser ablation followed by time-of-flight mass spectrometry [9, 43]. This approach allows for simultaneous quantitative profiling with up to 37 antibodies on a single tissue section at subcellular resolution. Moreover, computational tools have become available to extract single-cell information from highly multiplexed histological data [3, 24, 36]. In this proof-of-concept study, we applied IMC and single-cell analytics to two active MS lesions – one demyelinating and one post-demyelinating – to examine the cellular heterogeneity of myeloid cells and astrocytes based on thirteen markers known to be expressed by activated glial cells in MS lesions [2, 6, 13, 15, 31, 42, 44, 45]. We demonstrate that multiplexed tissue imaging, in combination with the appropriate computational tools, can extract previously unattainable information from histological sections, including definition of cellular subpopulations, their distribution within the lesion environment, specific cell-cell interactions, phenotypic transitions and the impact of spatial sources on marker expression.

Materials and methods

MS lesions

Human CNS tissue was obtained at autopsy from two patients with relapsing-remitting MS according to institutional review board-approved protocols. After autopsy, brain tissue was fixed in 10% formalin, and lesions were cut based on MRI. Lesion tissue was subsequently embedded in paraffin and sectioned at 5 μm thickness.

A highly inflamed active lesion was chosen for analysis from each patient: the demyelinating lesion was selected from a 42-year-old male with 5 years disease duration (5.5 h post mortem interval), while the post-demyelinating lesion was chosen from a 32-year-old female with 6 years disease duration (8 h post mortem interval). Lesions from both patients have been characterized in previous studies [10, 28].

Brightfield histology

For basic characterization, lesions were stained against CD68, myelin basic protein (MBP) and MAP2, and examined via brightfield microscopy. De-paraffinized and rehydrated sections were subjected to antigen retrieval in pH 6, 10 mM citrate buffer at 96 °C for 20 min, cooled, quenched in 0.3% peroxide, and blocked with FC receptor binding inhibitor and normal serum before incubation with primary antibody (Additional file 1: Table S1) overnight at 4 °C. Sections were subsequently incubated with appropriate biotinylated secondary antibodies, processed with an avidin/biotin staining kit with 3,3-diaminobenzidene (DAB) as the chromogen (Vector ABC Elite Kit and DAB Kit, Vector Laboratories), then counterstained with hematoxylin [28]. Adequate controls using isotype control antibodies were performed for each primary antibody. Sections were rinsed with distilled water, dehydrated, and cover-slipped with Permount (Vector Laboratories). Images were acquired using a Leica DM5000 B microscope with a Leica color camera DFC310 Fx and Leica Application Suite (version 4.9.0) imaging software. Images were processed with Panoramic Viewer (3DHISTECH) and Photoshop (Adobe) software.

Antibody validation and conjugation to metal isotopes for IMC

Lanthanide-conjugated antibodies were purchased from Fluidigm. Antibodies not available in metal-conjugated form were purchased in carrier-free solution and validated by brightfield immunohistochemistry using the appropriate isotype control antibodies. Subsequently, antibodies were conjugated to lanthanide metal isotopes following the Maxpar® Antibody Labeling Kit protocol (Fluidigm). Briefly, carrier-free antibodies were partially reduced with Bond-Breaker™ TCEP buffer (Thermo Scientific) at 37 °C before incubation with purified, lanthanide-loaded Maxpar® X8 polymer at 37 °C for 90 min. The percent yield of metal-conjugated antibodies was determined using the Pierce™ BCA Protein Assay Kit (Thermo Scientific). Metal-conjugated antibodies were stored at 0.5 mg/mL in PBS-based Antibody Stabilizer (Candor Bioscience) with 0.05% sodium azide at 4 °C. Working concentrations for all metal-conjugated antibodies were optimized by IMC (Additional file 1: Table S2) on MS lesion tissue.

Imaging mass cytometry

For IMC histology, tissue sections were de-paraffinized and rehydrated, and antigen retrieval was performed in pH 8, 1 mM EDTA buffer at 96 °C for 20 min. Sections were cooled at room temperature and rinsed in tap water and TBS (20 mM Tris with 500 mM NaCl, pH 7.5). Tissue was blocked for 3 h at room temperature with 0.3% BSA, 25% FBS and 0.1 mg/mL FC receptor binding inhibitor in TBS-T (TBS + 0.05% Tween-20). All antibodies (Additional file 1: Table S2) were diluted in 0.3% BSA in TBS-T and applied to the tissue for overnight incubation at 4 °C. Sections were then rinsed in TBS-T and TBS, and counterstained with 125 nM Maxpar® Intercalator-Ir (Fluidigm) in PBS for 1 h at room temperature. Sections were rinsed in TBS-T, TBS, and two washes of distilled water before air-drying at room temperature. Antibody-labeled tissue areas (1000 × 1000 μm) were raster-ablated using a Hyperion™ Laser Scanning Module (Fluidigm) with a 1 μm diameter spot size at 200 Hz. This process was coupled to a Helios™ Mass Cytometer (Fluidigm) for lanthanide metal detection [43]. Images for each antibody channel were acquired on CyTOF Software (Fluidigm, version 6.7). MCD Viewer (Fluidigm, version 1.0) was used to export raw 16-bit tiff images for computational analyses on histoCAT (version 1.75) [36]. For visualization purposes, images were processed in MCD Viewer and ImageJ [38].

Computational analyses

Single-cell segmentation

Merged CD68 (macrophages/microglia), S100B (astrocytes), and CD3 (T cells) antibody channel images were first processed on MCD Viewer and ImageJ to diminish unspecific staining noise that could interfere with segmentation. These image adjustments helped distinguish closely neighboring CD68, S100B and CD3 cell bodies from one another. The resulting images served as the staining template to segment single cell objects on CellProfiler (version 3.0.0) [16]. The IdentifyPrimaryObjects module was used for segmentation with three-class Adaptive Otsu thresholding, shape or signal intensity-based declumping, and the propagation method for drawing dividing lines between clumped cell objects. The typical object diameter was assigned based on the approximate range of cell sizes present in an image. For three-class thresholding, the middle class was assigned to either foreground or background based on the intensity of residual noise in an image. These optimal parameters were determined based on the following criteria: each segmented cell had one nucleus associated with it, the complex morphology of CD68⁺ myeloid cells and S100B⁺ astrocytes were reflected in the cell outlines, and co-segmentation of the different cell type markers was minimized to the highest extent possible. Fulfillment of these requirements was checked by visualizing the segmentation masks over merged CD68, S100B, CD3 and nuclear counterstain images on histoCAT. Moreover, perivascular CD68⁺ and CD3⁺ cells in the early lesion were too densely packed to separate them by segmentation, and were eliminated in CellProfiler with the EditObjectsManually module. The resulting segmentation mask images with outlined cell borders were exported from CellProfiler as 16-bit unsigned integer (uint16) images and loaded into histoCAT with corresponding IMC antibody channel images.

Identification of cellular phenotypes

On histoCAT, mean single-cell marker intensity values were extracted via segmentation masks from raw, 16-bit tiff images for each antibody channel and Z-score normalized per marker. Based on the expression intensities of thirteen markers (Additional file 1: Table S2), cell clusters were defined using the PhenoGraph algorithm [19] integrated into histoCAT. Default parameters with 75 nearest neighbors for the early lesion and 50 nearest neighbors for the late lesion were used. These nearest neighbor values were chosen such that over- and under-clustering of phenotypes were avoided. Additional normalization steps were performed internally, as previously described [36].

Analysis of cellular phenotypes

To visualize clusters, the Barnes-Hut t-SNE algorithm implemented in histoCAT was executed with the same image and marker inputs used in PhenoGraph, as well as default parameters (initial dimensions, 110; perplexity, 30; theta, 0.5) and internal normalization [1, 36]. t-SNE plots were colored to highlight cell clusters or lesion samples, or to show relative marker expression intensity. Images of cell phenotypes visualized in the tissue, as well as segmentation masks overlaid with histology images, were generated in histoCAT. For the remaining analyses, “.csv” files containing single-cell parameters were exported from histoCAT and appropriately processed for their application. To produce an expression heatmap for clusters, Z-score normalized marker intensity values were processed using the R ComplexHeatmap package, which hierarchically clusters single cells within clusters using Ward’s method [37]. Violin plots showing single-cell marker expression variability for each cluster were generated using the R ggplot2 package [12].

To study phenotype transitions, Potential of Heat-diffusion Affinity-based Transition Embedding (PHATE) mapping and Monocle 2 Pseudotime analyses were performed in R [24, 33, 34, 40]. For these analyses, files containing single-cell marker expression values along with metadata files indexing the phenotype of each cell served as input, as detailed in the online user guides. Moreover, relevant markers were selected for analyses of myeloid cells (all markers except for CD3, S100B and vimentin) or astrocytes (all markers except for CD3 and CD68).

PHATE mapping was performed using Z-score normalized marker intensity values, to be consistent with the input used in histoCAT for cell clustering and t-SNE plots. With this input, we tested different values of the adjustable, nearest neighbor parameter k that is built into the PHATE algorithm. Larger values of k make transition states less distinct from one another, while smaller values increase the influence of any artefacts on the analysis [19]. Since k is a nearest neighbor metric and large values could obscure biologically meaningful transitions, we used values that were significantly less than the number of cells analyzed, but large enough to avoid misleading results due to noise. We found that different values of k above 30 did not significantly change the results and ultimately chose k = 100 for our analyses. Other parameters were left as their default specifications.

Monocle 2 Pseudotime analysis was performed with internal, negative binomial normalization of raw marker intensity values, as suggested by the user guide. A test for differential marker expression among phenotypes was performed as part of the analysis, generating q values for each marker that indicated how significant expression differences were among phenotypes. Provided that Pseudotime was developed to accommodate full-transcriptome datasets, this test would allow genes with similar expression levels across phenotypes to be filtered out based on a chosen q value threshold. We broadly set q < 1 as the criteria for using a marker in the ordering_genes function after reviewing the test results, to ensure that all parameters would be used as they were in PHATE mapping.

Figures showing phenotype cell size and quantity as well as correlation matrices and plots were produced in Prism (version 7). FlowJo software (version 10.5.3) was used to visualize single-cell marker data on flow cytometry plots. Images and figures were recolored when necessary in Photoshop.

Analysis of cell spatial relationships

To study the spatial relationships of cell clusters, neighborhood analysis was performed on histoCAT using the PhenoGraph-generated clusters pertaining to each lesion. Significant pairwise phenotype interactions and avoidances were determined by an unbiased permutation test which compared the frequency of one cell type neighboring another to the frequency in 999 random permutations of cell cluster labels. Neighbors were identified within 4 μm of each cell, and neighboring frequency was normalized to the number of interacting cells [36]. Analyses were performed for different degrees of significance (p < 0.05 and p < 0.01), and the results were reported on a heatmap. To identify sources of single-cell marker variation, spatial variance component analysis (SVCA) [3] was performed on Python for different lesion zones, using standardized marker intensity values and the spatial coordinates of each cell in tissue images. SVCA plots were generated in R and Prism.

Statistical analyses

In plots showing phenotype cell size, data represent mean cell areas + standard deviation. Comparisons of phenotype cell sizes were analyzed by one-way ANOVA followed by the Tukey-Kramer multiple comparison test. Comparisons of two samples were carried out by unpaired Student’s t-tests. For correlation analyses, Pearson correlation coefficients were computed. *p < 0.0001.

Results

Histology and cell clustering overview

We analyzed one demyelinating and one post-demyelinating MS lesion according to a classification by Kuhlmann et al. [17], referred to throughout this report as the “early” and “late” lesion, respectively. Both lesions were located in the brain stem and characterized by complete loss of myelin, hypercellularity with the highest cellular density at the lesion rim, and diffuse infiltration with foamy macrophages. The laser-scanned area of the early lesion involved predominantly white matter (WM), but also interspersed gray matter (G/WM), while the scanned area of the late lesion consisted of only WM (Additional file 1: Figure S1). Consistent with demyelinating activity, macrophages at the rim of the early lesion contained myelin basic protein (MBP)-positive myelin debris, which was absent in macrophages from the late lesion (Additional file 1: Figure S1d, j) [17]. Moreover, foamy macrophages were numerous and large in size in the rim of the early lesion, whereas macrophages in the rim of the late lesion were smaller and less dense (Additional file 1: Figure S1e, k). Perivascular infiltrates in the early lesion contained mostly lymphocytes and only a few undifferentiated monocytes, while the perivascular cuffs in the late lesion consisted predominantly of lipid-laden macrophages, as described previously (Additional file 1: Figure S1i, l) [21, 39].

We immunolabeled both lesions with antibodies against the cellular and activation markers CD68, S100B, CD3, PLP, CD45, CD86, ferritin, HLA-DR, LAMP1, Mac2, MerTK, TIM-3 and vimentin; all known to be expressed in MS lesions [2, 6, 13, 15, 31, 42, 44, 45]. We defined single cells by segmenting cell bodies outlined by the markers CD68 (macrophages/microglia), S100B (astrocytes), and CD3 (T cells) using CellProfiler (Fig. 1; see Methods for a detailed description of segmentation parameters) [16]. Our segmentation methods captured the complex morphologies of myeloid cells and astrocytes, and allowed for clear delineation of cell types in the vast majority of segmented cells. In a small fraction of these cells, overlap between cellular markers could not be avoided with our segmentation pipeline. This overlap occurred between myelin fibers and juxtaposed microglial processes in normal-appearing white matter (NAWM), and between astrocyte processes and closely neighboring macrophages (Fig. 1c, d). We did not exclude these cells in order to avoid bias in our spatial analysis. However, we excluded perivascular lymphocytes from the early lesion, as they were too densely packed to be identifiable as individual cells.

Fig. 1

Single-cell segmentation. (a-c) Segmentation of myeloid cells (CD68, magenta), astrocytes (S100B, green) and T cells (CD3, cyan) performed on CellProfiler. (b) shows cells in (a) at higher magnification. (c) shows examples of S100B⁺ astrocyte processes (white arrows) that closely neighbor CD68⁺ cells, which resulted in the co-segmentation of these markers in a small fraction of cells. (d) Gating of CD68⁺ and S100B⁺ populations on a flow cytometry plot as a quality control for segmentation (late lesion). The overlap population consists of cells with co-segmented CD68 and S100B like those in (c). Scale bar in (a) = 30 μm; (b) = 15 μm; (c) = 30 μm, inset = 10 μm. Z-score normalized expression intensities are shown (d)

Single-cell segmentation masks were overlaid with IMC images from all thirteen antibody channels (Additional file 1: Figures S2, S3) in histoCAT [36] to extract mean single-cell marker expression intensity values from the images, and to cluster myeloid cells and astrocytes into phenotypic subpopulations based on normalized expression intensities. We obtained the optimal number of phenotype clusters, based on separation of the principal cell types (i.e. myeloid cells, astrocytes and T cells) and distinctive expression profiles of phenotypes on expression heatmaps and t-SNE plots, resulting in a total of twelve phenotypes in each lesion (Figs. 2, 4). Given differences in their cellular composition, this required different nearest neighbor values for the early and late lesion (75 and 50 nearest neighbors, respectively). t-SNE plots were generated in histoCAT with the same image and normalized marker expression inputs were used to generate phenotype clusters and phenotype locations in the tissue were automatically labeled in images following clustering.

Fig. 2

Early lesion cell phenotype profiles. (a) Marker expression heatmap for the myeloid, astrocyte and T cell phenotypes, identified by PhenoGraph clustering on histoCAT using segmented cells (n = 4397). The heatmap displays relative expression levels based on Z-score normalized marker intensity values, and single cells are hierarchically clustered within each phenotype group. Labels at the bottom of the heatmap indicate the area of the lesion to which each phenotype localizes. (b) Violin plot representation of the data in (a). (c) t-SNE plot showing the distinct phenotype clusters. (d) t-SNE plot colored by marker intensity, confirming the separation of CD68⁺, S100B⁺ and CD3⁺ cell types. G/WMoR = gray and white matter outer rim; WMoR = white matter outer rim; iR = inner rim; WM R/C = white matter rim/center; C = center; NAWM = normal-appearing white matter

Phenotype heterogeneity and distribution in early demyelination

In the early lesion, we analyzed a total of 4397 cells, of which 66.3% were myeloid cells and 32.5% were astrocytes (Additional file 1: Figure S4a). This ratio was higher at the rim than the lesion center. The cells clustered into six myeloid and five astrocyte subtypes (Fig. 2). The myeloid phenotypes were spatially sequestered into four lesional regions, i.e. lesion center (phenotypes m2 and 4), inner lesion rim (m3), outer lesion rim (m1 and 5) and NAWM (m6) (Fig. 3a). Most myeloid cell activation markers were highly expressed in the outer lesion rim and decreased in intensity towards the lesion center (Fig. 2a, b). The m1 phenotype at the G/WM interface showed the highest activation profile. m5 cells in the outer WM rim were the largest in size, whereas m2 cells in the lesion center were the smallest (Additional file 1: Figure S4b), which may reflect the amount of myelin fragments phagocytosed by foamy macrophages at the advancing lesion edge, and the degradation of myelin in phagocytes within the lesion center [32].

Fig. 3

Early lesion phenotype spatial distributions and transition analyses. (a, b) Spatial separation of (a) myeloid cell and (b) astrocyte phenotypes into NAWM, rim and center lesion zones. (c) T cells are primarily located in the lesion center. (d) PHATE mapping of myeloid cells, indicating that the G/WM outer rim (m1) and lesion center phenotypes (m2) are on a transition continuum (white arrow). (e) Pseudotime analysis of myeloid cells shows that phenotypes transition along independent trajectories. Phenotype color schemes on the PHATE and Pseudotime plots reflect the color palettes specific to each analysis. Scale bar for a-c = 200 μm. G/WMoR = gray and white matter outer rim; WMoR = white matter outer rim; iR = inner rim; WM R/C = white matter rim/center; C = center; NAWM = normal-appearing white matter

Astrocyte phenotypes were each defined by one distinct, highly expressed marker (Fig. 2a, b). Analogous to the myeloid phenotypes, the five astrocyte phenotypes were also spatially stratified into the NAWM, lesion rim, and lesion center (Fig. 3b). Furthermore, a5 astrocytes located in the outer rim within WM were larger than all other phenotypes (Additional file 1: Figure S4b). Unlike in myeloid cells, marker expression in astrocyte phenotypes did not follow a gradient from the rim to the lesion center, but was uniform throughout the lesion. T cells made up the smallest population of all immune infiltrates, and were concentrated within the lesion center (Fig. 3c; Additional file 1: Figure S4a). These cells uniformly expressed CD45 and HLA-DR (Fig. 2a, b), and did not separate into different clusters.

To determine possible transitions between phenotypes, we applied Potential of Heat-diffusion Affinity-based Transition Embedding (PHATE) mapping and Monocle 2 Pseudotime (referred to as Pseudotime) to the myeloid cell and astrocyte populations. PHATE mapping improves on t-SNE by visualizing phenotype transitions based on differential marker expression, where smooth continua from one phenotype to another indicate a transition trajectory [24]. Pseudotime creates a trajectory graph by computing a minimum spanning tree onto which the cells are projected [24, 40]. PHATE mapping of myeloid cells showed a linear transition continuum from the G/WM to the WM outer rim phenotypes (m1 and m5) and the lesion center phenotype (m2). (Fig. 3d). This transition did not include the inner rim phenotype (m3). Pseudotemporal ordering of myeloid cells with Pseudotime did not result in a linear transition trajectory of myeloid phenotypes but suggested several independent fates (Fig. 3e). Similarly, PHATE mapping and Pseudotime analysis of astrocyte phenotypes indicated independent phenotypic fates that did not transition into one another (Additional file 1: Figure S5a, b).

Low phenotype heterogeneity and random phenotype distributions in the late demyelinating lesion

In the late, post-demyelinating lesion we analyzed 6698 cells, with myeloid cells far outnumbering astrocytes (91.1% myeloid cells; Additional file 1: Figure S6a), particularly at the lesion rim. The same clustering criteria used for the early lesion resulted in nine myeloid phenotypes and two astrocyte phenotypes in this late lesion (Fig. 4). Myeloid phenotypes separated into lesion rim (m3), perivascular space (m4), and NAWM (m7) zones (Fig. 5a). In contrast to the early lesion, the six other myeloid phenotypes were intermixed throughout the lesion center. These phenotypes showed a low degree of separation on the t-SNE plot, indicating similar marker expression profiles (Fig. 4c). The phenotypes in the lesion rim and the perivascular space (m3, m4) were characterized by high expression of the majority of markers, and shared a similar expression profile with the G/WM rim phenotype in the early lesion (m1) (Fig. 4a, b). As in the early lesion, myeloid phenotypes in the lesion rim (m3) and perivascular space (m4) were significantly larger than those in the lesion center (Additional file 1: Figure S6b), but were overall smaller than in the early lesion (Additional file 1: Fig. S6c). Astrocytes were clustered into two phenotypes, with one phenotype localizing primarily to the lesion rim, and the other to the lesion center (Fig. 5b). The rim phenotype (a2) exhibited a marker expression profile similar to the rim phenotype in the early lesion (a5), (Fig. 4a, b). Finally, as in the early lesions, T cells were few (Fig. 5c; Additional file 1: Figure S6a), and expressed the activation markers CD45 and HLA-DR (Fig. 4a, b). To directly compare cell populations in both lesions, we mapped cells from both lesions on the same t-SNE plot. Cell populations overlapped moderately, highlighting differences between the phenotypes in each lesion (Additional file 1: Figure S7).

Fig. 4

Late lesion cell phenotype profiles. (a) Marker expression heatmap for the myeloid, astrocyte and T cell phenotypes, identified by PhenoGraph clustering on histoCAT using segmented cells (n = 6698). The heatmap displays relative expression levels based on Z-score normalized marker intensity values, and single cells are hierarchically clustered within each phenotype group. Labels at the bottom of the heatmap indicate the area of the lesion to which each phenotype localizes. (b) Violin plot representation of the data in (a). (c) t-SNE plot showing the phenotype clusters. Compared to the early lesion, myeloid cell phenotypes show a low degree of separation. (d) t-SNE plot colored by marker intensity, confirming the separation of CD68⁺, S100B⁺ and CD3⁺ cell types. R = rim; C = center; PVS = perivascular space; NAWM = normal-appearing white matter

Fig. 5

Late lesion phenotype spatial distributions and transition analyses. (a) Spatial organization of myeloid cell phenotypes in the lesion. Phenotypes in the rim (m3), perivascular space (m4) and NAWM (m7) separate into distinct zones, while lesion center phenotypes (m1, m2, m5, m6, m8, m9) are uniformly distributed. (b) Spatial distribution of astrocyte phenotypes. One phenotype (a1) predominantly occupies the lesion center and the other (a2) occupies the rim. (c) T cells are primarily distributed in the lesion center. (d) PHATE mapping of myeloid cells, showing that two lesion center phenotypes (m2, m6) are on a continuum with perivascular space cells (m4, white arrows). (e) Pseudotime analysis of myeloid cells shows a similar trajectory to PHATE mapping. Phenotype color schemes on the PHATE and Pseudotime plots reflect the color palettes specific to each analysis. Scale bar for a-c = 200 μm. PVS = perivascular space; NAWM = normal-appearing white matter

PHATE mapping demonstrated a linear transition continuum between two lesion center myeloid phenotypes (m2, m6) and the perivascular space phenotype (m4) (Fig. 5d), which was confirmed by Pseudotime (Fig. 5e), supporting a lesion center-to-perivascular phenotype trajectory, but not a continuum where all phenotypes align along a rim-to-center transition axis. The same analyses for astrocytes showed overlap between both phenotypes but not a linear transition (Additional file 1: Figure S5c, d).

Finally, in the early lesion, we found no correlation between the expression intensities of different markers at a single-cell level (Additional file 1: Figure S8a, b). In the late lesion, we found strong correlations between the M2 markers MerTK and Mac2, and MerTK and LAMP1 in both myeloid cells and astrocytes (Additional file 1: Figure S8c, d), resulting from high and continuous dynamic ranges of marker expression (Fig. 6).

Fig. 6

Single-cell marker correlations in the late lesion. (a, b) Co-expression of Mac2 and MerTK, and LAMP1 and MerTK in (a) myeloid cells (n = 6100), and (b) astrocytes (n = 528). Co-expression plots are shown with linear expression values and regression lines with Pearson correlation coefficients, and in flow cytometry contour plot form with log10-transformed expression values

Phenotypes in the early and late acute lesions engaged in specific cell-cell interactions

We next investigated the spatial relationships between different phenotypes with a computational tool integrated into histoCAT that performs an unbiased, systematic analysis of pairwise phenotype interactions and avoidances [36]. After excluding interactions between cells of the same or spatially adjacent phenotypes, our analysis demonstrated distinct interaction signatures for both lesions, (Fig. 7a, b). At a significance cut-off of p < 0.01, these included interactions between inner rim myeloid phenotype m3 (MerTK and CD45 high) and astrocyte phenotype a1 (Mac2 high), as well as interaction of highly activated myeloid rim phenotype m1 and the rim/center phenotype m4 (HLA-DR and ferritin high) with astrocyte phenotype a2 (HLA-DR high) in the early lesion. In the late lesion, highly activated perivascular macrophages (m4) interacted with most myeloid cell phenotypes and both astrocyte phenotypes. There were also significant interactions between myeloid and astrocyte phenotypes m6 and a1, and among lesion center myeloid phenotypes (m6 with m7 and m8). At a significance cut-off of p < 0.05, we found that T cells in the late lesion interacted with HLA-DR-expressing myeloid phenotypes in the perivascular space (m4) and in the lesion center (m8).

Fig. 7

Neighborhood and spatial variance component analyses for the early and late lesions. (a, b) Neighborhood analysis heatmaps of all significant pairwise phenotype interactions (red) and avoidances (blue) in the (a) early and (b) late lesions. White represents no significant spatial relationship. Dark boxes are highly significant spatial relationships (p < 0.01). Lightly shaded boxes are less significant relationships (p < 0.05) and interactions between cells of the same or spatially adjacent phenotypes. Rows visualize the significance of a phenotype surrounded by other phenotypes, and columns visualize the significance of a phenotype surrounding other phenotypes. (c, d) Spatial variance component analysis (SVCA) for the (c) early and (d) late lesions, showing the proportion of marker expression variance attributable to neighboring cell influences, environmental stimuli, cell-intrinsic factors and residual noise in different lesion zones. Additional plots show the average proportion of marker variance attributable to each factor in different lesion zones. G/WM = gray and white matter; WM = white matter; PVS = perivascular space; NAWM = normal-appearing white matter

Influence of lesion environment on marker expression

Finally, we used spatial variance component analysis (SVCA) to model the effects of extrinsic factors (neighboring cells and unobserved, non-cellular environmental stimuli) and cell-intrinsic factors on variations in cell marker expression, irrespective of phenotype [3]. This analysis was performed for each marker, using standardized single-cell marker expression values as well as the coordinate location of each segmented cell in the tissue. In the rim of both lesions, the expression of several markers was highly influenced by neighboring cells, including CD86, PLP, and Mac2 in the early lesion, and CD86, PLP and LAMP1 in the late lesion. Other markers, such as ferritin and vimentin (early lesion), and ferritin, vimentin, and TIM-3 (late lesion), were driven predominantly by non-cellular environmental stimuli (Fig. 7c, d). The relative influence of these factors changed toward the inner lesion rim and center, leading to an overall increased impact of cell-intrinsic factors and a decrease in the influence of external factors. In the lesion center, the primary agents influencing marker expression were cell-intrinsic factors, and to a lower degree environmental stimuli, while neighboring cells exerted no influence (Fig. 7c, d).

Discussion

Our study examines the landscape of myeloid and astrocyte phenotypes in early and late acute MS brain lesions using IMC. To our knowledge, this is the first application of highly multiplexed imaging to MS tissue. We applied thirteen markers that are known to be expressed by activated glial cells during MS lesion development. Clustering resulted in eleven myeloid cell and astrocyte phenotypes that localized to distinct lesion areas. Moreover, individual phenotypes interacted selectively with other cell types, and marker expression was driven by different factors in cells located at the lesion rim compared to the center. Thus, our approach provides a wealth of data on cellular spatial organization that is not accessible with standard histology.

The alignment of myeloid cell phenotypes with different lesional layers suggests functional specificity and validates our clustering approach. This spatial separation was most pronounced in the early lesion, and was reduced in the center of the late lesion where multiple phenotypes were intermixed. In addition, marker expression was the highest in myeloid phenotypes located at the lesion rim and diminished substantially towards the lesion center in both lesions. Consistent with the different stages of myelin phagocytosis and degradation, the myeloid phenotypes in the rim were larger than those in the lesion center. An additional feature of the late lesion was the presence of numerous highly activated macrophages in perivascular spaces throughout the lesion. Since these macrophages are believed to transition into the vasculature [21], this may indicate that they exit the CNS in a highly activated state. In contrast to myeloid cells, marker expression in astrocyte phenotypes did not follow a rim-to-center gradient, but was consistent throughout the lesion.

Our findings argue that macrophages/microglia in MS lesions do not transition from a pro- to anti-inflammatory state, as previously suggested [6], but convert from a highly activated to a less activated state as they move from the active edge to the lesion center. This is consistent with immunohistological results by Vogel and colleagues demonstrating that pro- and anti-inflammatory markers were simultaneously expressed by macrophages/microglia in MS lesions [42], and with single nucleus/cell RNA sequencing data of microglial cells in MS and neurodegenerative diseases, which do not produce categories related to M1 or M2 marker expression [22, 27]. Thus, our results add to the increasing evidence that activated macrophages and microglia in inflamed tissue do not follow a M1/M2 polarization dichotomy.

Using PHATE mapping, we found that myeloid cell but not astrocyte phenotypes followed a linear transition continuum from the G/WM outer rim to the WM outer rim and the lesion center (early lesion), and from lesion center phenotypes to the perivascular phenotype (late lesion). In contrast, phenotype trajectories on Monocle 2 Pseudotime showed no definite transition patterns. Although PHATE and Pseudotime provide biologically accurate transitions when applied to data sets with comparable parametric depth as ours, both methods have previously been shown to produce discrepant results, which may be attributed to their different computational approaches [24]. Our results deviate from the predicted transition of myeloid phenotypes from the outer to the inner rim and lesion center. Based on the myeloid states defined by our marker panel, myeloid cells develop along several independent fates, rather than one linear phenotype trajectory. We can however, not exclude that inclusion of more or different activation markers may yield different results.

The neighborhood analysis demonstrated distinct cellular interaction signatures for both lesions, e.g. between phagocytic inner rim macrophages and center astrocytes in the early lesion, and between T cells and two myeloid phenotypes in the late lesion. This indicates that cellular interactions in this hypercellular lesion environment are not random, but occur between specific subpopulations and cell types such as lymphocytes. The low parametric depth of our study does not allow us to identify the functional implications of these interactions; however, they may represent nodal points of cellular communication critical for lesion formation and maintenance of low-grade inflammation.

Finally, spatial variance component analysis (SVCA) suggests that cell-extrinsic factors drive marker expression to a higher degree in the lesion rim than in the center. Conversely, cell-intrinsic factors have a more prominent influence on marker expression in the lesion center. This suggests that glia cells in the lesion rim respond to cues from the microenvironment, such as cytokines or receptor-ligand interactions, while glial activation in the lesion center is the result of cell-intrinsic programs set in motion e.g. by myelin phagocytosis.

Myeloid cell/microglial heterogeneity has recently been examined by us and others with single-cell RNA sequencing in the healthy CNS, MS lesions, and other neurological diseases such as Alzheimer’s disease, Parkinson’s disease and temporal lobe epilepsy [22, 27]. These efforts have identified multiple myeloid cell/microglial phenotypes, comparable with our results. One of the microglia clusters, which was enriched for genes associated with MS susceptibility and characterized by high expression of CD74, was also enriched for genes that were highly expressed in our rim phenotypes (m1 and 5), suggesting that the MS-related CD74⁺ phenotype corresponds to our rim myeloid phenotypes. We confirmed this congruency by staining our MS lesions with anti-CD74, which was expressed predominantly by myeloid cells occupying the lesion rim (Additional file 1: Figure S9). Other attempts to cluster myeloid cells in experimental autoimmune encephalomyelitis (EAE), a mouse model of MS, using single-cell cytometry [25], and in MS lesions using single nuclear RNA sequencing [14], have yielded substantially less myeloid cell heterogeneity.

Our study is limited by the small sample size and the low number of markers, which may result in inaccurate phenotype clustering. Moreover, we acknowledge that no definite conclusions can be drawn from a comparison of two lesions from different individuals. Nevertheless, as a proof-of-concept study it demonstrates the ability of multiplexed tissue imaging and appropriate single-cell analytics to reveal the heterogeneity and spatial properties of glial cell phenotypes in MS lesions. Future applications may combine cell clustering based on single-nucleus RNA sequencing data with highly multiplexed imaging to obtain maximal parametric depth and spatial resolution of phenotypes. This will help define the phenotypes and key interaction networks that drive acute demyelination and chronic low-grade inflammation in established lesions. This may ultimately provide novel targets for therapeutic intervention in relapsing-remitting and progressive MS.

Conclusions

In summary, we found that phenotypic clustering based on differential expression of thirteen glial activation markers produced multiple myeloid cell and astrocyte phenotypes that occupied specific lesion zones. Myeloid cells were activated along a rim-to-center axis, and specific myeloid cell-astrocyte-lymphocyte interactions were present in both lesions. Our study highlights the potential of imaging mass cytometry, paired with novel computational tools, to provide insight into lesion-forming phenotypes and their spatial organization in MS lesions.