Introduction

Multiple sclerosis (MS) is a chronic autoimmune disease of the central nervous system (CNS) characterized by a relapsing remitting and/or progressive course of demyelination, axonal injury, and neurologic dysfunction1. Highly efficacious disease modifying therapies have revolutionized the prevention and treatment of MS relapses, but are less effective at stopping hallmarks of MS progression such as brain atrophy2. Accumulating evidence points to a pivotal role for leptomeningeal inflammation (LMI) in contributing to this pathology3,4. LMI is found in all subtypes of MS, ranging histologically from disorganized collections of leukocytes to highly organized ectopic lymphoid follicles, and correlates with the presence of cortical grey matter demyelination, neurite loss, and decreased volume5. Grey matter pathology (GMP) has been linked to debilitating symptoms such as cognitive impairment and depression6, and tends to occur in spatial relation to areas of LMI. Indeed, the most common grey matter lesion location is directly sub-pial7, and there is a gradient of increased pathology towards the surface of the brain in MS patients8,9. Interestingly, GMP occurs without local blood-brain barrier disruption or robust infiltration of peripheral immune cells into the brain parenchyma10,11. LMI is therefore speculated to be a source of pro-inflammatory molecules that contribute to GMP12, but the pathway(s) involved remain unknown.

Several putative mechanisms linking LMI to GMP have been proposed. Magliozzi and colleagues identified increased expression of numerous cytokines and chemokines, including interferon gamma (IFNɣ), tumor necrosis factor (TNF), interleukin (IL)-2, IL-22, CXCL13, and CXCL10, in the meninges and CSF of postmortem MS cases with high levels of meningeal inflammation and GM demyelination13. Microarray analysis of cortical lesions from MS cases with LMI revealed a shift in TNF signaling from TNFR1/TNFR2 and I-mediated anti-apoptotic pathways towards TNFR1- and RIPK3-mediated pro-apoptotic/pro-necroptotic pathways14. In marmoset and rat models of experimental autoimmune encephalomyelitis (EAE), sub-pial cortical lesions were found with prominent microglial activation and immunoglobulin deposition on myelin sheaths15,16. Interestingly, complement deposition, one possible mechanism of immunoglobulin-related cellular injury, is not found in purely cortical grey matter lesions, in contrast to white or grey/white matter lesions17. Numerous other mediators of injury, including reactive oxygen and nitrogen species, metabolic stress, and excitotoxicity have been proposed to drive neurodegeneration in MS and may be at play in GMP.

Previous attempts to characterize the relation between LMI and GMP have been limited by the absence of spatially resolved data, which therefore lacks critical information about the anatomic relationship between LMI and the underlying brain parenchyma. Here, we present a spatial transcriptomic dataset and analysis in a mouse model of relapsing/remitting CNS autoimmunity and meningeal inflammation, Swiss Jim Lambert (SJL) mouse EAE18. Our prior work in SJL EAE demonstrated meningeal areas of gadolinium contrast enhancement on magnetic resonance imaging (MRI) that correspond histologically to collections of B cells, T cells, and myeloid cells19. In the parenchyma adjacent to areas of meningeal inflammation we identified astrocytosis, activated microglia, demyelination and evidence of axonal stress and damage19. This work characterizes the spatially resolved transcriptome of LMI in the SJL EAE model system, revealing the subset of genes and gene sets active in the LMI that extend into the adjacent parenchyma and providing insights into immune pathways that could underlie sub-pial neurodegeneration.

Materials and methods

Animals

SJL/J mice were purchased from Jackson Laboratories for all experiments. All mice were maintained in a federally approved animal facility at Johns Hopkins University in accordance with the Institutional Animal Care and Use Committee. Female mice aged 7–8 weeks were used in all experiments and were housed in the animal facility for at least 1 week prior to the start of experiments.

Induction of SJL EAE

Female SJL/J mice were immunized subcutaneously at two sites over the lateral abdomen with 100 μg PLP139-151 peptide with complete Freund’s adjuvant containing 4 μg/ml Mycobacterium tuberculosis H37RA (Difco Laboratories). Mice were weighed and scored serially to document disease course. Scoring was performed using the following scale: 0, normal; 1, limp tail; 2, hind limb weakness; 3, hind limb paralysis; 4, hind limb and forelimb weakness; and 5, death. CFA only controls received CFA without PLP139-151 peptide.

MRI imaging

At weeks 6, 8, and 10 post immunization, a horizontal 11.7 T scanner (Bruker BioSpin) with a triple-axis gradient system (maximum gradient strength = 740 mT/m), 72 mm volume transmit coil and 4-channel receive-only phased array coil was used to image the mouse forebrain. During imaging, mice were anaesthetized with isoflurane together with mixed air and oxygen (3:1 ratio) and respiration was monitored via a pressure sensor and maintained at 60 breaths/min. Before imaging, 0.1 ml diluted Magnevist® (gadopentetate dimeglumine, Bayer HealthCare LLC, 1:10 with PBS) was injected. Scans were then analyzed to identify areas of meningeal contrast enhancement by two independent examiners (P.B., S.K.). We counted the number of areas of meningeal contrast enhancement on each individual MRI slice and used the cumulative number to represent the amount of meningeal contrast enhancement. All quantifications were performed by at least two independent examiners and their scores were averaged.

Immunofluorescence preparation and analysis

At 11 weeks post-injection, CFA only and age-matched naïve animals were perfused with 1x PBS followed by 4% PFA. Brains were dissected, post-fixed for 48 hours, and then cryopreserved in 30% sucrose for 48 hours. Tissue was frozen and then sectioned (12µm thick) directly onto slides using a cryostat. For immunofluorescent staining, slices were first blocked in blocking buffer (5% normal goat serum in 0.4% Triton™ X-100 in PBS) for 1 hour at room temperature, followed by overnight incubation in primary antibody. The next day, samples were washed three times in PBS and then incubated with secondary antibody for 1 hour at room temperature. Primary antibodies used in this study were Iba-1 at 1:300 dilution (cat no. 019-19741, Wako Chemicals), CD68 at 1:300 dilution (cat no. 14-0681-82, Invitrogen), and GFAP at 1:400 dilution (cat no. GA52461-2, Agilent Dako). Secondary antibodies were Alexa™ fluorophores (cat no. A-21244 and A-11006, Invitrogen), all at 1:1000 dilution.

Images were acquired on a Zeiss Axio Observer Z1 epifluorescence microscope and analyzed using FIJI (NIH). Images were centered at the interpeduncular and quadrigeminal cisterns, where LMI typically develop in SJL-EAE, and included the lateral thalamus. Images were acquired and analyzed by an experimenter blinded to the group of each sample and data presented in the form of mean fluorescence intensity (MFI).

Tissue preparation and spatial gene expression assay

At 11 weeks post-immunization, animals were euthanized in CO2 chamber and perfused with cold PBS. Brains were dissected and one hemisphere drop fixed in isopentane cooled on dry ice. Fresh frozen brain samples were then cut coronally at a thickness of 10μm and placed on the capture area of Visium Gene Expression slides (v1; 10x Genomics). Each slide contained four 6.5 mm x 6.5 mm capture areas. Sample preparation was carried out according to manufacturer’s instructions. After fixation with methanol at −20°C, Hematoxylin and Eosin (H&E) staining was performed for morphological analysis and spatial alignment of the sequencing data. After the enzymatic permeabilization, mRNA was captured by probes and cDNA generated. Barcoded cDNA was isolated using SPRIselect-cleanup (Beckman Coulter) and amplified. Amplified cDNA was fragmented and subjected to end-repair, poly-A-tailing, adaptor ligation, and 10x-specific sample indexing as per the manufacturer’s instructions. Following assessment of RNA quality, sequencing was performed using a Novaseq S2 100. Brains from four naïve and four EAE mice were used to prepare 5 and 6 individual slices per group, respectively, with once naïve mouse contributing two slices and two EAE mice contributing two slices.

Spatial transcriptomics data processing

Each sample went through identical quality control processing steps. SpaceRanger software (v.1.3.1) was used to pre-process the sequencing data, aligning to the mm10-2020A reference transcriptome. Feature barcoded expression matrixes were used as input for downstream spatial transcriptomics analysis using Seurat (v.4.3.0) and SPATA2 (v.0.1.0)20. Data was loaded and analyzed using the Seurat, and all spots that were determined to not be over tissue were discarded using the filter.matrix option in Seurat’s Load10XSpatial function. Spots with less than 250 measured genes and less than 500 unique molecular identifiers (UMIs) were filtered out. Data normalization and stabilization of sequence depth variance was performed on each sample using SCTransform with default parameters21. Sample data was then annotated and combined into a merged object for downstream quality control and analysis. Dimensionality reduction was performed using principal component analysis (PCA), followed by computation of shared nearest neighbors of the first 10 principal components and cluster identification (resolution 0.3). To visualize all spots in a two-dimensional plot, a UMAP was created with Seurat’s RunUMAP function using the first 10 principal components. Cluster and sub-cluster enriched genes were identified with Wilcoxon tests as implemented in the FindMarkers function within Seurat. Differential gene expression between groups was attained using DESeq2 (v.1.38.2) on samples pseudobulked by biological replicate. Gene set enrichment analysis and visualization was performed using the gene ontology database22,23 and the clusterProfiler package (v.4.6.0)24. Estimated signaling pathway activities were calculated for each spot with the top 500 genes of each pathway on SCTransformed data using the PROGENy package (v.1.20.0)25. For subclustering analysis, select clusters were subsetted before undergoing dimensionality reduction, neighbor calculation, cluster identification, marker identification, and gene set enrichment as described above. Trajectory gene and gene set modeling analysis was performed with the SPATA2 package. Spatial trajectories were drawn in EAE slices from the center of cluster 11 to the centro-medial nucleus of thalamus. Gene and gene-sets were analyzed along these trajectories using the assessTrajectoryTrends function within SPATA2.

RNAscope

Brain tissue was collected from N = 3 naïve or EAE (11 weeks post immunization) SJL mice after transcardial perfusion with 4% PFA. Brains were post-fixed for 48 hours and then dehydrated via sucrose gradient over 48 hours before freezing in OCT. Fourteen µm thick tissue sections were collected directly onto Superfrost+ slides (Fisherbrand, Cat. No. 22-037-246). The RNAscope assay was then carried out using the Multiplex Fluorescent Reagent Kit v2 (Advanced Cell Diagnostics, Cat. No. 323100) per the manufacturer’s instructions. Briefly, tissue sections underwent dehydration in an ethanol gradient, peroxide blocking, and antigen retrieval steps. Antigen retrieval was performed using a steamer at 100°C and the epitope retrieval solution provided by the manufacturer. The mRNA in the tissue was hybridized to RNAscope probes (Advanced Cell Diagnostics) against: C3 (Cat. No. 417841-C3), Fcgr3 (Cat. No. 587241), B2m (Cat. No. 415191-C2), Cd74 (Cat. No. 437501-C3), and Gfap (Cat. No. 313211) for 2 h at 40°C and stored overnight in 5x saline sodium citrate. Next, amplification steps were performed according to the manufacturer’s instructions. Fluorescent labeling was performed with TSA Vivid Fluorophore 520, TSA Vivid Fluorophore 570, and TSA Vivid Fluorophore 650 (all at 1:2000 dilution). Slides were washed and mounted with ProLong Gold Antifade Mountant (ThermoFisher, Cat. No. P36830). Slides were imaged via confocal microscopy (Ziess LSM900).

Positive cells were manually labeled in ImageJ (NIH). The x-y coordinate of each cell was used to determine the shortest distance to areas of leptomeningeal inflammation in R. Distances were batched at intervals of 50μM. Non-linear regression of the batched data was computed using the stats::nls function and the formula percentage ∼ a * distance ^ b.

RNAseq statistics and data visualization

Plots were generated and statistics calculated using R (v.4.2.2) and Rstudio (v.2022.07.2). Bar plots, box plots, MA plots, Venn diagram, and dot plots were produced with ggplot2 (v.3.4.0); heat maps were produced with pHeatmap (v.1.0.12) and enrichplot (v.1.18.4)26; spatial feature plots and dimensionality reduction plots were produced with Seurat; tree plots were produced with enrichplot, and trajectory heatmaps were produced with SPATA2.

All reported P-values adjusted for multiple comparisons were corrected using the Benjamini-Hochberg method27 unless otherwise specified. The number of samples in the EAE and naïve groups was chosen based on expected levels of variability in prior experiments.

Results

We tracked the development of LMI during SJL EAE using contrast enhanced serial MRI imaging. MRI data was then used to target areas of LMI for spatial transcriptomics (Figure 1A). Mice developed a characteristic relapsing pattern of neurologic impairment, and contrast-enhanced MRI was performed at weeks 6, 8, and 10 post immunization (Figure 1B-C). In the SJL EAE model, contrast enhancing meningeal lesions are most frequently found in the intrapeduncular cistern, quadrigeminal cistern, and the cleft between the hippocampus and medial geniculate nucleus19,28 (Figure 1D). Lesion number remained stable throughout the disease course independent of EAE score (Supplemental Figure 1A–B).

MRI guided spatial transcriptomics of meningeal-based inflammation in SJL EAE.

(A) Schematic describing the experimental paradigm. SJL mice underwent brain MRI 6-, 8-, and 10-weeks post immunization with PLP139-151. Brain slices from regions with meningeal inflammation were collected and processed for spatial transcriptomics on the 10x Genomics platform. (B–C) Behavior scores (B) and mouse weights (C) of the EAE cohort. Red arrows indicate MRI time points, black arrow indicates time of tissue harvesting (N = 4). (D) Representative post-contrast MRI brain images, white arrows indicate areas of meningeal-based inflammation. (E) Representative images of H&E-stained tissue sections mounted on spatial transcriptomics slides (left, naïve; right, EAE). (F) Spatial feature plots from naïve (top row) and EAE (bottom row) representative samples demonstrate altered expression of genes related to complement (i), immune infiltration (ii), antigen presentation (iii), dendritic spines, and astrocyte activation (iv).

Brain slices were collected from four naïve mice and four EAE mice 11 weeks post-immunization and used to prepare five and six samples for spatial transcriptomics, respectively. H&E staining confirmed the presence of meningeal inflammation in the areas of contrast enhancement, as we had previously observed19 (Figure 1E). Since naïve animals were used as controls, we confirmed that CFA alone does not produce lasting glial reactivity or LMI formation. Groups of animals were given CFA only or left naïve. Neither group developed neurologic signs, and after 11 weeks the brains were processed for IHC analysis. There was no evidence of LMI development, and no difference in glial reactivity as measured by GFAP, IBA1, or CD68 intensity (Supplemental Figure 1C–E).

Spatial Transcriptomic data from all samples were of high quality and read depth when assessed by treatment group (Supplemental Figure 2A–C) or sample (Supplemental Figure 2D– E). There was expected anatomic variability in number of read counts and number of features per spot, with relatively low counts and features in white matter as compared to the cortex or hippocampus (Supplemental Figure 2F–G), and spots had similar degrees of average complexity between EAE and naïve samples (Supplemental Figure 2H). UMAP clustering revealed no significant independent effect of sample (Supplemental Figure 2I) or slide (Supplemental Figure 2J). Numerous genes were differentially expressed between EAE and naïve slices as estimated by DESeq2 on samples pseudobulked by biological replicate (Supplementary Table 1; Supplemental Figure 2K), including genes associated with the complement cascade, immune infiltration, antigen presentation, and astrocyte activation (Figure 1F; Supplementary Figure 2L).

We next explored the activity of a broad range of pathways using the pathway responsive genes (PROGENy) method25 (Figure 2A). Inflammatory pathways related to TNF, JAK-STAT, and NFκB signaling were upregulated in EAE compared to naïve, with peak activity around sites of meningeal inflammation (Figure 2B–C). Pathways related to Trail and PI3K were downregulated in EAE compared to naïve (Figure 2D–E), and TGFβ pathway activity was unchanged between groups (Figure 2D–E).

PROGENy analysis reveals spatially restricted pathway activity differences between naïve and EAE.

(A–B) Heat map displaying averaged PROGENy pathway analysis results. (C) Representative spatial plot showing activity of the JAK-STAT, NFkB, and TNFa signaling pathways. (D) Comparison of JAK-STAT, NFkB, and TNFa pathway activities between groups. (E) Representative spatial plot showing activity of the Trail, PI3K, and TGFb signaling pathways. (F) Comparison of Trail, PI3K, and TGFb pathway activities between groups. (Naïve mouse N = 4, sample N = 5; EAE mouse N = 4, sample N = 6; multiple T tests corrected for multiple comparisons with the Benjamini, Krieger, and Yekutieli method; * p < 0.05, ** p < 0.01, *** p < 0.001).

To focus our analyses on foci of meningeal inflammation specifically, we performed unbiased UMAP clustering on the spatial transcriptomic dataset and identified 12 distinct clusters (Figure 3A). Grouping the dimensional reduction UMAP plot by EAE and naïve revealed that cluster eleven (C11) was restricted to EAE samples (Figure 3B). C11 makes up 1–5% of the total spots in those samples (Figure 3C), or about 20–120 total spots, and was significantly enriched in EAE relative to naïve samples (Figure 3D). Visualizing the clusters using spatial feature plots confirmed that most clusters map consistently to specific anatomic regions and were similar between naïve and EAE (Figure 3E; Supplementary Figure 3A). C11 overlapped with previously noted areas of meningeal MRI enhancement, suggesting that C11 represents areas of meningeal inflammation (Figure 3E; Supplementary Figure 3B–D). We compared the gene expression in C11 to other clusters and found 132 upregulated genes and 70 downregulated genes (p < 0.05; Log2 fold change > 1; Figure 3F). Inflammation-related genes were prominently represented in C11 relative to other clusters (Supplementary Table 2), with top genes including Cd74, C3, and Gfap (Figure 3G; Supplementary Figure 3E). The presence of glial genes such as Gfap within cluster 11 likely represents areas of subpial brain parenchyma included within the cluster. Immunoglobulin genes, which are highly expressed in areas of human grey matter pathology,14 were also upregulated in C11 (Supplementary Figure 3F). We next performed gene set enrichment analysis of C11 spots using the Gene Ontology (GO) database22,23, finding 538 variable gene sets (adjusted P-value < 0.05; Supplementary Table 3). Among the most prominently enriched gene sets were those involved in antigen processing and presentation, complement activation, lymphocyte activation, and cytokine production and response (Figure 3H).

Unbiased clustering reveals a group of spots enriched in inflammatory genes.

(A–B) UMAP dimensionality reduction plots colored by (A) cluster or (B) group. (C) Bar plot showing the proportion of spots in each cluster by sample. (D) Number of spots in cluster 11 by group (N = 11; Student’s two-tailed T test). (E) Representative spatial feature plots of naïve and EAE samples showing the spatial distribution of each cluster. (F) MA plot comparing differences in gene expression between cluster 11 and all other clusters averaged across samples. Red and blue spots represent genes in cluster 11 that are significantly increased or decreased, respectively (adjusted p-value < 0.05, log 2 fold change > 1). (G) Bar plot of top 15 genes enriched in cluster 11 compared to other clusters. (H) Tree plot displaying gene set enrichment results using the gene ontology (GO) database. Spots in cluster 11 were compared to other spots and gene set sizes ranging from 10-500 were included (adjusted p-value < 0.05).

After establishing the defining transcriptomic features of meningeal inflammation in our model, we next sought to characterize inflammatory changes in the adjacent CNS parenchyma. We performed unbiased sub-clustering of each cluster in turn and found EAE-specific subclusters and differentially expressed genes within clusters 1 and 2. These subclusters were labelled 1_3, 1_4, and 2_6 (Figure 4A–B; Supplemental Figure 4) and represent regions of the thalamus and hypothalamus (Figure 4C).

Subclustering of spots adjacent to meningeal immune follicles reveals a subset of active immune patterns.

(A) UMAP dimensionality reduction plots showing subclustering of cluster 1 colored by (left) group or (right) cluster. (B) UMAP dimensionality reduction plots showing subclustering of cluster 2 colored by (left) group or (right) cluster. (C) Representative spatial feature plot showing the locations of cluster 1 and 2 subclusters. (D) Distance from the center of indicated subclusters to the nearest point of cluster 11 (N = 11; Student’s two-tailed T-Test). (E–G) Tree plot displaying gene set enrichment results using the gene ontology (GO) database for subcluster 1_3 (E), 1_4 (F), and 2_6 (G) compared to other spots in their respective clusters. (H) Venn diagram shows overlap of significantly enriched GO gene sets between cluster 11 and subclusters 1_3, 1_4, and 2_6, with (I) 31 gene sets elevated in all. GO gene set of size ranging from 10-500 were included (adjusted p-value < 0.05).

Next, we assessed whether EAE-specific subclusters were physically closer to meningeal inflammation than other related subclusters. The distance between the average location in each subcluster to the nearest point of C11 was calculated. In the subclusters of cluster 1, we found no difference in proximity to C11 between EAE-specific subclusters vs. subclusters present in EAE and naïve. However, in subclusters of cluster 2, the EAE-specific subcluster was significantly closer to C11 on average compared to other subclusters (Figure 4D). To explore which pathways were activated in these subclusters, we performed gene set enrichment analysis using the GO database (Figure 4E–G; Supplementary Tables 4–6). All subclusters displayed enrichment of GO genesets related to inflammation, antigen processing and presentation, and humoral immunity. Interestingly, pathways related to neuron death, cellular stress, and negative regulation of cellular processes were also enriched (Figure 4E–G; Supplementary Figure 5). Interestingly, other pathways related to cell death, including positive regulation of apoptotic process or programmed cell death, were upregulated in subcluster 1_3 but downregulated in subcluster 1_4, suggesting spatially variable survival signals (Supplementary Tables 4-5). Overall, 31 upregulated pathways were conserved between cluster 11 and subclusters 1_3, 1_4, 2_6 (Figure 4H). We found prominent conserved upregulation of pathways related to adaptive/B-cell mediated immunity, antigen processing and presentation, cell killing, and others (Figure 4I).

Our pathway analysis of meningeal inflammation and areas of inflamed adjacent brain parenchyma suggested that inflammatory signals increased in meningeal inflammation could have variable ‘penetration’ into the adjacent brain. We sought to test this using spatial trajectory gene/gene set expression modelling available within the SPATA2 software package20. Trajectories were drawn in EAE samples from the largest region of meningeal inflammation to the central thalamus (Figure 5A). Gene and gene set expression levels were then evaluated along the length of these trajectories and compared to ideal patterns of expression, as demonstrated for representative genes B2m and C3 (Figure 5B–E). The difference between gene or gene set expression and the ideal patterns, here “logarithmic descending” or “gradient descending”, is represented by the residual line (Figures 5B, 5E). The area under the residual curve (residual AUC) is therefore inversely proportional to fit for the given gene or gene set and ideal pattern. C3 expression declines rapidly along the trajectory and the logarithmic descending residual AUC is lower than gradient descending residual AUC (Figure 5E), while B2m follows a less steep decline and fits the two patterns similarly (Figure 5C). Trajectory analysis of gene sets enriched in meningeal inflammation and adjacent brain parenchyma was performed in this way, and the average residual AUC was calculated for gradient descending/ascending and logarithmic descending/ascending patterns (Figure 5Fi). As expected, all gene sets fit descending patterns better than ascending ones. There was variability in fit to the gradient descending pattern of expression, representing a more gradual decline in pathway enrichment along the trajectory, with gene sets related to antigen processing and presentation, cell killing, interleukin 6 production, and interferon gamma response having the best fit (Figure 5Fi). Enrichment score trajectory heatmap of a representative sample corroborates this, showing increased activity farther along the spatial trajectory for these gene sets (Figure 5Fii).

Trajectory analysis reveals gradients of gene expression originating from meningeal lymphoid follicles.

(A) Trajectories were drawn based on spatial cluster plot (i) from C11 to C2 (ii). (B) Representative plot of B2m relative expression along the trajectory length. Green line: B2m expression; black line: ideal model fit, “logarithmic descending” (top) or “gradient descending” (bottom); red line: residual area under the curve (AUC) representing the difference between B2m expression and the ideal model. (C) Barplot showing residual AUC of B2m relative expression along the trajectory direction compared to “logarithmic descending” or “gradient descending” (Student’s two-tailed T-test). (D) Representative plot of C3 relative expression along the trajectory length. Green line: C3 expression; black line: ideal model fit, “logarithmic descending” (top) or “gradient descending” (bottom); red line: residual area under the curve (auc) representing the difference between C3 expression and the ideal model. (C) Barplot showing residual AUC of C3 relative expression along the trajectory direction compared to “logarithmic descending” or “gradient descending” (Student’s T-test). (F) Genesets that were previously identified as significantly enriched in C11 were selected for trajectory analysis. Residual AUCs were calculated for “logarithmic descending”, “gradient descending”, “logarithmic ascending”, and “gradient ascending” ideal fits and displayed on (i) a heatmap sorted by “gradient descending”. (ii) Representative feature plot demonstrating deeper penetration of upper genesets (related to antigen presentation and processing, microglial activation, IL-6 production, interferon gamma) response relative to other gene sets (B cell activation, T cell activation, TNF production, complement, humoral immune response). (G, I) Representative images of RNAscope labeling for (G) Fcgr3, B2m, Cd74, and (I) Gfap, C3 in SJL mice 11 weeks after EAE induction. Yellow dashed lines indicate the areas of leptomeningeal inflammation, scale bars represent 100μM. (H, J) Barplots representing the percent of marker-positive cells present at distances from leptomeningeal inflammation. Lines represent best fit curves from exponential regression. N = 3 animals per group; bars represent mean, error bars represent standard error.

We next validated this variability in pattern of expression in a separate cohort of SJL mice with EAE using RNAscope to label selected transcripts related to glial activation (Fcgr3 and Gfap) antigen presentation (B2m and Cd74) and complement (C3) (Figure 5G, I). Consistent with our spatial transcriptomics data, each of these transcripts was substantially induced in EAE as compared to Naïve (Supplemental Figure 6). We quantified the shortest distance of each target-positive cell from the region of LMI (Figure 5H, J). For all transcripts, the number of positive cells decreased with distance exponentially. Exponential regression was applied to find the best fit curves of the data, modeled by the equation y = a * xb where y is percentage of cells, x is distance from LMI. The exponential constant b reflects the rate of decrement and had higher absolute value in C3 (b = −1.13) relative to B2m (b = −0.59) (Table 7).

Discussion

Meninges-restricted inflammation in MS is intimately related to subpial grey matter demyelination3, atrophy, and neurocognitive symptoms, but therapeutically targeting this aspect of the disease is challenging due to poorly understood pathologic mechanisms. Here, we present spatial transcriptomics analysis in a mouse model of LMI, finding a broad swath of inflammatory pathways upregulated at foci of LMI and a subset of them upregulated in the nearby brain parenchyma. Notably, genes related to B cell mediated responses and antigen processing and presentation were upregulated in inflammatory parenchymal subclusters. Variable penetrance of inflammatory pathway activity was evident when analyzing the expression pattern of genesets along a linear trajectory from LMI into the CNS parenchyma, where antigen processing and presentation, interleukin 6 production, and the IFNɣ response pathway activity followed a more gradual decline as compared to other pathways.

An important limitation of the spatial transcriptomics method used in this study is the spatial resolution of each spot. Since each spot is ∼55um in diameter it is likely that multiple cells are captured in each one. The simultaneous addition of single cell analysis would facilitate deconvolution of cell types within each spot. Furthermore, some spots of the borders of anatomic structures are likely to include both regions. This is exemplified by the spots within cluster 11, which in this study mainly represent meningeal inflammation, also being enriched in glial genes such as Gfap. Other methods of spatial RNA analysis exist that allow single cell resolution, but to date are probe-based and therefore limited by the number of genes assessed.

Prior work supports roles for pathways identified in our dataset, including B cell and IFNɣ mediated responses, in contributing to neurodegeneration in MS. LMI in MS is rich in B cells, and the critical role of B cells in MS has been underscored by the success of B cell depletion in relapsing and progressive MS. Recent studies have proposed numerous mechanisms whereby B cells could contribute to cortical pathology, including indirectly through activation and inflammatory polarization of T cells, myeloid cells, and astrocytes or directly through production of neurotoxic cytokines or antibodies29. B cell culture supernatants from MS patients, but not healthy controls, are toxic to rat and human neurons and oligodendrocytes with this effect being mediated by the extracellular vesicle (EV) fraction of the supernatants30,31. B cells are also sources for inflammatory cytokines, such as IL-6 and GM-CSF, and antibodies, which are speculated to contribute to GMP.29

The role of IFNɣ in the pathogenesis of MS and EAE is complex, and likely has stage specific protective and pathologic effects32. Its prominent upregulation in EAE was initially considered evidence of its pathogenic nature, but subsequent experiments showed that it was pathogenic during the initiation phase but protective later in the course32. IFNɣ signaling and subsequent upregulation of antigen processing and presentation on glia has been identified as a mechanism that could lead to remyelination failure and subsequent neuronal loss33. Oligodendrocyte precursor cells (OPCs) upregulate antigen presentation and cross-presentation pathways in response to IFNɣ, promoting inflammation and making them susceptible to CD8+ T cell killing33. Upregulation of genes for self-antigen presentation has also been noted in neurons and oligodendrocyte lineage (OL) cells in post-mortem MS single nucleus RNAseq34, and OL cells in EAE3538. Recent work also assessed subcortical white matter lesions from postmortem MS cases via spatial transcriptomics, finding prominent elevation of TNF signaling, in agreement with presented results in the SJL EAE model.39

Another important consideration in these experiments is our choice of naïve, rather than CFA only, controls. While often used as the control in EAE studies focused on mechanisms of autoimmunity, CFA only can independently induce systemic inflammation. Since this study seeks to describe transcriptomic changes in neuroinflammation more broadly, we chose to use a healthy comparison group to maximize our ability to find genes enriched in neuroinflammation. Ultimately, however, the choice of naïve or CFA only controls is unlikely to have affected our conclusions. SJL-EAE, unlike the more common C57Bl6-EAE, does not require pertussis toxin during the induction. The only difference between naïve and CFA only controls is the subcutaneous CFA delivered at time of immunization (11 weeks prior to experiment endpoint). Indeed, when we compared CFA only and healthy animals at 11 weeks there was no difference in glial reactivity by GFAP, IBA1, or CD68 MFI. There was also no evidence of neurologic symptoms or LMI development in CFA only controls.

While SJL EAE models many features of LMI in MS18, there are important differences that limit the direct translation of these results to MS. The majority of LMI identified in mice with SJL EAE in our experiments occurred in subarachnoid cisterns, and as a result prominently affected areas of deep grey matter including thalamus and hypothalamus as opposed to cortical lesions also seen in MS. Notably, neuronal loss in the thalamus and other deep grey nuclei does occur in MS with a similar ‘surface-in’ gradient of neuronal injury, thought to also be related to toxic CSF-derived factors40,41. Demyelination, microglia and astrocyte activation, and neuronal loss are evident in the parenchyma adjacent to LMI in SJL EAE19,42. We also noted some variability in the extent and location of LMI, which is unavoidable in EAE. While other animal models of LMI, such as directly injecting inflammatory cytokines into the meninges/cortex, produce more predictable regions of LMI, they involve traumatic injury to the brain and typically lack follicle-like structures43. Notably, novel rodent models of cytokine-mediated LMI have recently been developed, including virus-mediated meningeal overexpression of TNF and IFNɣ or sub-arachnoid injection of recombinant Lymphotoxin-α, were recently developed and exhibit lymphoid like structures.44,45 Spatially resolved analyses of these models could provide additional insights into sub-pial pathology during neuroinflammation.

This work is the first to characterize a mouse model of LMI and grey matter injury using spatial transcriptomics, and in addition to the analysis presented here contributes a publicly available dataset for future research. We highlight the importance of antigen processing and presentation and complement signaling, which are prominently upregulated in sub-pial grey matter, in our model. Future studies should focus on spatial transcriptomics in post-mortem or biopsied human tissue. While access to appropriate samples remains a significant barrier, recent advances in RNA extraction from formalin fixed paraffin embedded tissues46,47 will allow for the use of large banks of historically collected and preserved samples and could dramatically improve availability.

Acknowledgements

We thank members of the Calabresi lab for their valuable comments during multiple discussions of this work. We thank the Johns Hopkins Medicine Single Cell and Transcriptomics Core for their assistance with planning and data acquisition. We acknowledge the contribution of animals used in this research study.

Availability of data and material

Data availability: spatial transcriptomic data presented in this study is available on the Gene Expression Omnibus data (series record GSE236963).

Competing interests

The authors have no competing interests to disclose.

Funding

This work was supported in part by an investigator-initiated grant from EMD-Serono to PB, a Harry Weaver Neuroscience Scholar award to PB, and a fellowship grant from the National Multiple Sclerosis Society and the American Brain Foundation (FAN-2106-37832) to SPG.

Authors’ contributions

SPG, SS, SK, MS, PAC, and PB conceived and designed the experiments. SS, SK, JH, and MS performed the experiments and collected the data. SPG and MS analyzed the data. SPG, SS, PAC, and PB wrote the manuscript. All authors reviewed and suggested improvements to the manuscript.

Figure Legends

Supplemental Figure 1

Contrast enhancing meningeal inflammation in SJL EAE does not change at chronic time points or with clinical disease scores. (A) Box plot of lesion number versus week post immunization. ANOVA with Tukey’s HSD post-hoc test. (B) lesion number plotted against EAE score demonstrates no strong correlation (R2 < 0.4 in all cases). N = 16, data is representative of three pooled experiments. (C–E) Immunofluorescence analysis of glial reactivity in naïve and CFA only groups, 11 weeks after injection. Representative staining at low and high magnification and MFI (Mean fluorescence intensity) calculations of (C) GFAP, (D) IBA1, (E) CD68. N = 5, Mann Whitney Tests.

Supplemental Figure 2

Quality control analysis of spatial transcriptomic data. (A) Density of read counts, (B) number of unique molecular identifiers (UMI), (C) and the ratio of read count/UMI per spot for naïve and EAE samples. Most spots had > 500 read count (A), > 250 UMI (B), and > 0.80 log(nCount/UMI) (C). (D) Dot plot shows most spots having linear correlation between UMI and read count. (E–F) Violin plots of read counts (E) and UMI (F) per sample. (G-H) Representative spatial feature plots of read count (G) and UMI (H) demonstrate expected anatomic variability in in transcript amount and diversity. (I–J) UMAP dimensionality reduction plots, colored by mouse (I) or slide (J). (K) MA plot of genes enriched in EAE compared to naïve samples. Data was pseudobulked by sample and groups were compared using DESeq2 for gene enrichment (adjusted p-value < 0.05, log 2 fold change > 1). (L) Comparison of average expression for select transcripts in naïve and EAE samples. N = 4, Student’s T test.

Supplementary Figure 3

Spatial and transcriptional properties of inflammatory clusters. (A) Representative spatial feature plots of naïve (top row) and EAE (bottom row) samples colored by cluster show consistent labeling of anatomic regions and the distribution of meningeal inflammation (cluster 11). (B–D) Representative images of gadolinium enhanced MRI scans collected 10 weeks after immunizations and samples collected from the same animal. Areas of meningeal enhancement (red arrows) correlate with cluster 11. (E) Dot plot of the top 100 significantly enriched genes (sorted by adjusted P value) in cluster 11 compared to other clusters. Dot color represents expression value and dot size represents the percent of spots within that cluster expressing the gene.

Supplementary Figure 4

Subcluster analysis reveals EAE-specific subclusters and gene enrichment. (A-K) Data from each individual cluster was extracted and unbiased subclustering performed based on gene expression; (A) cluster 0, (B) cluster 1, (C) cluster 2, (D) cluster 3, (E) cluster 4, (F) cluster 5, (G) cluster 6, (H) cluster 7, (I) cluster 8, (J) cluster 9, (K) cluster 10. UMAP plots are shown colored by subcluster (left panel) and group (middle panel). Data from each cluster was isolated, pseudobulked by sample, and groups were compared using DESeq2 for gene enrichment (right panel; adjusted p-value < 0.05, log 2 fold change > 1).

Supplementary Figure 5

Gene contributions to enriched genesets in subclusters 1_3, 1_4, and 2_6. (A-C) Heatmaps of GO genesets enriched in subcluster 1_3 (A), 1_4 (B), and 2_6 (C). Individual genes contributing to genesets are displayed on the X axis. Data from each subcluster was pseudobulked by sample and groups were compared using DESeq2 for gene enrichment (adjusted p-value < 0.05, log 2 fold change > 1). Heatmaps of select genesets were generated using the enrichPlot package.

Supplementary Figure 6

Representative images of RNAscope labeling in naïve and EAE tissues. RNAscope labeling of (A) Fcgr3, B2m, Cd74, and (B) Gfap, C3. Images are representative of N = 3 animals per group. Yellow dashed lines outline the subarachnoid space. Scale bars represent 100μM (main image) and 10μM (inset).