Abstract
The cuprizone (CPZ) model is widely used for modeling demyelination in multiple sclerosis (MS) and for testing potential remyelination therapies. To better understand the underlying pathology of the CPZ model and evaluate its translatability, we integrated single-cell and spatial transcriptomics (ST) to investigate spatial cellular and molecular interactions during de- and remyelination in multiple brain regions. ST revealed global demyelination and neuroinflammation in the brain beyond the corpus callosum (CC), with region-specific differences. We identified oligodendroglia and microglia as two major cell types with significant transcriptomic changes in the model. CPZ-associated subclusters of oligodendroglia (marker genes Arap2, Dock10, Tenm4, Pex5l and Dock1) and microglia (marker genes ApoE, Axl, Cd9 and Lpl) were mapped to the CC by ST. During remyelination, while mature oligodendrocytes (MOL) nearly reversed their phenotype back to the control state, microglia remained associated with the demyelination phenotype. Ligand‒receptor (LR) pairing analyses predicted growth factor and phagocytic pathway enrichment during demyelination, which is consistent with changes in MS lesions, and microglia were predicted to be the major sender cells. LR pairing also predicted a high likelihood of interaction between oligodendroglia and microglia, and a novel interaction between MOL and oligodendrocyte precursor cells (OPC), underscoring their roles during de- and remyelination. Finally, astrocytes in the CPZ model had the greatest preservation of disease-associated modules in MS lesions, while MOL, OPC, and microglia showed moderate to low preservation, which overall suggests that the CPZ model has moderate translatability to chronically active MS lesions.
Introduction
Multiple sclerosis (MS) is an inflammatory demyelinating disease of the central nervous system (CNS) with pathological hallmarks of gray and white matter demyelination, neurodegeneration, inflammation, and glial reactions (Lassmann, 2018; Reich et al., 2018). A genome-wide association study (GWAS) previously highlighted the role of peripheral immune cells in MS disease risk (). Indeed, approved MS therapies targeting peripheral inflammation have shown high efficacy in dampening disease activities. However, despite well-controlled acute inflammation, disease progression can continue unabated, resulting in increased disability that currently lacks effective treatments (). A recent GWAS revealed a locus associated with MS progression (). Notably, the locus does not overlap with previously identified immunological susceptibility genes. The nearest gene to the locus, dysf, encodes DYSFERLIN, a protein specifically expressed by oligodendrocytes and neurons, underscoring CNS tissue involvement in MS disease progression. Therefore, it is critical to understand the mechanisms underlying the progressive biology at the cellular and molecular levels in the CNS and pave the paths for identifying potential drug targets to treat MS disease progression.
The cuprizone (CPZ) model is widely used to model the MS demyelination process and/or test potential new therapies for remyelination (Vega-Riquer et al., 2019). CPZ is a copper chelator that inhibits mitochondrial functions, causing the selective degeneration of mature oligodendrocytes (MOLs) (). The model has a stereotypical spatial-temporal demyelination-remyelination course, and the corpus callosum (CC), which is the most studied brain region in the CPZ model, shows a highly predictable demyelination pattern. During CPZ intoxication, demyelination and remyelination events may coexist at the same time as the lesion evolves until peak demyelination is reached. This feature has been considered a strength of the model to mimic early-stage MS lesions when damage and repair can occur concurrently (; Ransohoff, 2012). The CNS spatial phenotypes in the CPZ model, such as demyelination, neurodegeneration, and gliosis, are traditionally characterized by histological staining, immunohistochemistry (IHC), and electron microscopy (Stidworthy et al., 2003; ; ). While histopathological examination preserves findings relating to anatomical architecture, the disadvantage is that only a few selected markers can be examined, and the overall cellular and molecular response is not fully captured.
In recent years, both the CPZ model and MS have been investigated by omics-based efforts to better understand the underlying pathology of neuroinflammation and demyelination (; ; Lerma-Martin et al., 2024). However, these studies usually focused on one type of omics and investigated either the CPZ model or MS with the technology. In our current study, we integrated multiple types of omics data with several methodologies to analyze the CPZ model and evaluate the CPZ model and MS datasets together when appropriate. We aimed to unbiasedly unwrap cellular and molecular changes during de- and remyelination in the CPZ model by employing spatial transcriptomic (ST) technology (Liu et al., 2021; Rao et al., 2021) to preserve the anatomical context. We have resolved spatial gene expression in the model by integrating the ST data with the single-nucleus RNA-seq (snRNA-seq) reference in the same study. We found that demyelination and neuroinflammation were prominent and widespread in the CPZ brain with brain region-specific gene expression changes. Ligand‒receptor (LR) pairing analyses revealed an enrichment of phagocytic and growth factor pathways in the CPZ group from both the snRNA-seq and ST data. Finally, we compared the snRNA-seq data from our study with a human MS dataset () and found a moderate level of preservation of MS-associated disease modules in the CPZ model. In the future, the methodology developed in this study can be applied to human MS samples to map expression patterns and changes in the heterogeneous lesion environment ().
Results
We collected brain tissues (approximately bregma −1.5 mm) from three experimental groups of mice (∼10 weeks old): (1) the cuprizone group (CPZ), in which mice were fed with 0.3% CPZ for 4 weeks to induce peak demyelination, (2) the recovery group (RCV), in which mice received 4 weeks of CPZ followed by 2 weeks of normal chow to allow for spontaneous remyelination, and (3) the control group (CTL), which received normal chow for 6 weeks. Three different sequencing platforms—snRNA-seq, bulk RNA-seq, and 10X Visium spatial gene expression profiling (Figure 1A, Methods)—were used to analyze the gene expression patterns.
FIGURE 1
ST reveal global demyelination and region-specific transcriptomic changes induced by CPZ
We first compared the Visium spatial maps of Mbp, Cd68, and Gfap genes at the spot level with IHC staining of proteins and showed that the expression patterns were comparable between the Visium and IHC data (Supplementary Figure S1). These genes were chosen based on their well-established expression patterns during de- and remyelination in the CPZ model. The ST dataset contained a total of 57,217 spots (average = 2,043) from all the sections analyzed, with an average of 75,141 reads per spot and a median of 4,157 genes per spot. To look at the overall structure of the ST dataset, we performed unsupervised spot-level pseudobulk principal component analysis (PCA) by aggregating spot counts for each gene across all spots within each section to form a gene-by-section count matrix. Our results showed that the three experimental groups were distinctly separated along PC1 (Figure 1B). We also observed concordance of cell abundance estimates in the same region between animals after ST cell type deconvolution (Supplementary Figure S2), which further suggested that the separation shown by PCA was driven by experimental groups.
We performed unbiased clustering of the ST spots using the SpaGCN package () and identified 17 spatial clusters closely resembling the reference anatomical regions (http://www.brain-map.org) (Figure 1C). The identification of clusters was based on spatially variable genes enriched in anatomical domains and enabled by coherent expression and histology mapping. We identified 484, 1495, 1086 and 2762 differentially expressed genes (DEGs) between CPZ and CTL groups in Cluster 4 (occupying cortical upper layers, roughly layers I-III), Cluster 0 (occupying cortical deeper layers, roughly layers IV-VI), Cluster 8 (roughly occupying white matter tract), and Cluster 5 (roughly occupying hippocampus), respectively. Cluster 5 (hippocampus) had the greatest number of DEGs and highest fold changes among the 4 clusters (Figure 1D). We curated a list of genes that represent MS gene signature (Method) (; ) and found the upregulation of proinflammatory genes (B2m, Cxcl10, C4b, and Trem2) and the downregulation of myelin-related genes (Mbp and Mog) in all 4 clusters (Figure 1D). While the CC has been viewed as the major anatomical location for evaluating demyelination and remyelination in the CPZ model (Xie, 2010; Praet et al., 2014), our analysis showed that inflammatory and demyelination gene signatures were prominent and widespread in various brain regions. Further analysis of biological pathways revealed that neuroinflammation and cell death pathways were prominently upregulated in all 4 clusters/regions. PTEN signaling and the superpathway of cholesterol biosynthesis were downregulated in all 4 clusters, and they have both been linked to myelin integrity and demyelination (; ) (Figure 1E).
Clusters 4 and 0 occupied different layers of the cortex. Higher number of genes from the MS gene signature reached statistical significances in Cluster 0 compared with Cluster 4, suggesting a higher inflammatory state in deeper cortical layers. This observation parallels a previous finding that the upper cortical layers remyelinated more robustly compared with the deeper layers in the CPZ model, and the difference in the repair capacity was attributed to a hypothesized increase in inflammatory activity, myelin debris, and gliosis in the deeper layers (Orthmann-Murphy et al., 2020).
Oligodendroglia and microglia are the major cell types associated with the CPZ model
To resolve cell types at each ST spot, we dissected samples from dorsal brain regions, including the cortex, CC, and hippocampus, to perform snRNA-seq. We obtained a total of 294,447 nuclei with a mean of 3,392 detected genes per nucleus after filtering out low-quality cells and potential doublets and removing ambient RNAs. By transferring cell labels from publicly available datasets, we identified a total of 32 annotated cell types and validated cell purity (Figure 2A; Supplementary Figure S8; Supplementary Table S3) (Marques et al., 2016; ; Scheyltjens et al., 2022).
FIGURE 2
To identify cell types that are associated with treatment conditions, we took two complementary approaches. First, we integrated both bulk and single-cell data using the Scissor package to unbiasedly identify cell populations associated with treatment conditions by asking which single-cell expression profiles best explain the phenotype differences based on the bulk data labels of CTL, CPZ, and RCV (Sun et al., 2022). In the CPZ and CTL contrast, we found that 2.99% of cells (8,802 cells among 294,447 nuclei) were positively associated with the CPZ phenotype (ScissorCPZ+), 14.9% of cells (43,891 cells) were associated with the CTL phenotype (ScissorCPZ−), and the remaining cells were considered background cells (Figures 2B,C). Among ScissorCPZ+ cells, 95.9% were microglia and within the microglial population, 83.3% were classified as ScissorCPZ+. This result suggests that microglia are highly correlated with CPZ phenotype, if not the driving force of the phenotype. On the other hand, the cell types were more diverse among ScissorCPZ− cell. We found 13.1% were MOL, 2.14% were MFOL (myelin-forming oligodendrocytes), 21.2% were EN_Hpc, 20.6% were EN_L23, 19.4% were EN_L5 and 6.9% were EN. Furthermore, 70.8% of MFOL, 40.6% of EN_Hpc, 22% of MOL, 23.8% of EN_L23, 21.9% of EN_L5 and 18% of EN were classified as ScissorCPZ− within their respective cell type populations (Figure 2C). In the CPZ and RCV contrast, 4.13% of the cells (12,182 cells) were positively associated with the RCV phenotype (ScissorRCV+ cells), and 1.95% of the cells (5,729 cells) were associated with CPZ phenotype (ScissorRCV− cells). The ScissorRCV+ cells comprised 25% MOL, 48% EN_Hpc, and 7.8% MFOL, and the ScissorRCV− cells comprised 96.6% microglia. This result suggests that oligodendroglial populations were closely associated with the RCV phenotype. While some microglia had adopted the RCV phenotype, there was still a considerable microglia population that remained associated with the CPZ phenotype. Among the Scissor cell types, microglia had the greatest number of DEGs, and neurons generally had lower numbers of DEGs (Supplementary Figure S3).
Secondly, we performed differential expression analysis using NEBULA across all annotated cell types in our snRNA-seq dataset (22 cell types passing quality filters) for all three contrasts (CPZ vs. Control, Recovery vs. Control, CPZ vs. Recovery), using FDR <0.05 and |log2FC| > 0.25 as significance thresholds (Supplementary Figure S7; Supplementary Table S3). In the CPZ vs. Control contrast, MFOL showed the highest DEG count, followed by MOL and microglia despite a small population of MFOL (n = 1,331). In the CPZ vs. RCV contrast, MFOL and MOL had highest numbers of changes (2,085 and 1,125 DEGs, respectively), while microglia showed a more moderate change (418 DEGs), consistent with our finding that microglia maintained a sustained demyelination-associated phenotype during recovery (Figures 2C, 3G).
FIGURE 3
Since CPZ selectively targets the vulnerable oligodendrocyte population, it is reassuring to see them identified by both methods. Microglial involvement in the model is also consistent with previous publications using histology or omics methods (; ; ). The results of Scissor further suggests that microglia exhibited a sustained demyelination phenotype, while oligodendroglia exhibited a phenotype closer to that of the CTL group after 2 weeks of remyelination. The similarity of oligodendroglia between the RCV and CTL phenotypes was further supported by pseudotime analysis. That is, the CPZ phenotype was distant from the RCV and CTL phenotypes, while the RCV and CTL phenotypes were indistinguishable from each other (Figure 2D). While we recognize some technical limitations of Scissor, the convergence of phenotype association by Scissor and the direct differential expression analysis by NEBULA provides orthogonal support for our conclusion that oligodendroglia and microglia are the principal cell types affected in the CPZ model.
Identification of oligodendroglial and microglial subclusters associated with CPZ phenotype
Following the Scissor finding, we focused on oligodendroglia and microglia in the snRNA-seq dataset. We identified 12 subclusters from oligodendroglia, including subclusters 0, 1, 3, 4, 6, 9, 10, and 11 from MOLs; subclusters 3, 4, and 6 from MFOLs; subclusters 5 and 7 from NFOLs (newly formed oligodendrocytes); subcluster 5 from COPs (differentiation-committed oligodendrocyte precursors); and subclusters 2 and 8 from OPCs (oligodendrocyte precursors) (Figure 3A) using the Leiden clustering algorithm (see Material and Method). Based on T-test, the proportion of cells in subcluster 3 increased significantly in the CPZ condition (CPZ vs. CTL, p value = 6.223e-10) and decreased significantly in the RCV condition (RCV vs. CPZ, p value = 5.861e-10) (Figures 3B,C). This result supports that subcluster 3 is exclusively associated with the CPZ condition and expresses the marker genes Arap2, Dock10, Tenm4, Pex5l and Dock1 (Figure 3D). Most of these genes are involved in small GTPase signaling and cytoskeleton remodeling (Thurnherr et al., 2006), which suggests extensive oligodendrocyte structural remodeling associated with the CPZ condition. Mutations of some of these marker genes have been linked with dysregulation of myelination (; ). On the other hand, subclusters 0, 1, and 4 showed a marked reduction in proportion in the CPZ condition (CPZ vs. CTL, p value = 0.006, 0.0004, and 0.0475, respectively). The proportions of these 3 subclusters increased in the RCV condition, suggesting a phenotype reversal, i.e., remyelination (RCV vs. CPZ, p value = 0.03, 0.001, and 0.002, respectively) (Figure 3C).
We also identified 12 subclusters from microglia (Figures 3E–H). More than 88% of the CTL microglia were assigned to subcluster 0, which expressed homeostatic genes, such as Tmem119, P2ry12, Tgfbr1, and Siglech. Microglia expanded substantially during demyelination, and the proportions of microglia in the CTL, CPZ and RCV groups were 19%, 46.5% and 34.5%, respectively. CPZ-associated subcluster 1 showed high expression of the ApoE, Axl, Cd9 and Lpl genes, which have been implicated as neurodegenerative disease-associated genes (). The proportion of subclusters 1 and 5 exhibited a significant increase during demyelination (CPZ vs. CTL, p value = 1.567e-08 and 0.029, respectively), followed by a significant decrease during remyelination (RCV vs. CPZ, p value = 0.009 and 0.035, respectively) (Figure 3G). Unlike those in oligodendroglia, the changes in the proportions of microglial subclusters from the CPZ to the RCV condition are far from normalization compared to those in the CTL condition. Surprisingly, we did not find any subclusters of astrocytes associated with the CPZ phenotype (Figure 3I).
CellChat identified phagocytosis, growth factor, and chemokine pathways as novel cell-cell interactions from the CPZ condition in the snRNA-seq dataset
We analyzed the snRNA-seq dataset by CellChat () to uncover novel enriched ligand-receptor (LR) pairs. We were particularly interested in “novel” LR pairs that appeared and were significant in the CPZ condition but either did not exist or were not significant in the CTL condition.
We identified enhanced phagocytosis, growth factor, and chemokine pathways in microglia, either as receiver or emitter cells under the CPZ condition (Figures 4A,B). When microglia were analyzed as receiver cells, they paired with multiple emitter cell types with significant novel LR pairs, such as GAS (Gas6-Mertk), NRG (Nrg1-(Itgav + Itgb3)), CX3C (Cx3cl1-Cx3cr1), VEGF (Vegfa-Vegfr1), ANGPT (Angpt2-(Itga5 + Itgb1)), and ANGPTL (Angptl2-(Itga5 + Itgb1)) (Figure 4A). The Gas6-Mertk pathway is considered an anti-inflammatory pathway and has been shown to be critical for modulating microglial phagocytic activity to remove myelin debris, therefore promoting remyelination in the CPZ model (Shen et al., 2021). Interestingly, pathways involved in blood vessel remodeling, such as the VEGF and ANGPT pathways, were predicted to be present in microglia (Proescholdt et al., 2002; Thurston and Daly, 2012). Increasing evidence has suggested that interactions between microglia and blood vessels could play roles in MS and other diseases (; Zhao et al., 2018). Other novel LR pairs in the CPZ condition included FGF (Fgf1-Fgfr2, astrocytes to microglia), GDN (Gdf15-Tgfbr2, MOL to microglia), and BMP (Bmp7-(Acvr1 + Bmpr2), OPC to microglia) pairs. Growth differentiation factor 15 (GDF15) is a neurotrophic factor of the TGFβ superfamily, and elevated levels of GDF15 in serum are associated with MS disease stability and are presumed to be anti-inflammatory (). Gdf15 was recently identified as a marker of demyelination-associated oligodendrocytes and was specific to the CPZ model but not the 5XFAD model of Alzheimer’s disease (). Here, we further predicted that Gdf15 expression in MOLs may exert its function through Tgfbr2 on microglia.
FIGURE 4
When microglia were analyzed as emitter cells, interactions between IGF (Igf1-Igf1r), GRN (Grn-Sort1), FGF (Fgf2-Fgfr2), VISFATIN (Nampt-Insr), and PROS (Pros1-Tyro3) were predicted in multiple receiver cell types (Figure 4B). PSAP (Psap-Gpr37l1, microglia to astrocytes), EGF (Tgfa-Egfr, microglia to astrocytes), PDGF (Pdgfb-Pdgfra, microglia to OPC), and IL1 (Il1a-(Il1r1 + Il1rap), microglia to EN-Hpc) pathway interactions were predicted between microglia and selective cell types. Finally, TGFb (Tgfb1-(Tgfbr1 + Tgfbr2)), ApoE (ApoE-Trem2) and Csf1-Csf1r (CSF) were identified as pathway interactions between microglia. Among these pathways, most of them have been demonstrated to be neuroprotective, while the CSF1-CSF1R interaction is considered proinflammatory. For example, the IGF pathway could play a role in decreasing the severity of neuroinflammation (), and Gpr37l1 expression in astrocytes modulates neuronal N-methyl-D-aspartate (NMDA) receptor activation during ischemia to exert protective effects (). So far, LR analysis by CellChat revealed that microglia interacting pairs may provide protective functions in addition to proinflammatory functions under CPZ conditions.
Most of the LR pairs identified from MOL were involved in growth factor signaling pathways (Figures 4C,D), such as the FGF, EGF, IGF, and VEGF pathways, regardless of whether the MOL was a receiver or emitter cell type. Interestingly, our analysis also revealed MOL and OPC interactions under the CPZ condition. When MOLs were analyzed as receiver cells, they received FGF, PTN and NRG signaling from OPCs, whereas OPCs received PDGF signaling from MOLs (Pdgfc-Pdgfra, OPC to MOL) when the MOLs were analyzed as emitter cells. To our knowledge, the cell-cell interaction between MOLs and OPCs is novel and could be important for understanding the repair capacity and/or maintenance of homeostasis of oligodendrocytes in diseases.
The CellChat analysis revealed that while microglia have a high likelihood of interacting with other microglia and astrocytes, they also interacted with MOLs and OPCs. Conversely, MOLs have a high likelihood of interacting with microglia. The data so far suggests that not only could oligodendroglia and microglia play major roles during de- and remyelination, but they may also interact directly with each other throughout the process ().
ST maps CPZ-associated cell subclusters to CC
Since the 10X Visium Spatial Gene Expression platform captures multiple cells per spot, we used the snRNA-seq data as a reference to deconvolute cell types at each spot by cell2location (Figure 5A) (). As expected, EN-Hpc neurons were mapped to the hippocampal region after deconvolution (Figure 5A). To take advantage of ST to spatially resolve gene expression in each cell type, we manually annotated regions of interest (the CC, cortex, and hippocampus) on each Visium section, using adjacent histology section images as references (Figure 5A).
FIGURE 5
CPZ-associated microglial subclusters 1 and 5 were mapped mostly to the medial CC on the representative CPZ section, and these subclusters were almost nonexistent on the CTL or RCV sections (Figure 5B; Supplementary Figure S6). This finding is consistent with the IHC data showing that microglial activation is most prominent in the CC under CPZ (Supplementary Figure S1). In contrast, homeostatic microglia subcluster 0 did not concentrate in any brain region. CPZ-associated oligodendroglia subcluster 3 was mapped mostly onto the lateral CC on the representative CPZ section and was almost nonexistent on the CTL or RCV sections. Oligodendroglia subcluster 0 was mapped to the entire CC region in CTL and to the medial CC region in the RCV section. Interestingly, MFOLs from the oligodendrocyte lineage mostly mapped to the RCV section (Supplementary Figure S4A), suggesting an active remyelination process.
We next investigated the abundance of each cell type by region under each condition. In general, there was a substantial increase in the number of cells per Visium spot in the CC (Figure 5C) compared with other brain regions. Microglia exhibited a sharp increase in cell abundance in the CC and hippocampus in the CPZ group. The abundance of microglia remained high in the RCV group, suggesting incomplete reversal of the phenotype to that of the CTL group. On the other hand, the MOL abundance in the CC, cortex, and hippocampus was reduced in CPZ but recovered to the CTL level in RCV, suggesting a robust remyelination process after CPZ withdrawal. The OPC number in the CC continued to remain high under the RCV condition.
Finally, we investigated region-specific cellular responses to CPZ. Microglial populations in the cortex and CC had their own gene signatures under control conditions. Under the CPZ condition, both cortical and CC microglia responded by increasing gene expression. In the microglia CC vs. cortex contrast, 11% of genes were unique to the CPZ condition (Figure 5D). MOLs located in the CC and cortex also had different gene signatures in the CTL condition, and both responded mildly to CPZ (Figure 5E). Overall, we found region-specific and cell type-specific responses under the CPZ condition, which is further supported by astrocytes (Supplementary Figure S4B) and OPC populations. OPCs located in the CC and cortex showed similar gene expression signatures, and they did not seem to respond to CPZ by changing their expression even though the number of OPCs markedly increased under CPZ conditions (Figure 5C; Supplementary Figure S4C).
SpaTalk predicts heterogeneous cell‒cell communication in the cortex and CC
We used SpaTalk to infer spatially resolved cell‒cell interactions from the deconvoluted Visium data (Shao et al., 2022) and focused on interactions among 4 cell types: microglia, astrocytes, MOL, and OPC. In the cortex, the interaction between astrocytes and microglia in the CPZ had the greatest number of identified receptor pathways, including death receptor, TNF, and interleukin family signaling pathways. The interaction between microglia and MOL in the CPZ condition was associated with the second highest number of receptor pathways identified, including interleukin family signaling, ERK1/2 activation, and the Wnt signaling pathway. Most receptor pathways in the cortex were identified in the CPZ group, while the MOL-to-OPC and OPC-to-MOL-enriched receptor pathways were identified in the CTL group (Supplementary Table S1).
In the CC, fewer significant receptor pathways were identified than in the cortex, and astrocyte-to-microglia interactions continued to be associated with the CPZ condition (Supplementary Table S2). Selected receptor pathways of interest from SpaTalk with a representative LR pair are plotted in Supplementary Figure S5. Finally, we applied the same filters to CellChat and SpaTalk and identified 20 consensus LR pairs in the CC under the CPZ condition (Table 1). Growth factor and phagocytic pathways were enriched in 11 out of the 20 consensus LR pairs, and microglia were the major sender cell type.
TABLE 1
| Sender | Receiver | Ligand | Receptor |
|---|---|---|---|
| Astrocyte | Microglia | Fgf1 | Fgfr2 |
| Astrocyte | Microglia | Gas6 | Axl |
| Astrocyte | Microglia | Pros1 | Axl |
| Astrocyte | Microglia | Vegfa | Flt1 |
| Astrocyte | Microglia | Vegfb | Flt1 |
| Microglia | Astrocyte | Fgf2 | Fgfr1 |
| Microglia | Astrocyte | Fgf2 | Fgfr2 |
| Microglia | Astrocyte | Fgf2 | Fgfr3 |
| Microglia | Astrocyte | Gas6 | Axl |
| Microglia | Astrocyte | Igf1 | Igf1r |
| Microglia | Astrocyte | Pdgfa | Pdgfrb |
| Microglia | mol | Fgf2 | Fgfr2 |
| Microglia | mol | Igf1 | Igf1r |
| Microglia | opc | Fgf2 | Fgrf2 |
| Microglia | opc | Igf1 | Igf1r |
| Microglia | opc | Pdgfa | Pdgfra |
| mol | Microglia | Gdf15 | Tgfbr2 |
| mol | Microglia | Vegfa | Flt1 |
Common LR pairs in CC under the CPZ condition predicted by both CellChat and SpaTalk.
The CPZ model exhibits moderate translatability to MS chronic active lesions
Lastly, we explored the translatability of the CPZ model by performing an integrated analysis of our dataset with a published MS RNA-seq dataset that profiled chronic active lesions in humans (). MS white matter lesions can be classified into various types based on the infiltrating immune cells, the activity of monocyte-derived macrophages and resident microglia, and the degree of demyelination (). Among white matter lesions, disease disability correlates with chronic active lesions that have iron-laden microglia present at the lesion rim, which can be visualized by specific magnetic resonance imaging (MRI) sequences in the clinic (; ). We reanalyzed both human and mouse datasets using the same methods presented in this study and calculated DEGs between the disease and control groups to determine the level of DEG conservation between the two datasets. To take advantage of the lesion segmentation defined by Absinta et al., we calculated DEGs of both “lesion core vs. control” and “chronic active lesion edge vs. control” for each glial cell type (Figure 6, left). We then calculated the overlap in DEGs between mice and humans and subsequently performed KEGG analysis on DEGs common to both species. We found pathways that were conserved between the CPZ mouse and human MS samples depending on both the cell type and the lesion architecture investigated. For example, gene expression changes in microglia in the CPZ condition were more in concordance with gene expression changes in the MS lesion core, whereas changes in oligodendrocyte in CPZ were more in concordance with changes at the lesion edge.
FIGURE 6
We also used an alternative approach to investigate the translatability by using a variation of weighted correlation network analysis designed for single-cell data (hdWGCNA) (Morabito et al., 2021; Morabito et al., 2023). We used hdWGCNA to define disease-associated modules that were both correlated with the disease state and preserved between the two species (). We found disease-associated modules in the CPZ model that showed moderate to strong preservation in humans depending on the cell type investigated (Figure 6, middle column). Astrocytes had the highest preservation scores of disease-associated modules, while MOL, OPC, and microglia showed moderate to low preservation between the CPZ and the Absinta et al. datasets. These disease-associated and preserved modules showed KEGG pathway enrichment that was shared among cell types, such as metabolic pathways (in astrocytes, microglia, and the MOL), protein processing in the endoplasmic reticulum (in astrocytes and microglia), and ubiquitin-mediated proteolysis (in astrocytes, microglia, and the MOL) (Figure 6, right column). In the future, this workflow can be applied to compare different MS lesion types with the CPZ model or others to understand the validity of mouse models.
Discussion
In this study, we profiled the CPZ model by “space” and “time” using ST along with an array of bioanalytical tools. We found that demyelination and neuroinflammation in the CPZ model were more prominent and widespread across the dorsal to ventral axis than previously identified by IHC. We identified oligodendroglia and microglia as two major cell types associated with de- and remyelination in the model and we mapped CPZ-associated subclusters of oligodendroglia and microglia in the CC, which shows stereotypical demyelination. In the RCV group, while MOL almost normalized to the control state, microglia continued to be associated with the CPZ phenotype, i.e., demyelination. We inferred cell‒cell interactions using both snRNA-seq and ST data and found that the phagocytosis, growth factor, and chemokine pathways were enriched under the CPZ condition. ST is a fast-emerging field that enables high-throughput investigations of spatially resolved transcripts (Liu et al., 2021; Rao et al., 2021). Prior to the development of ST technology, region-specific gene expression studies in MS brains relied heavily on manual dissections (; ; ). Manual dissection of lesions is sufficient to identify major biological heterogeneity but lacks fine resolutions. Based on our study of the CPZ model, ST appears to be an appropriate tool for investigating detailed cellular and molecular substrates underlying MS lesion evolution as well as predicting cellular interactions.
Our findings identified oligodendroglia and microglia as cell types of interest during de- and remyelination. This result is consistent with the IHC results showing that demyelination (oligodendrocytes) and gliosis (microglia) are two prominent features of the CPZ model, particularly in CC lesions (Praet et al., 2014). However, considering proposed functions of astrocytes in the experimental autoimmune encephalomyelitis (EAE) model and MS (Ponath et al., 2018; Wheeler et al., 2020; Wheeler et al., 2023), it is unexpected that astrocytes were not significantly associated with any experimental conditions, and no disease-associated subclusters were identified in our dataset. It is possible that (1) the recovery time point (2 weeks after CPZ withdrawal) is too early to capture a full astrocytic response in our CPZ model, or (2) astrocyte involvement in demyelination pathology is animal model dependent, as the EAE model is more inflammatory in nature. (3) It has been suggested that reactive astrocytes are highly heterogeneous and the responses are context/region dependent. Moreover, there is co-existence of many pathology-related reactive sub-states of astrocytes (Patani et al., 2023). Thus, it is possible that the astrocyte population in the CPZ model undergoes a widespread and substantial gene expression change that is moderately consistent with the changes in MS without forming identifiable CPZ-associated subclusters. Nonetheless, CellChat and SpaTalk predicted numerous potential LR pairs between astrocytes and microglia, highlighting the interactions between the 2 cell types and the neuroinflammatory pathways associated with these interactions. Interestingly, astrocytes had the highest preservation scores of disease-associated modules with MS among all the cell types analyzed by hdWGCNA for model translatability. Although we did not identify astrocytes as “pathogenic” cells in the CPZ model, this finding suggests that astrocyte transcriptomic changes in the model are still in line with the changes in MS.
The growth factor and phagocytic pathways were predicted by both CellChat and SpaTalk as potential LR interactions in the CC under the CPZ condition. The role of growth factors in MS has long been recognized (Webster, 1997), and approximately three-quarters of the LR pairs in our study are related to growth factor signaling. While experimental validation of the identified LR interactions is out of scope in this study, we propose growth factor pathways as a high priority direction for future functional investigations, particularly the VEGF and FGFR pathways. In addition to blood vessel remodeling, the VEGF pathway is involved in pathogenic glial responses (Rothhammer et al., 2018) and meningeal lymphatic dysfunction () that are relevant to MS. Based on our findings, further investigation of the contributions of VEGFA-FLT1 and VEGFB-FLT1 cross-talks across space and time will help elucidate the disease mechanisms at the CNS borders. Through FGFR signaling, FGFs are involved in various biological and pathological processes, and FGF1 and FGF2 are the main mitogens in the CNS (). Evidence has shown that FGF1 and FGF2 are expressed in MS lesions with a remyelination phenotype (; Mohan et al., 2014), and it has also been shown that FGF signaling is critical for regulating oligodendrocyte and myelin regeneration (). Thus, recombinant FGFs have been proposed to be a potential therapeutic approach for MS. On the other hand, oligodendrocyte-specific deletion of FGFR1 and FGFR2 results in a protective effect in the EAE model; therefore, FGFR inhibitors have also been proposed for MS (Rajendran et al., 2023). Disentangling the complex biology of FGF signaling will be critical for MS therapy development and utilizing an ST approach can provide much needed information in terms of cell-cell interaction in situ.
The remaining LR pairs are related to phagocytic pathway signaling. Myelin debris have been shown to inhibit remyelination (), and it has also been suggested that a dysfunctional myelin sheath under inflammatory attack could pose a greater risk of axonal degeneration (Schaffner et al., 2023). Therefore, a precise and efficient phagocytic program to clean myelin debris is critical for successful remyelination. In nondemented individuals from the Swedish BioFINDER-2 cohort, increased levels of microglial markers in the CSF, such as soluble TREM2, AXL, MERTK, GAS6, LPL, CST7, SPP1 and CSF1, are associated with slower cognitive decline (Pereira et al., 2022). In our dataset, Gas6-Axl and Pros1-Axl were in the consensus LR pair list from CellChat and SpaTalk, while Gas6-Mertk and ApoE-Trem2 were predicted from the CellChat analysis. Gas6 and Pros1 are ligands for the TAM (Tyro3, Axl, and Mertk) receptor tyrosine kinase family, and their signaling regulates inflammation, cell growth, and clearance of debris (). The role of phagocytosis-related proteins in MS and whether they can be potential therapeutic targets are less clear than those in Alzheimer’s disease research (Weinger et al., 2009). However, studies in the CPZ model have suggested that both Mertk and Trem2 play roles in myelin debris clearance and subsequently promote remyelination (; Shen et al., 2021). A detailed profiling of phagocytic proteins in MS cerebrospinal fluid to see how they correlate with MS disease severity and progression could be a worthwhile effort. Overall, the LR pairs prediction in our study showed similarities to cell‒cell interactions identified from human MS cortical ST data (). For example, growth factor pathway interactions, such as PSAP-GPR37L1, CSF1-CSF1R, and PDGFC-PDGFRA, which represent approximately one-third of the pathways identified in the MS cortical ST data, were also present in our CPZ dataset. One-fourth of the pathways in the MS cortical ST data were grouped as anti-inflammatory, such as the PROS1-AXL interaction, which was also predicted in our dataset. These findings further support the moderate translatability of CPZ model to human MS.
The advantage of our study is that we employed a hypothesis-free whole-genome approach, but we also encountered some study limitations. First, while cell‒cell interaction investigations in our study provided novel findings, the LR pairs did not immediately predict therapeutic directionalities. Due to the scope of the study, functional validation will be performed in a separate investigation. Second, it is important to acknowledge that 10X Visium technology lacks single-cell resolution, which is a key limitation. To address this issue, we used cell2location, which integrates snRNA-seq data with spatial data to deconvolute cell types, providing insights into the abundances of different cell types within the 10X Visium spots. Although this technique does not achieve true single-cell resolution, it represents a crucial step toward a more detailed characterization of the cellular landscape. Third, the sample size in our study is small. Some potential results may not reach statistical significance due to the small sample size. Conversely, significant findings from our study should also be validated in a larger study or using orthogonal methods. Finally, while we believe that the CPZ model captures some of the salient demyelination features of MS, we cannot exclude direct effects of CPZ on cells other than the oligodendrocyte population or effects from CPZ on biological processes unrelated to demyelination, which could introduce noise into the analyses (Skripuletz et al., 2011).
Materials and methods
Cuprizone model
All procedures involving animals were reviewed and approved by the Biogen Institutional Animal Care and Use Committee and conducted under protocol 801 in accordance with relevant guidelines and regulations. All methods are reported in accordance with ARRIVE guidelines (https://arriveguidelines.org). Male wild-type C57BL/6J mice aged about 10 weeks were purchased from Charles River Laboratories (stock number 027C57BL/6). The research animals were housed under specific-pathogen-free (SPF) conditions in an AAALAC accredited facility according to Biogen’s institutional animal care and use committee (IACUC) protocol with a 12 h–12 h light–dark cycle and environmental conditions controlled at 72°F/22 °C and 40%–60% humidity. To establish the CPZ model, customized 0.3% (w/w) cuprizone chow was made by mixing cuprizone powder (Sigma‒Aldrich 14690) with control base chow (AIN-76A) at Research Diets Inc. All food chow (CPZ and control) was stored in vacuum bags at 4 °C, and fresh food chow was replaced weekly during the study. In the CPZ group, the animals were fed cuprizone chow for 4 weeks. For the recovery group, animals were fed cuprizone chow for 4 weeks, followed by 2 weeks of control base chow for recovery. For the control group, animals were provided with control base chow for 6 weeks. The body weight of each mouse was closely monitored weekly throughout the study, and mice with ≥30% weight loss were excluded.
Tissue collection for sequencing studies
For snRNA-seq, spatial transcriptional sequencing, and bulk RNA-seq studies, the CPZ-treated mice were sacrificed, and brain tissues were collected at either week 4 (n = 3, CPZ group) or week 6 (n = 3, recovery group). The control mouse brain tissues were collected at week 6 (n = 4, control group). At each time point, the mice were euthanized by CO2 inhalation, followed by transcardiac PBS perfusion and brain tissue dissection. Each brain was cut sagittally to separate the left and right hemispheres; one brain hemisphere was trimmed and embedded in OCT for snRNA-seq or spatial transcriptomic study, and the other hemisphere was snap frozen for bulk RNA-seq analysis (Figure 1). To enrich for cells from the corpus collosum and cortex, only the top half of the coronal section, which included the dorsal cortex, CC, and hippocampus, was collected for snRNA-seq.
Immunohistochemistry (IHC)
A separate cohort of animals was used for the IHC study. Brain tissues were dissected and fixed for 24 h in 10% NBF at 4 °C and then transferred to ice-cold PBS for coronal sectioning on a vibratome. Coronal brain slices of 40 µm thickness from bregma −1 mm to −2.5 mm were collected and stored in 30% sucrose/PBS buffer at −20 °C. For staining, 3 slices spanning 400 µm apart were selected from each animal. For staining, brain slices were thawed in PBS at room temperature and then blocked with blocking solution (0.5% Triton, 5% goat serum, and 0.1% BSA in PBS) for 1 h at room temperature. Primary antibodies were diluted in blocking solution and incubated with the slices overnight at 4 °C. The corresponding secondary antibodies were applied for 1–2 h at room temperature. After incubation with secondary antibodies, the slices were washed and mounted onto glass coverslips with Prolonged Gold Mounting Solution (Thermo Fisher, P36931). The following primary antibodies were used: mouse anti-MBP (BioLegend, 88,408, 1:500), rat anti-CD68 (Bio-Rad, MCA1957, 1:300), and chicken anti-GFAP (Thermo Fisher, PA1-10004, 1:500). The following secondary antibodies were used: goat anti-chicken IgY 647 (Thermo Fisher, A32933, 1:500) and donkey anti-rat IgG488 (Thermo Fisher, A21208, 1:500).
Images were obtained using a Zeiss LSM 710 confocal microscope with a ×20 objective and Z-stack of 10 planes at 1 µm intervals. Individual fields of each slice were tiled into a large image.
Tissue sectioning and preparation for spatial transcriptomics
Frozen sections for ST profiling were collected using a Leica CM1520 cryostat with a Leica Surgipath high profile blade at −23 °C. Ten-micron-thick sections were collected from mouse brain hemispheres embedded in OCT (TissueTek) on Visium Spatial Gene Expression slides (10X Genomics Slide Kit, Ref No. 1000185). Tissue sections were carefully collected to allow accurate positioning within individual capture areas of Visium slides following a randomized collection scheme across multiple Visium slides to help control for (1) slide-to-slide batch effects by ensuring that a section from each treatment group was present on each slide and (2) the potential impact section location on the slide by ensuring that sections from a given treatment group were not in the same position or order on every slide. The frozen slides were processed for spatial transcriptomic sequencing according to the manufacturer’s instructions. Briefly, the sections were heated at 37 °C for 1 min and then fixed in methanol at −20 °C for 30 min. Following fixation, the samples were stained with hematoxylin and eosin. Brightfield images were collected on an Olympus VS120 Virtual Slide Microscope equipped with a VC50 camera using a UPLSAPO ×20 objective. After imaging, tissue permeabilization was carried out using the supplied permeabilization buffer at 37 °C for 6 min. Permeabilization time was selected based on a prior optimization experiment in mouse brain tissues using Visium Spatial Tissue Optimization slides (10X Genomics, Ref. No. 1000193). Following permeabilization, reverse transcription and second strand synthesis reactions were carried out on slides before elution of the resulting cDNA into separate sample PCR tubes according to Visium guidelines. qPCR was performed on 1 µL of eluted cDNA as described in the Visium Spatial Gene Expression Reagent Kits User Guide (Revision-B) and was run on an Applied Biosystems QuantStudio 12K Flex qPCR instrument (Applied Biosystems™ 4471087). Individual sample cDNA PCR amplification cycles were selected and run based on sample cDNA Ct. Total yields of amplified cDNA ranged from 300 to 550 ng. Sequencing library generation was performed using the 10X Genomics Library Construction Kit (PN-1000190) following the manufacturer’s instructions with 13 cycles of PCR based on the cDNA input used. Library quality was measured using a PerkinElmer LabChip GXII with a Hi Sens Lab Chip (Perkin Elmer 760517) and DNA High Sensitivity Reagents (Perkin Elmer CLS760672), which showed an average spatial gene expression library size of 473 bp. Libraries were normalized and pooled before loading at 300 pM on an Illumina NovaSeq 6000 S2 flow cell with paired-end sequencing parameters of Read1-28 bp × i7-10 bp × i5-10 bp × Read2-120 bp.
Bulk RNA sequencing and data processing
Bulk RNA-seq libraries were prepared using the Kapa mRNA HyperPrep Kit (Kapa Biosystems KK8581) according to the manufacturer’s instructions. Briefly, 100 ng of RNA was fragmented for 6 min at 94 °C, followed by RT, second strand synthesis, A-tailing, adaptor ligation, and amplification with 14 cycles of PCR. Libraries were quantified on a LabChip GXII as described above. Libraries were normalized, pooled, and sequenced on an Illumina HiSeq2500 platform, with an average depth of 9 million 50 bp paired-end reads. Run quality metrics were assessed using Illumina’s BaseSpace run summary tool. Bulk RNA analysis was performed using the RNASequest pipeline (Zhu et al., 2023) following the tutorial provided in the manuscript.
Bulk RNA-seq integration with scRNA-seq
We applied Scissor (Sun et al., 2022) (https://github.com/sunduanchen/Scissor) to integrate the bulk and snRNA-Seq datasets. It utilizes the phenotypes (CPZ, RCV, and control) collected from bulk assays to identify the most highly phenotype-associated cell subpopulations from single-cell data.
Nuclei isolation
Nuclei isolation was performed as described previously (). Tissue was dounced in 0.5 mL Hyman Lysis Buffer (0.32 M Sucrose, 5 mM CaCl2, 3.0 mM MgAc2, 0.1 mM EDTA, 10 mM Tris-HCl, pH 8.0, 1 mM DTT, 0.1% Triton X-100, 1/1000th volume Promega RNasin Plus RNase Inhibitor) 10X with loose pestle, followed by 10X with tight pestle. Following douncing, 1 mL of Hyman Low Sucrose Buffer (0.32 M Sucrose, 5 mM CaCl2, 3.0 mM MgAc2, 0.1 mM EDTA, 10 mM Tris-HCl, pH 8.0, 1 mM DTT, 1/1000th volume Promega RNasin Plus RNase Inhibitor) was added to the homogenate, and passed through a 70 μm filter, followed by a 1 mL Hyman Low Sucrose Buffer rinse of the dounce and filter. Homogenates were spun at 600 × g for 5–10 min at 4 °C. After supernatant aspiration, pellets were gently resuspended in 5 mL Complete Buffer HB (0.25 M sucrose, 25 mM KCl, 5 mM MgCl2, 20 mM Tricine-KOH, pH 7.8, 1 mM DTT, 0.15 mM spermine, 0.5 mM spermidine, 1/500th volume Promega RNasin Plus RNase Inhibitor), followed by addition of 5 mL Working Solution (50% Iodixanol: 5 volumes OptiPrep Density Gradient Medium, 1 volume Diluent Buffer (150 mM KCl, 30 mM MgCl2, 120 mM Tricine-KOH, pH 7.8), 1/1,000th volume Promega RNasin Plus RNase Inhibitor). In a 12 mL Corex tube, 1 mL 40% iodixanol (Working Solution diluted in Complete Buffer HB, 0.042% BSA), and 30% iodixanol (Working Solution diluted in Complete Buffer HB, 0.021% BSA) were layered, followed by 10 mL of resuspended sample. The samples were spun in an ultracentrifuge at 10,000 × g for 18 min at 4 °C with the brake off. After spinning, the fatty myelin top layers were vacuum aspirated, followed by a 25% iodixanol layer and approximately 1/3 of the 30% iodixanol layer. Approximately 0.5 mL of nuclear sample band at the 30%–40% iodixanol interface was collected, transferred to an Eppendorf tube and mixed by pipetting. Nuclei were stained with DAPI and counted with a Countess II FL. Nuclei suspensions were diluted to 1000 nuclei/μL in RB (PBS + 2% BSA + 1/1000th volume Promega RNasin Plus RNase Inhibitor).
snRNA sequencing
Single-nuclei RNA-seq was performed using the 10X Genomics Chromium Next GEM Single Cell 3′ Kit v3.1 (PN1000268), Chromium Next GEM Chip G (PN1000127), and Dual Index Kit TT Set A (PN1000215). Nuclei suspensions were loaded on the 10X Chromium instrument for recovery of 5000 nuclei per sample. cDNA synthesis, PCR, and library preparation were performed according to the manufacturer’s instructions, with 12 cycles of PCR for cDNA amplification, and 15 cycles for library amplification. cDNA and libraries were quantified on a LabChip GXII as described above. Libraries were normalized, pooled, and loaded at 300pM on an Illumina NovaSeq 6000 S4 flowcell with paired end sequencing parameters of Read1-28bp × i7-10bp × i5-10bp × Read2-90bp targeting an average depth of 85,000 reads per cell.
snRNA-seq data preprocessing and quality control
FASTQ files were generated from Illumina bcl files using CellRanger mkfastq. CellRanger count was used to identify UMI and cell barcodes and to align reads to the mouse genome assembly mm10, with intronic counts included. To remove ambient RNA from the gene expression matrices, the ‘remove-background’ function of CellBender was used (; ). In addition, CellBender identifies and removes empty droplets. The CellRanger output file ‘raw_feature_bc_matrix.h5’ was used as the input file for CellBender. The ‘Estimated Number of Cells’ from CellRanger was used as the ‘expected-cells’ parameter, while a value in the plateau of the barcode rank plot was used as the ‘total-droplets-included’ parameter. Potential doublets were identified and removed using scDblFinder (). Filtered barcode matrices were concatenated using Scanpy (Wolf et al., 2018). QC metrics, including the number of detected genes, total counts, and percentage of counts from mitochondrial genes, were calculated for each barcode. Cell barcodes with fewer than 1000 detected genes, greater than 100,000 counts, or greater than 10% mitochondrial UMIs were filtered out from the dataset.
As an additional QC metric, intronic mapping rates were derived on a per-barcode basis using the BAM alignment files output by CellRanger count. Compared with cellular transcriptomes, nuclear transcriptomes have been demonstrated to show a greater proportion of reads mapping to intronic regions (), enabling the detection of nonnuclear cell barcodes representing contaminating cytoplasmic transcripts.
snRNAseq cell type annotation and subclustering
Cell type labels were assigned using label transfer from E-MTAB-11115 (), with region-specific astrocyte subtypes merged into a single astrocyte label. Oligodendrocyte-lineage cell types were further clarified by label transfer from GSE75330 (Marques et al., 2016) after merging major cell types from this reference. Filtered UMI counts for E-MTAB-11115 were downloaded from ArrayExpress along with cell type annotations. The UMI was then processed by Seurat (), LIGER (Welch et al., 2019) and Harmony () integration on different sample libraries. The kBET () evaluation revealed that LIGER achieved the best performance. Thus, the Seurat SCT and LIGER layout (PCA) was used to generate the Seurat azimuth reference. In parallel, the cells were sequentially subclustered via Harmony and Louvain clustering.
To investigate the heterogeneity within the microglial and oligodendrocyte populations, we subclustered these cell types using the Leiden clustering algorithm implemented in Scanpy. By applying a resolution parameter of 0.5, we identified 12 subclusters each for microglia and oligodendrocytes. Differential gene expression analysis was performed to identify marker genes that exhibited significant differences in expression between subclusters with the Wilcoxon rank-sum test. To visualize marker gene expression, we generated dot plots that depict the percentage of cells expressing the gene through dot size, and the color intensity represents the relative expression level.
snRNA-seq DE analyses
snRNA-seq data were analyzed per cell type using the NEBULA-HL single cell DE method (). For each cell type, 3 contrasts were considered: CPZ vs. Control, CPZ vs. Recovery, and Recovery vs. Control. Cell types with at least 10 cells in each of the biological replicates within the contrast were analyzed. The genes tested in each contrast group were expressed in at least 10% of the cells in both contrast groups. Mouse ID was used as a random effect in the model. Library size was used as an offset covariate.
Pseudotime trajectory analysis and velocity visualization
We used scFates () for trajectory analysis. scFates combines a tree inference approach similar to SimplePPT (Mao et al., 2015) and a principal graph learning algorithm named ElPiGraph () to calculate trajectories.
10x Visium spatial sequencing data analysis
The 10x Visium spatial sequencing data were processed using Space Ranger v1.1.0 (10x Genomics). The sequencing results (BCL files) were converted to fastq files with the Space Ranger mkfastq pipeline. The fastq files were then mapped to the mouse reference genome (GRCm38.p3) using the Space Ranger count pipeline, which produces a spot-barcoded expression matrix.
Next, we used SpaGCN (), a graph convolutional network-based method, to identify spatial clusters based on both gene expression profiles and spatial locations. All the slides were aligned and normalized using the PASTE algorithm (Zeira et al., 2022). The MS gene signature selected for this study was based on consensus genes in multiple sclerosis signaling, myelination signaling and neuroinflammation signaling pathways from QIAGEN Ingenuity Pathway Analysis (IPA) that were identified in Figure 1E, including: B2m, C1qa, C1qb, C1qc, C3, C3ar1, C4b, Ccl2, Cd86, Creb3l1, Cxcl10, Egr2, Eif4ebp1, Fos, Hck, H2-Q4, H2-Dmb1, H2-Oa, H2-Aa, H2-Ab1, H2-Eb1, Hmox1, Icam1, Idi1, Igf1, Il33, Irf7, Itgax, Itgb2, Itgb4, Itgb5, Mag, Mapk4, Mbp, Mobp, Mog, Mthfd2, Ncf2, Opalin, Parp14, Parp9, Plp1, Pycard, Rac2, Rnf213, Serping1, Slc6a12, Tgfbr1, Tlr2, Tlr3, Tlr7, Trem2, Tyrobp.
Then, we used Cell2location () to spatially map cell types in the Visium spatial transcriptomics data. This package integrates single-cell RNA-seq (scRNA-seq) or single-nucleus RNA-seq (snRNA-seq) with spatial transcriptomics data to deconvolute cell types in each spot of the Visium slide. A reference mouse brain cell type signature representing gene expression profiles for each annotated cell type was estimated from snRNA-seq data using a negative binomial regression model. The parameters used for the model were batch size = 2500, train_size = 1, learning rate = 0.002, and max_epochs = 250. Then, the reference cell type signature was integrated with the processed 10X Visium data to deconvolute the absolute spatial abundance of each cell type on the spatial transcriptomics slide. The deconvolution model requires the user to input two hyperparameters: 1) the expected cell abundance (N_cells_per_location), which was set to 5, and 2) the regularization strength of the detection efficiency effect (detection_alpha), which was set to 200 (default). The model was trained for 30,000 iterations. The default settings were used for the remaining parameters.
ST pseudobulk PCA
To perform pseudobulk PCA of the spatial transcriptomics data, spot counts for each gene were aggregated across all spots within each slice to form a gene-by-slice count matrix. After performing variance stabilizing transformation (vst), the top 500 most variable genes from the normalized matrix were used to generate PCs. The plot shown in Figure 1B was generated from data processed with DESeq2 (Love et al., 2014) after correcting for batch effects using ComBat-seq (Zhang et al., 2020).
ST DE and pathway analysis
ST spot-level data were analyzed per cell type using the NEBULA-HL single-cell DE method (). In this analysis, each spot was handled independently (spot-level information regarding spatial adjacency was not considered). Regional subsets of mouse brain regions, including the CC, cortex, and hippocampus, were manually annotated. H&E histology section adjacent to each ST section were collected and processed as an overlay image to guide the anatomical annotation. Regional subsets of mouse brain regions, including the CC, cortex, and hippocampus, were manually annotated by examining the ST spot images with histology images to include spots in appropriate brain regions. Allen Brain Atlas mouse reference atlas (https://mouse.brain-map.org/) was referenced to confirm the annotation. For each region, 3 contrasts were considered for DE analysis, namely, CPZ vs. Control, CPZ vs. RCV, and RCV vs. Control. Spots with fewer than 100 counts or greater than 45% mitochondrial gene counts were filtered out from each slice before analysis. Mouse ID was used as a random effect in the model. Library size was used as an offset covariate. DEGs from each cluster were analyzed by QIAGEN Ingenuity Pathway Analysis (IPA) for pathway discovery using Core Analysis. To compare pathways across clusters, all 4 clusters were subjected to Comparison Analysis.
Cell‒cell interaction analysis
We applied CellChat () to understand the cell-cell interactions from the snRNAseq dataset. For ligand-receptor (LR) analysis, the scRNA-seq count matrix and corresponding metadata were loaded into a Seurat object for use in CellChat. The count matrix was normalized, and cell type identities were assigned to cells. The mouse specific CellChat LR interaction database ‘Secreted Signaling’ was selected. CellChat with default settings was used to perform cell-cell interaction (CCI) analysis and infer the CCI network by computing cellular and signaling pathway communication probabilities (strength). The method truncatedMean was used to calculate the average gene expression per cell group. Each comparison (CPZ vs. Control, CPZ vs. Recovery, Recovery vs. Control) was analyzed separately and then merged into a single CellChat object. Differential LR combinations within the CCI network were identified in CellChat by comparison of the communication probabilities, and then those results were extracted for the cell type pairs of interest.
We applied SpaTalk (Shao et al., 2022) to infer cell‒cell interactions by enriching receptor pathways from the ST dataset. SpaTalk was run on three spatial sections per condition, with one section per biological replicate. Scaled cell type deconvolution weights precalculated by cell2location were used to deconvolve spots into cell types that were derived from the scRNA-seq reference data. Deconvoluted spots from two annotated regions, the cortex and CC, were analyzed in each section. Cell type assignments per spot were made using the default recommended size for 10X Visium 55-micron spots. Default SpaTalk spot-level analysis parameters were used to calculate cell‒cell LR interactions and to perform enriched LR receptor pathway analysis. Overlapping significant cell‒cell LR interaction pairs from both CellChat analysis based on snRNA-seq data and SpaTalk analysis based on ST data were used for downstream comparisons between CPZ and control conditions. The following set of 8 custom LR interaction pairs relevant to MS were added to the overlapping set: Il34-Csf1r, Apoe-Trem2, Pecam1-Pdgfra, Th-Drd1, Th-Drd2, Ddc-Drd1, Ddc-Drd2, and Spp1-Cd44.
Translatability and hdWGCNA
For comparative DEG analyses between CPZ and the human dataset, we first reprocessed the data from (see snRNA-seq data preprocessing and quality control) and called DEGs between control and MS tissues using Nebula (see snRNA-seq DE Analyses). For the human dataset, we performed DEG analysis on both the “control vs. lesion core” and control vs. “chronic active lesion edge” datasets to take advantage of the different lesion areas provided by Absinta et al. We identified genes that were upregulated in the disease state, with an FDR <0.05 for both datasets. We then converted mouse gene IDs into their human orthologs and counted the number of distinct and overlapping DEGs between the two datasets on a cell type basis. For a given cell type, DEGs that were common between CPZ and humans were used as input to the R package KEGGprofile for pathway analysis (https://github.com/slzhao/KEGGprofile/blob/master/DESCRIPTION).
Single-cell WGCNA was performed according to hdWGCNA (Morabito et al., 2021; Morabito et al., 2023) with modifications to regress out pseudobulk covariates, such as age and sex. The working principle of hdWGCNA is to first create pseudobulk samples of single-nucleus data based on kNN and subsequently identify groups of coexpressed genes (“modules”) for each cell type. Briefly, raw pseudobulk expression was calculated from specific cell populations by aggregating UMI counts from 50 cells using a k-nearest neighbor approach. We then normalized the raw count data by variance stabilizing transformation (VST) implemented by DESeq2. We regressed out confounding features such as total UMI count, mitochondrial percentage, and intronic rate. We then applied WGCNA to identify coexpressed genes (modules) using the normalized count data. Briefly, WGCNA starts by constructing a topological overlap (TO) matrix (Zhang and Horvath, 2005). To this end, we used weighted correlation with individual sample weights determined as described above and the “signed hybrid” network in which negatively correlated genes are considered unconnected. We used module eigengenes (MEs), the first principal component of the module’s gene expression matrix, to represent the overall expression patterns of the corresponding coexpression modules. The disease correlation score was defined as the Pearson correlation coefficient between the module eigengenes and disease status. Finally, we performed module preservation analysis using the modulePreservation function in WGCNA. The level of translatability between species is defined by the module preservation score (). The output consists of “modules” of genes with similar trajectories, including scores to assess both correlation with disease (module correlation) and preservation between mice and humans (Zsummary). We considered modules with a disease correlation score >0.7 and a Zsummary >2 to be associated with disease and preserved between species. For pathway analysis, we used the module as an input for the KEGGProfiler.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/, GSE255371.
Ethics statement
The animal study was approved by Biogen Institutional Animal Care and Use Committee. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
H-HT: Writing – review and editing, Supervision, Writing – original draft, Project administration, Conceptualization. SP: Formal Analysis, Writing – original draft, Writing – review and editing. JW: Writing – review and editing. JZ: Writing – review and editing. WH: Writing – review and editing. ARG: Writing – review and editing, Formal Analysis. ShC: Formal Analysis, Writing – review and editing. AJG: Data curation, Writing – review and editing, Formal Analysis. SJC: Writing – review and editing, Methodology, Data curation. MS: Writing – review and editing, Methodology. JC: Writing – review and editing, Formal Analysis. ZO: Formal Analysis, Writing – review and editing. MR: Writing – review and editing, Formal Analysis. ML: Formal Analysis, Writing – review and editing. WW: Writing – review and editing, Formal Analysis. EZ: Writing – review and editing, Formal Analysis. PC: Writing – review and editing, Methodology. RC: Writing – review and editing, Methodology. EM: Software, Writing – review and editing. WZ: Resources, Writing – review and editing. YJK-W: Methodology, Writing – review and editing. CE: Writing – review and editing, Methodology. LJ: Methodology, Writing – review and editing. HM: Resources, Writing – review and editing. TC: Writing – review and editing, Methodology. JG: Formal Analysis, Writing – review and editing. TR: Methodology, Writing – review and editing. ML: Writing – review and editing, Resources. KL: Writing – review and editing, Methodology, Formal Analysis, Conceptualization. BZ: Writing – review and editing, Supervision, Methodology, Writing – original draft, Conceptualization.
Funding
The author(s) declared that financial support was not received for this work and/or its publication.
Acknowledgments
The authors thank Bing Zhu and Robert Miller for critically reading the manuscript.
Conflict of interest
H-HT, SP, JW, JZ, WH, ARG, ShC, AJG, SJC, MS, JC, WW, EZ, PC, RC, EM, WZ, YJK-W, CE, LJ, HM, TC, JG, TR, KL, and BZ are current and former Biogen Inc. employees.
Author ZO was employed by Data Science, BioInfoRx Inc.
Author MR was employed by PharmaLex Inc.
The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbinf.2026.1832826/full#supplementary-material
References
1
AbsintaM.MaricD.GharagozlooM.GartonT.SmithM. D.JinJ.et al (2021). A lymphocyte-microglia-astrocyte axis in chronic active multiple sclerosis. Nature. 577 (7878), 709–714. 10.1038/s41586-021-03892-7
2
AlberganteL.MirkesE.BacJ.ChenH.MartinA.FaureL.et al (2020). Robust and scalable learning of complex intrinsic dataset geometry via ElPiGraph. Entropy (Basel)223), 296. 10.3390/e22030296
3
AlsemaA. M.WijeringM. H. C.MiedemaA.KotahJ. M.KosterM.RijnsburgerM.et al (2024). Spatially resolved gene signatures of white matter lesion progression in multiple sclerosis. Nat. Neurosci.27 (12), 2341–2353. 10.1038/s41593-024-01765-6
4
AmstadA.CorayM.FrickC.BarroC.OechteringJ.AmannM.et al (2020). Growth differentiation factor 15 is increased in stable MS. Neurol. Neuroimmunol. Neuroinflamm7 (2), e675. 10.1212/NXI.0000000000000675
5
BaakliniC. S.HoM. F. S.LangeT.HammondB. P.PandaS. P.ZirngiblM.et al (2023). Microglia promote remyelination independent of their role in clearing myelin debris. Cell Rep.42, 113574. 10.1016/j.celrep.2023.113574
6
BakkenT. E.HodgeR. D.MillerJ. A.YaoZ.NguyenT. N.AevermannB.et al (2018). Single-nucleus and single-cell transcriptomes compared in matched cortical cell types. PLoS One13, e0209648. 10.1371/journal.pone.0209648
7
BarkhofF.BrückW.De GrootC. J. A.BergersE.HulshofS.GeurtsJ.et al (2003). Remyelinated lesions in multiple sclerosis: magnetic resonance image appearance. Arch. Neurol.60, 1073–1081. 10.1001/archneur.60.8.1073
8
BeenkenA.MohammadiM. (2009). The FGF family: biology, pathophysiology and therapy. Nat. Rev. Drug Discov.8, 235–253. 10.1038/nrd2792
9
BellanM.PirisiM.SainaghiP. P. (2016). The Gas6/TAM system and multiple sclerosis. Int. J. Mol. Sci.17 (11), 1807. 10.3390/ijms17111807
10
BüttnerM.MiaoZ.WolfF. A.TeichmannS. A.TheisF. J. (2019). A test metric for assessing single-cell RNA-seq batch correction. Nat. Methods16, 43–49. 10.1038/s41592-018-0254-1
11
CiccarelliO.BarkhofF.CalabreseM.De StefanoN.EshaghiA.FilippiM.et al (2024). Using the progression independent of relapse activity framework to unveil the pathobiological foundations of multiple sclerosis. Neurology103, e209444. 10.1212/wnl.0000000000209444
12
CignarellaF.FilipelloF.BollmanB.CantoniC.LoccaA.MikesellR.et al (2020). TREM2 activation on microglia promotes myelin debris clearance and remyelination in a model of multiple sclerosis. Acta Neuropathol.140 (4), 513–534. 10.1007/s00401-020-02193-z
13
ClementeD.OrtegaM. C.ArenzanaF. J.De CastroF. (2011). FGF-2 and Anosmin-1 are selectively expressed in different types of multiple sclerosis lesions. J. Neurosci.31, 14899–14909. 10.1523/JNEUROSCI.1158-11.2011
14
International Multiple Sclerosis Genetics Consortium (2019). Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility. Science365 (6460), eaav7188. 10.1126/science.aav7188
15
International Multiple Sclerosis Genetics Consortium; MultipleMS Consortium (2023). Locus for severity implicates CNS resilience in progression of multiple sclerosis. Nature619 (7969), 323–331. 10.1038/s41586-023-06250-x
16
Das NevesS. P.DelivanoglouN.RenY.CucuzzaC. S.MakuchM.AlmeidaF.et al (2024). Meningeal lymphatic function promotes oligodendrocyte survival and brain myelination. Immunity57, 2328–2343 e2328. 10.1016/j.immuni.2024.08.004
17
DavalosD.Kyu RyuJ.MerliniM.BaetenK. M.Le MoanN.PetersenM. A.et al (2012). Fibrinogen-induced perivascular microglial clustering is required for the development of axonal damage in neuroinflammation. Nat. Commun.3, 1227. 10.1038/ncomms2230
18
DoanR. A.MonkK. R. (2025). Dock1 functions in Schwann cells to regulate development, maintenance, and repair. J. Cell Biol.224, e202311041. 10.1083/jcb.202311041
19
FaureL.SoldatovR.KharchenkoP. V.AdameykoI. (2023). scFates: a scalable python package for advanced pseudotime and bifurcation analysis from single-cell data. Bioinformatics39, btac746. 10.1093/bioinformatics/btac746
20
FlemingS. J.ChaffinM. D.ArduiniA.AkkadA.-D.BanksE.MarioniJ. C.et al (2023). Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender. Nat. Methods20, 1323–1335. 10.1038/s41592-023-01943-7
21
FurushoM.RouloisA. J.FranklinR. J.BansalR. (2015). Fibroblast growth factor signaling in oligodendrocyte-lineage cells facilitates recovery of chronically demyelinated lesions but is redundant in acute lesions. Glia63, 1714–1728. 10.1002/glia.22838
22
GermainP. L.LunA.Garcia MeixideC.MacnairW.RobinsonM. D. (2021). Doublet identification in single-cell sequencing data using scDblFinder. F1000Res10, 979. 10.12688/f1000research.73600.2
23
GudiV.Moharregh-KhiabaniD.SkripuletzT.KoutsoudakiP. N.KotsiariA.SkuljecJ.et al (2009). Regional differences between grey and white matter in cuprizone induced demyelination. Brain Res.1283, 127–138. 10.1016/j.brainres.2009.06.005
24
GudiV.GingeleS.SkripuletzT.StangelM. (2014). Glial response during cuprizone-induced de- and remyelination in the CNS: lessons learned. Front. Cell Neurosci.8, 73. 10.3389/fncel.2014.00073
25
GudiV.GaiL.HerderV.TejedorL. S.KippM.AmorS.et al (2017). Synaptophysin is a reliable marker for axonal damage. J. Neuropathol. Exp. Neurol.76, 109–125. 10.1093/jnen/nlw114
26
HaoY.HaoS.Andersen-NissenE.MauckW. M.ZhengS.ButlerA.et al (2021). Integrated analysis of multimodal single-cell data. Cell184, 3573–3587.e3529. 10.1016/j.cell.2021.04.048
27
HarringtonE. P.ZhaoC.FancyS. P.KaingS.FranklinR. J.RowitchD. H. (2010). Oligodendrocyte PTEN is required for myelin and axonal integrity, not remyelination. Ann. Neurol.68, 703–716. 10.1002/ana.22090
28
HeL.Davila-VelderrainJ.SumidaT. S.HaflerD. A.KellisM.KulminskiA. M. (2021). NEBULA is a fast negative binomial mixed model for differential or co-expression analysis of large-scale multi-subject single-cell data. Commun. Biol.4, 629. 10.1038/s42003-021-02146-6
29
HendrickxD. a.E.Van ScheppingenJ.Van Der PoelM.BossersK.SchuurmanK. G.Van EdenC. G.et al (2017). Gene expression profiling of multiple sclerosis pathology identifies early patterns of demyelination surrounding chronic active lesions. Front. Immunol.8, 1810. 10.3389/fimmu.2017.01810
30
HorH.FrancescattoL.BartesaghiL.Ortega-CuberoS.KousiM.Lorenzo-BetancorO.et al (2015). Missense mutations in TENM4, a regulator of axon guidance and central myelination, cause essential tremor. Hum. Mol. Genet.24, 5677–5686. 10.1093/hmg/ddv281
31
HouJ.ZhouY.CaiZ.TerekhovaM.SwainA.AndheyP. S.et al (2023). Transcriptomic atlas and interaction networks of brain cells in mouse CNS demyelination and remyelination. Cell Rep.42, 112293. 10.1016/j.celrep.2023.112293
32
HrvatinS.TzengC. P.NagyM. A.StroudH.KoutsioumpaC.WilcoxO. F.et al (2019). A scalable platform for the development of cell-type-specific viral drivers. Elife8, e48089. 10.7554/eLife.48089
33
HuJ.LiX.ColemanK.SchroederA.MaN.IrwinD. J.et al (2021). SpaGCN: integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network. Nat. Methods18, 1342–1351. 10.1038/s41592-021-01255-8
34
ItohN.ItohY.TassoniA.RenE.KaitoM.OhnoA.et al (2017). Cell-specific and region-specific transcriptomics in the multiple sclerosis model: focus on astrocytes. Proc. Natl. Acad. Sci. U. S. A.115 (2), E302–E309. 10.1073/pnas.1716032115
35
IvanD. C.BerveK. C.WalthertS.MonacoG.BorstK.BouilletE.et al (2023). Insulin-like growth factor-1 receptor controls the function of CNS-resident macrophages and their contribution to neuroinflammation. Acta Neuropathol. Commun.11, 35. 10.1186/s40478-023-01535-8
36
JäkelS.AgirreE.FalcãoA. M.Van BruggenD.LeeK. W.KnueselI.et al (2019). Altered human oligodendrocyte heterogeneity in multiple sclerosis. Nature566 (7745), 543–547. 10.1038/s41586-019-0903-2
37
JanssenP.KliesmeteZ.ViethB.AdiconisX.SimmonsS.MarshallJ.et al (2023). The effect of background noise and its removal on the analysis of single-cell expression data. Genome Biol.24, 140. 10.1186/s13059-023-02978-x
38
JinS.Guerrero-JuarezC. F.ZhangL.ChangI.RamosR.KuanC. H.et al (2021). Inference and analysis of cell-cell communication using CellChat. Nat. Commun.12, 1088. 10.1038/s41467-021-21246-9
39
JollyS.BazarganiN.QuirogaA. C.PringleN. P.AttwellD.RichardsonW. D.et al (2018). G protein-coupled receptor 37-like 1 modulates astrocyte glutamate transporters and neuronal NMDA receptors and is neuroprotective in ischemia. Glia66, 47–61. 10.1002/glia.23198
40
KaufmannM.SchauppA. L.SunR.CosciaF.DendrouC. A.CortesA.et al (2022). Identification of early neurodegenerative pathways in progressive multiple sclerosis. Nat. Neurosci.25, 944–955. 10.1038/s41593-022-01097-3
41
KaunznerU. W.KangY.ZhangS.MorrisE.YaoY.PandyaS.et al (2019). Quantitative susceptibility mapping identifies inflammation in a subset of chronic multiple sclerosis lesions. Brain142, 133–145. 10.1093/brain/awy296
42
Keren-ShaulH.SpinradA.WeinerA.Matcovitch-NatanO.Dvir-SzternfeldR.UllandT. K.et al (2017). A unique microglia type associated with restricting development of Alzheimer’s disease. Cell169, 1276–1290 e1217. 10.1016/j.cell.2017.05.018
43
KippM. (2024). How to use the cuprizone model to study de- and remyelination. Int. J. Mol. Sci.25, 1445. 10.3390/ijms25031445
44
KleshchevnikovV.ShmatkoA.DannE.AivazidisA.KingH. W.LiT.et al (2022). Cell2location maps fine-grained cell types in spatial transcriptomics. Nat. Biotechnol.40 (5), 661–671. 10.1038/s41587-021-01139-4
45
KorsunskyI.MillardN.FanJ.SlowikowskiK.ZhangF.WeiK.et al (2019). Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods16, 1289–1296. 10.1038/s41592-019-0619-0
46
KotterM. R.LiW. W.ZhaoC.FranklinR. J. (2006). Myelin impairs CNS remyelination by inhibiting oligodendrocyte precursor cell differentiation. J. Neurosci.26, 328–332. 10.1523/JNEUROSCI.2615-05.2006
47
KuhlmannT.LudwinS.PratA.AntelJ.BruckW.LassmannH. (2017). An updated histological classification system for multiple sclerosis lesions. Acta Neuropathol.133, 13–24. 10.1007/s00401-016-1653-y
48
LangfelderP.LuoR.OldhamM. C.HorvathS. (2011). Is my network module preserved and reproducible?PLoS Comput. Biol.7, e1001057. 10.1371/journal.pcbi.1001057
49
LassmannH. (2018). Multiple sclerosis pathology. Cold Spring Harb. Perspect. Med.8, a028936. 10.1101/cshperspect.a028936
50
Lerma-MartinC.BadiaI. M. P.Ramirez FloresR. O.SekolP.SchaferP. S. L.RiedlC. J.et al (2024). Cell type mapping reveals tissue niches and interactions in subcortical multiple sclerosis lesions. Nat. Neurosci.27, 2354–2365. 10.1038/s41593-024-01796-z
51
LiuB.LiY.ZhangL. (2021). Analysis and visualization of spatial transcriptomic data. Front. Genet.12, 785290. 10.3389/fgene.2021.785290
52
LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15, 550. 10.1186/s13059-014-0550-8
53
MaoQ.YangL.WangL.GoodisonS.SunY. (2015). Simpleppt: A Simple Principal Tree Algorithm. Vancouver, Canada: Society for Industrial and Applied Mathematics Publications.
54
MarquesS.ZeiselA.CodeluppiS.Van BruggenD.Mendanha FalcãoA.XiaoL.et al (2016). Oligodendrocyte heterogeneity in the mouse juvenile and adult central nervous system. Science352, 1326–1329. 10.1126/science.aaf6463
55
MohanH.FrieseA.AlbrechtS.KrumbholzM.ElliottC. L.ArthurA.et al (2014). Transcript profiling of different types of multiple sclerosis lesions yields FGF1 as a promoter of remyelination. Acta Neuropathol. Commun.2, 168. 10.1186/s40478-014-0168-9
56
MorabitoS.MiyoshiE.MichaelN.ShahinS.MartiniA. C.HeadE.et al (2021). Single-nucleus chromatin accessibility and transcriptomic characterization of Alzheimer's disease. Nat. Genet.53, 1143–1155. 10.1038/s41588-021-00894-z
57
MorabitoS.ReeseF.RahimzadehN.MiyoshiE.SwarupV. (2023). hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data. Cell Rep. Methods3, 100498. 10.1016/j.crmeth.2023.100498
58
Orthmann-MurphyJ.CallC. L.Molina-CastroG. C.HsiehY. C.RasbandM. N.CalabresiP. A.et al (2020). Remyelination alters the pattern of myelin in the cerebral cortex. Elife9, e56621. 10.7554/eLife.56621
59
PataniR.HardinghamG. E.LiddelowS. A. (2023). Functional roles of reactive astrocytes in neuroinflammation and neurodegeneration. Nat. Rev. Neurol.19 (7), 395–409. 10.1038/s41582-023-00822-1
60
PereiraJ. B.JanelidzeS.StrandbergO.WhelanC. D.ZetterbergH.BlennowK.et al (2022). Microglial activation protects against accumulation of tau aggregates in nondemented individuals with underlying Alzheimer's disease pathology. Nat. Aging2, 1138–1144. 10.1038/s43587-022-00310-z
61
PonathG.ParkC.PittD. (2018). The role of astrocytes in multiple sclerosis. Front. Immunol.9, 217. 10.3389/fimmu.2018.00217
62
PraetJ.GuglielmettiC.BernemanZ.Van Der LindenA.PonsaertsP. (2014). Cellular and molecular neuropathology of the cuprizone mouse model: clinical relevance for multiple sclerosis. Neurosci. Biobehav Rev.47, 485–505. 10.1016/j.neubiorev.2014.10.004
63
ProescholdtM. A.JacobsonS.TresserN.OldfieldE. H.MerrillM. J. (2002). Vascular endothelial growth factor is expressed in multiple sclerosis plaques and can induce inflammatory lesions in experimental allergic encephalomyelitis rats. J. Neuropathol. Exp. Neurol.61, 914–925. 10.1093/jnen/61.10.914
64
RajendranR.RajendranV.BottigerG.StadelmannC.ShirvanchiK.Von AuL.et al (2023). The small molecule fibroblast growth factor receptor inhibitor infigratinib exerts anti-inflammatory effects and remyelination in a model of multiple sclerosis. Br. J. Pharmacol.180, 2989–3007. 10.1111/bph.16186
65
RansohoffR. M. (2012). Animal models of multiple sclerosis: the good, the bad and the bottom line. Nat. Neurosci.15, 1074–1077. 10.1038/nn.3168
66
RaoA.BarkleyD.FrancaG. S.YanaiI. (2021). Exploring tissue architecture using spatial transcriptomics. Nature596, 211–220. 10.1038/s41586-021-03634-9
67
ReichD. S.LucchinettiC. F.CalabresiP. A. (2018). Multiple sclerosis. N. Engl. J. Med.378, 169–180. 10.1056/NEJMra1401483
68
RothhammerV.BoruckiD. M.TjonE. C.TakenakaM. C.ChaoC. C.Ardura-FabregatA.et al (2018). Microglial control of astrocytes in response to microbial metabolites. Nature557, 724–728. 10.1038/s41586-018-0119-x
69
SchaffnerE.Bosch-QueraltM.EdgarJ. M.LehningM.StraussJ.FleischerN.et al (2023). Myelin insulation as a risk factor for axonal degeneration in autoimmune demyelinating disease. Nat. Neurosci.26, 1218–1228. 10.1038/s41593-023-01366-9
70
ScheyltjensI.Van HoveH.De VlaminckK.KanchevaD.BastosJ.Vara-PerezM.et al (2022). Single-cell RNA and protein profiling of immune cells from the mouse brain and its border tissues. Nat. Protoc.17, 2354–2388. 10.1038/s41596-022-00716-4
71
ShaoX.LiC.YangH.LuX.LiaoJ.QianJ.et al (2022). Knowledge-graph-based cell-cell communication inference for spatially resolved transcriptomic data with SpaTalk. Nat. Commun.13, 4429. 10.1038/s41467-022-32111-8
72
ShenK.ReicheltM.KyaukR. V.NguH.ShenY. A.ForemanO.et al (2021). Multiple sclerosis risk gene Mertk is required for microglial activation and subsequent remyelination. Cell Rep.34, 108835. 10.1016/j.celrep.2021.108835
73
SkripuletzT.GudiV.HackstetteD.StangelM. (2011). De- and remyelination in the CNS white and grey matter induced by cuprizone: the old, the new, and the unexpected. Histol. Histopathol.26, 1585–1597. 10.14670/HH-26.1585
74
StidworthyM. F.GenoudS.SuterU.ManteiN.FranklinR. J. (2003). Quantifying the early stages of remyelination following cuprizone-induced demyelination. Brain Pathol.13, 329–339. 10.1111/j.1750-3639.2003.tb00032.x
75
SunD.GuanX.MoranA. E.WuL. Y.QianD. Z.SchedinP.et al (2022). Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat. Biotechnol.40, 527–538. 10.1038/s41587-021-01091-3
76
ThurnherrT.BenningerY.WuX.ChrostekA.KrauseS. M.NaveK. A.et al (2006). Cdc42 and Rac1 signaling are both required for and act synergistically in the correct formation of myelin sheaths in the CNS. J. Neurosci.26, 10110–10119. 10.1523/JNEUROSCI.2158-06.2006
77
ThurstonG.DalyC. (2012). The complex role of angiopoietin-2 in the angiopoietin-tie signaling pathway. Cold Spring Harb. Perspect. Med.2, a006550. 10.1101/cshperspect.a006650
78
Vega-RiquerJ. M.Mendez-VictorianoG.Morales-LuckieR. A.Gonzalez-PerezO. (2019). Five decades of cuprizone, an updated model to replicate demyelinating diseases. Curr. Neuropharmacol.17, 129–141. 10.2174/1570159X15666170717120343
79
WebsterH. D. F. (1997). Growth factors and myelin regeneration in multiple sclerosis. Mult. Scler.3, 113–120. 10.1177/135245859700300210
80
WeingerJ. G.OmariK. M.MarsdenK.RaineC. S.Shafit-ZagardoB. (2009). Up-regulation of soluble Axl and Mer receptor tyrosine kinases negatively correlates with Gas6 in established multiple sclerosis lesions. Am. J. Pathol.175, 283–293. 10.2353/ajpath.2009.080807
81
WelchJ. D.KozarevaV.FerreiraA.VanderburgC.MartinC.MacoskoE. Z. (2019). Single-Cell multi-omic integration compares and contrasts features of brain cell identity. Cell177, 1873–1887.e1817. 10.1016/j.cell.2019.05.006
82
WheelerM. A.ClarkI. C.TjonE. C.LiZ.ZandeeS. E. J.CouturierC. P.et al (2020). MAFG-driven astrocytes promote CNS inflammation. Nature578 (7796), 593–599. 10.1038/s41586-020-1999-0
83
WheelerM. A.ClarkI. C.LeeH. G.LiZ.LinnerbauerM.RoneJ. M.et al (2023). Droplet-based forward genetic screening of astrocyte-microglia cross-talk. Science379, 1023–1030. 10.1126/science.abq4822
84
WolfF. A.AngererP.TheisF. J. (2018). SCANPY: large-scale single-cell gene expression data analysis. Genome Biol.19, 15. 10.1186/s13059-017-1382-0
85
XieM. (2010). Rostrocaudal analysis of corpus callosum demyelination and axon damage across disease stages refines diffusion tensor imaging correlations with pathological features. J. Neuropathol. Exp. Neurol.69 (7), 704–716. 10.1097/NEN.0b013e3181e3de90
86
ZeiraR.LandM.StrzalkowskiA.RaphaelB. J. (2022). Alignment and integration of spatial transcriptomics data. Nat. Methods19, 567–575. 10.1038/s41592-022-01459-6
87
ZhangB.HorvathS. (2005). A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol.4, Article17. 10.2202/1544-6115.1128
88
ZhangY.ParmigianiG.JohnsonW. E. (2020). ComBat-seq: batch effect adjustment for RNA-seq count data. Nar. Genom Bioinform2, lqaa078. 10.1093/nargab/lqaa078
89
ZhaoX.EyoU. B.MuruganM.WuL.-J. (2018). Microglial interactions with the neurovascular system in physiology and pathology. Dev. Neurobiol.78, 604–617. 10.1002/dneu.22576
90
ZhuJ.SunY. H.OuyangZ.LiK.NegiS.PiyaS.et al (2023). RNASequest: an end-to-end reproducible RNAseq data analysis and publishing framework. J. Mol. Biol.435, 168017. 10.1016/j.jmb.2023.168017
Summary
Keywords
cell-cell interaction, cuprizone model, demyelination, multiple sclerosis, oligodendrocytes, remyelination, snRNA-seq, spatial transcriptomics (ST)
Citation
Tsai H-H, Piya S, Wang J, Zhu J, Hu W, Gehrke AR, Cao S, Guise AJ, Chan SJ, Sheehan M, Chu J, Ouyang Z, Ryals M, Lee M, Wang W, Zhao E, Cullen P, Challa R, Marshall E, Zeng W, Kaeser-Woo YJ, Ehrenfels C, Jandreski L, McLaughlin H, Carlile TM, Gagnon J, Reynolds TL, Li M, Li K and Zhang B (2026) Spatial transcriptomics reveal heterogeneous cell‒cell interactions among brain regions in cuprizone model consistent with multiple sclerosis lesions. Front. Bioinform. 6:1832826. doi: 10.3389/fbinf.2026.1832826
Received
17 March 2026
Revised
08 June 2026
Accepted
29 June 2026
Published
06 August 2026
Volume
6 - 2026
Edited by
Zhaokai Zhou, Central South University, China
Reviewed by
Zidong Zhang, Icahn School of Medicine at Mount Sinai, United States
Sachin Gadani, University of Pittsburgh, United States
Updates
Copyright
© 2026 Tsai, Piya, Wang, Zhu, Hu, Gehrke, Cao, Guise, Chan, Sheehan, Chu, Ouyang, Ryals, Lee, Wang, Zhao, Cullen, Challa, Marshall, Zeng, Kaeser-Woo, Ehrenfels, Jandreski, McLaughlin, Carlile, Gagnon, Reynolds, Li, Li and Zhang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hui-Hsin Tsai, huihsin.tsai@biogen.com; Baohong Zhang, baohongz@gmail.com; Kejie Li, likejie@gmail.com
† These authors share first authorship
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.