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 dysfunction¹. 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 atrophy². Accumulating evidence points to a pivotal role for leptomeningeal inflammation (LMI) in contributing to this pathology^{3, 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 volume⁵. Grey matter pathology (GMP) has been linked to debilitating symptoms such as cognitive impairment and depression⁶, and tends to occur in spatial relation to areas of LMI. Indeed, the most common grey matter lesion location is directly sub-pial⁷, and there is a gradient of increased pathology towards the surface of the brain in MS patients^{8, 9}. Interestingly, GMP occurs without local blood-brain barrier disruption or robust infiltration of peripheral immune cells into the brain parenchyma^{10, 11}. LMI is therefore speculated to be a source of pro-inflammatory molecules that contribute to GMP¹², 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 demyelination¹³. Microarray analysis of cortical lesions from MS cases with LMI revealed a shift in TNF signaling from TNFR1/TNFR2 and NFkB-mediated anti-apoptotic pathways towards TNFR1- and RIPK3-mediated pro-apoptotic/pro-necroptotic pathways¹⁴. 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 sheaths^{15, 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 lesions¹⁷. 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 EAE¹⁸. 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 cells¹⁹. In the parenchyma adjacent to areas of meningeal inflammation we identified astrocytosis, activated microglia, demyelination and evidence of axonal stress and damage¹⁹. 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.

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.

Tissue preparation and spatial gene expression assay

At 12 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 -20C, 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 samples per group, respectively.

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)²⁰. 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 parameters²¹. 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 database^{22, 23} and the clusterProfiler package (v.4.6.0)²⁴. 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)²⁵. 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.

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); spatial feature plots and dimensionality reduction plots were produced with Seurat; tree plots were produced with clusterProfiler, and trajectory heatmaps were produced with SPATA2.

All reported P-values adjusted for multiple comparisons were corrected using the Benjamini-Hochberg method²⁶ 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 nucleus^{19, 27} (Figure 1D). Lesion number remained stable throughout the disease course independent of EAE score (Supplemental Figures 1A-B).

Figure 1:

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 observed¹⁹ (Figure 1E). 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).

We next explored the activity of a broad range of pathways using the pathway responsive genes (PROGENy) method²⁵ (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).

Figure 2

To focus our analyses on foci of meningeal inflammation specifically, we performed unbiased UMAP clustering on the spatial transcriptomic dataset and identified 11 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). We next performed gene set enrichment analysis of C11 spots using the Gene Ontology (GO) database^{22, 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).

Figure 3:

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).

Figure 4:

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). 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 package²⁰. 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).

Figure 5

Discussion

Meninges-restricted inflammation in MS is intimately related to subpial grey matter demyelination³, 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 genesets’ expression pattern 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.

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 antibodies²⁸. 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 supernatants^{29, 30}. B cells are also sources for inflammatory cytokines, such as IL-6 and GM-CSF, and antibodies, which are speculated to contribute to GMP²⁸.

The role of IFNɣ in the pathogenesis of MS and EAE is complex, and likely has stage specific protective and pathologic effects³¹. 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 course³¹. 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 loss³². 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 killing³². Upregulation of genes for self-antigen presentation has also been noted in neurons and oligodendrocyte lineage (OL) cells in post-mortem MS single nucleus RNAseq³³, and OL cells in EAE^{34, 35}.

While SJL EAE models many features of LMI in MS¹⁸, it has several important limitations. 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 parenchyma including thalamus and hypothalamus as opposed to cortical lesions typically seen in MS. Interestingly, 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 factors^{36, 37}. 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 structures³⁸.

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 B cell mediated inflammation, 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 tissues^{39, 40} 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.

Ethics approval and consent to participate

All animal studies followed national and institutional guidelines for human animal treatment in compliance with the Johns Hopkins ACUC.

Availability of data and material

Data availability: spatial transcriptomic data presented in this study will be made publicly available on the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) at the time of publication.

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, SK, SS, MS, PAC, and PB conceived and designed the experiments. SK, SS, 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.

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 (R² < 0.4 in all cases). N = 16, data is representative of three pooled experiments.

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 (F) and UMI (G) 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).

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. 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). EAE specific subclusters and enriched genes were identified in Clusters 1 (B) and 2 (C).