METHOD DETAILS

Isolation of nuclei was performed similarly as previously described.³⁰ Briefly, frozen brain cortex and white matter were Dounce homogenized in 5 mL of lysis buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl2, and 0.025% NP-40) for 15 min. Following initial Dounce homogenization, the solution was then filtered through a 30-μm cell filter and pelleted at 500g for 5 min at 4°C. Nuclei were washed and filtered twice with nuclei wash buffer (1% BSA in PBS with 0.2 U μl⁻¹ RNasin (Promega)). Nuclei pellets were resuspended in 500 μL nuclei wash and 900 μL 1.8 M sucrose. This 1,400 μL mixture was carefully layered on top of 500 μL 1.8 M sucrose and centrifuged at 13,000g for 45 min at 4°C to separate the nuclei from myelin debris. The nuclei pellet was resuspended in nuclei wash buffer at 1,200 nuclei μl⁻¹ and filtered through a 40-μm FlowMi Cell Strainer.

Droplet-based snRNA-seq was performed using the Chromium Single Cell 5^’ Reagent Kits (v1 for WT C57BL/6J snRNA-seq dataset and v2 for Trem2^−/− snRNA-seq dataset) per the manufacturer’s instructions (10x Genomics). Nuclei were resuspended to a concentration of 1,200 nuclei per μl before loading according to the manufacturer’s protocol. The libraries were sequenced using Illumina sequencers (NovaSeq instrument) with 150bp paired-end sequencing at the McDonnell Genome Institute. Sample demultiplexing, barcode processing and single-cell counting were performed using the Cell Ranger Single-Cell Software Suite (10x Genomics). To include Clec7a gene, we built a custom reference genome by adding Clec7a cDNA to the pre-built mouse reference genome (mm10). Cellranger count (v.3.0.2) was used to align reads to the custom reference genome, quantify reads and filter reads with a quality score below 30.

The Seurat package in R was used for subsequent analysis.¹⁰⁵ For quality control, nuclei with mitochondrial content >5% were removed. Nuclei that are doublets or multiplets were filtered out by two steps. First, nuclei with high UMI and gene number per nucleus were filtered out. Cutoffs for UMI and gene number were determined on the basis of histograms showing cell density as a function of UMI per gene counts. Then nuclei with more than one cell type marker gene expressed were removed. For WT C57BL/6J mouse brain snRNA-seq analysis, a cutoff of 500–10,000 UMI and 500–4,000 genes was applied. After filtering, a total of 58,079 individual nuclei across all conditions remained, with a median of 1,808 UMIs and 1,250 genes per nucleus for downstream analysis. For Trem2^−/− and littermate WT mouse brain snRNA-seq analysis, a cutoff of 500–15,000 UMI and 400–5,000 genes was applied. After filtering, a total of 117,729 individual nuclei across all genotypes remained, with a median of 3,082 UMIs and 1,951 genes per nucleus for downstream analysis. Data were log normalized and regressed on mitochondrial gene percentage during data scaling, using the mitochondria ratio as the argument for ‘vars.to.regress’. Samples were batch corrected using FindIntegrationAnchors function and Canonical Correlation Analysis (CCA). Principal component analysis was performed using the top 3,000 most variable genes and UMAP analysis was performed with the top 20 PCAs. Clustering was performed using a resolution of 0.3. For identifying markers for each cluster, we performed differential expression of each cluster against all other clusters, identifying positive markers for that cluster, using FindAllMarkers() function in the Seurat package with default parameters. Cluster marker genes are listed in Tables S1 and S4. Nuclei from broad cell types (neurons, oligodendrocytes, astrocytes, microglia and vascular cells) were taken and re-clustered to further analyze the sub-clusters in each cell type. For data visualization, BBrowser (version 2) was also used.

To identify subsets in each cell type, nuclei from all samples belonging to a given major cell type were extracted for downstream analysis. Doublets were filtered out as described above. Re-integration was applied to microglia sub-clustering (both C57BL/6J and Trem2^−/− datasets), and oligodendrocyte sub-clustering (Trem2^−/− dataset). Principal component analysis was performed prior to clustering and the first ten principal components were used based on the ElbowPlot. Clustering was performed using the FindClusters() function. For C57BL/6J dataset, clustering resolution at 0.3 was used for oligodendrocyte lineage, 0.6 for astrocytes, and 0.6 for microglia. For vascular cells, the first 15 principal components were used, and clustering was done with resolution 1.1. For Trem2^−/− dataset, clustering resolution at 0.3 was used for oligodendrocyte lineage, and 0.6 for microglia.

Differential expression of genes between conditions was done using FindMarkers() function with default parameters of the Seurat package in R,¹⁰⁶ which is based on the non-parametric Wilcoxon rank-sum test. Log2(fold change) of average expression and the percentage of cells (pct) expressing the genes in each condition were generated. The adjusted p value was calculated using Bonferroni correction. Genes with log2(Fold Change) > 0.5, adjusted p < 0.05 were considered as significant DEGs. Lists of DEGs were generated by filtering all genes for log2 (fold change) > 0.25 (Tables S2 and S5). These gene lists were used as inputs for downstream pathway analyses.

Average expression per cluster per sample was calculated using the AverageExpression() function in the Seurat package. For heat maps of relative gene expression across cell types, average expression of each gene was Z score transformed across samples and plotted using heat.map2 function in the gplots packge (v3.1.1) in R.

Ligand-receptor interactions were mined from the snRNA-seq data using the NicheNet algorithm (nichenetr 1.0.0).⁷⁶ Briefly, sender and receiver cells were defined and indicated in Figure 6. NicheNet analysis was performed using upregulated DEGs in the defined receiver cells between demyelination versus normal (log2(fold change) > 0.5, adjusted p value <0.05) as input with default settings. Circos plot was computed based on calculation of interaction scores between possible ligand-receptor combinations. Interaction pairs were filtered for those bona fide interactions that were documented in the literature and publicly available datasets.

Single-cell pseudotime trajectories were inferred using the R package slingshot (v2.1.0).¹⁰⁰ All C57BL/6J microglia nuclei except for the proliferating microglia were used for the analysis.

RNA velocity analysis was performed using the package scVelo (v.0.2.4) with stochastic modeling, as instructed.³⁸ Briefly, we used the.bam and.bam.bai files from the CellRanger output and generated the.loom files for the representative samples in each condition (demyelination: A3, and remyelination: B7) using the velocyto.py version 0.17 with the velocyto run10x command. Next, the.loom files were loaded into Python (v3.9.13) using scv.read function to generate count tables in JupyterLab interface. The Seurat objects from each genotype were converted into anndata files as described (https://smorabit.github.io/tutorials/8_velocyto/) and then loaded separately into JupterLab interface with Python via the function sc.read.h5ad. RNA velocity was computed via scv.tl.velocity(), and the final velocity streams and pseudotime were generated using the commands scv.pl.velocity_embedding_stream() and scv.pl.velocity_pseudotime(), respectively. For visualization, the RNA velocity stream was projected onto each predefined UMAPs.

To compare oligodendrocyte gene signatures in the cuprizone model to that observed in mouse juvenile CNS²⁶ as well as DOL signature characterized in the 5XFAD model,³⁹ we re-analyzed the raw³⁹ or filtered²⁶ matrix data of the two aforementioned datasets. Filtering parameters were the same as in the original manuscript. After filtering, each dataset was processed using functions implemented in the Seurat package (SCTransform, RunPCA, RunUMAP, FindNeighbors, FindClusters). Oligodendrocyte lineage cells were extracted using subset() function in the Seurat package. Samples were integrated using FindIntegrationAnchors() function and CCA.¹⁰⁷ CPZ oligodendrocyte dataset was integrated with Marques et al.²⁶ dataset and Kenigsbuch et al.³⁹ dataset, respectively. Principal component analysis was performed using the top 3,000 most variable genes and UMAP analysis was performed with the top 20 PCAs. Annotation of subclusters were based on original reference or clusters characterized in the current manuscript. Dot plot of marker gene expression was generated using scanpy (v1.9.3) and dendrogram was built by performing hierarchical clustering on predefined subclusters.

To reanalyze publicly available datasets on human patients with MS, raw or filtered matrix of UMI counts and associated metadata were downloaded from GEO under accession number {"type":"entrez-geo","attrs":{"text":"GSE124335","term_id":"124335"}}GSE124335, {"type":"entrez-geo","attrs":{"text":"GSE118257","term_id":"118257"}}GSE118257 and {"type":"entrez-geo","attrs":{"text":"GSE180759","term_id":"180759"}}GSE180759, and analyzed using Seurat package. For oligodendrocyte signature comparisons, oligodendrocyte lineage cells were extracted from the two datasets, respectively.^6,8 Differential expression analysis was performed between each pathological state and the control as described above. DEGs were considered as genes with log2 (fold change) > 0.25, adj. p value <0.05 for {"type":"entrez-geo","attrs":{"text":"GSE124335","term_id":"124335"}}GSE124335 dataset, with log2 (fold change) > 0.25, adj. p value <0.05 for {"type":"entrez-geo","attrs":{"text":"GSE118257","term_id":"118257"}}GSE118257 dataset, and with log2 (fold change) > 0.5, adj. p value <0.05 for {"type":"entrez-geo","attrs":{"text":"GSE180759","term_id":"180759"}}GSE180759 dataset. For astrocyte signature comparisons, top 100 marker genes of the AIMS were extracted from Absinta et al., 2021.⁸ For microglia comparisons, top 94 marker genes of the MIMS-foamy isolated from the chronic lesions of patients with MS,⁸ or top 100 feature genes of microglia from acute lesions of patients with MS²² were extracted. Gene sets listed in Table S3 were used as inputs for this analysis. Human genes were converted to mouse genes using biomaRt package (2.48.3) in R. Venn diagrams were plotted using BioVenn package (v1.1.3) in R.

Gene set scores were calculated using the AddModuleScore() function in Seurat package. Briefly, for each cell, the log-transformed average expression of all genes in each gene set was calculated and subtracted by the log-transformed expression of control features sets. The white matter-associated microglia (WAM) score was calculated by using the top 100 signature genes of WAM identified from aged mice.⁵⁶ The disease-associated microglia (DAM) score was determined by using the top 71 signature genes of DAM identified in the 5XFAD model.⁵⁵ The reactive oligodendrocyte signature was defined by the cluster markers of DOLs cluster. Gene signatures of myelination was obtained from the “Myelination” pathway from Metascape pathway analysis. Gene sets listed in Table S3 were used as inputs for this analysis.

For pathway enrichment analysis, enrichment for Gene Ontology terms was obtained through Metascape.³⁰ Gene set enrichment analysis (GSEA) was performed using the hallmark pathways from the Molecular Signature Database.

Animals were anesthetized by 10 mg/mL ketamine and 1 mg/mL xylazine solution i.p. and perfused transcardially with PBS containing 1U/ml heparin. Fixation of brain samples was done in 4% paraformaldehyde (PFA) overnight. The brain tissues were then incubated in 30% sucrose in PBS for 24 h. After freezing the tissue on dry ice using Tissue-Tek O.C.T and 30% sucrose at a ratio of 1:2, 20-μm coronal sections were cut by a Leica CM 1900 cryostat. Floating mouse brain sections were blocked with 3% BSA and 0.25% Triton X-100 in PBS for 1h at room temperature, and stained with anti-IBA1 (1:500; rabbit, #17198S, Cell Signaling Technology), anti-IBA1 (1:500; goat, #5076, Abcam), anti-CD74 (1:200; Alexa Fluor 647-labeled rat IgG2b, #151004, Biolegend), anti-OLIG2 (1:500; Rabbit, #AB9610, EMB Millipore), anti-CA2 (1:200, rat, #MAB2184, R&D), anti-CC1 (1:50, Mouse, #OP80, Calbiochem), anti-SER-PINA3N (1:300; Goat, #AF4709, R&D), anti-IL33 (5 μg/mL; Goat, #AF3626, R&D), anti-GFAP (1:1000; Alexa Fluor 488-labbeled mouse IgG1, #53–9892-82, ThermoFisher Scientific), anti-APOE (1:300; biotinylated, mouse monoclonal: HJ6.3, a gift from Dr. David M. Holtzman), anti-CD11c (1:500; Rabbit, #97585, Cell Signaling Technology), anti-MAFB (1:400; Rabbit, #A700–046, BETHYL), anti-Vimentin (1:500; Mouse, #sc-373717, Santa Cruz Biotechnology), anti-STAT1 (1:300; Rabbit, #14994, Cell Signaling Technology), anti-HMGCS1 (1:200; Rabbit, #PA5–29488, Invitrogen) and anti-Synaptophysin (1:200; Rabbit, #36406, Cell Signaling Technology) overnight at 4°C, followed by staining with anti-rabbit IgG Alexa Fluor 555 (1:2,000; #ab150074, Abcam), anti-rat IgG Alexa Fluor 647 (1:2,000; #A21472, Invitrogen), anti-mouse IgG2b Alexa Fluor 488 (1:2,000; #A21141, Invitrogen), anti-goat IgG Alexa Fluor 488 (1:2,000; #ab150129, Abcam), anti-Streptavidin, conjugate (1:2,000; #s868, ThermoFisher Scientific), and DAPI (1 μg/ml; #D9542, Sigma) for 1 h at room temperature. All antibodies were used in blocking buffer, and between all incubations, sections were washed for 10 min in PBS 3 times. Images were collected using a Nikon A1Rsi confocal microscope or Zeiss LSM880 microscope. z-Stacks with 1-μm steps in the z direction, 1,024 × 1,024-pixel resolution, were recorded. Extraction of parameters were performed in Imaris (version 9.8.0), and further processing was performed using automated scripts in MATLAB.⁸⁸

To quantify staining signals in the confocal images by volume in MCC or LCC, we created surfaces using the “skip automatic function, edit manually” option in Imaris and manually selected the MCC or LCC area on the top and bottom sections; we then used the “masked channel” function to create new channels that selectively contain signals of protein of interest inside the selected regions. Then “batch colocalization” of MATLAB function was run on the newly generated channels and colocalization channels were created. The volumes of signal of interest in each image were determined using “Surface” function in Imaris. Ratios of double-positive volumes over single positive volumes (IBA1⁺CD74+/IBA1⁺, IBA1⁺CD11C/IBA1⁺, IBA1⁺APOE⁺/IBA1⁺ and SERPINA3N⁺GFAP⁺/GFAP⁺) were then calculated. To analyze SERPINA3N⁺ volume in Trem2^−/− and WT mice, total SERPINA3N⁺ volume was measured within 2-μm radius of OLIG2, then ratios of total SERPINA3N⁺ volume over OLIG2⁺ number were calculated. To quantify VIM⁺ and Synaptophysin⁺ area in the LCC and cortex, images were first converted to a 2D image and then were converted to a binary image with the ImageJ “Make Binary” tool. The area of LCC or cortex were selected manually on the binary image using the ImageJ “Freehand Selections” tool. The mean signal intensity value was measured by the ImageJ “Measure” tool, using “mean gray value” as the readout.

Brains from 4 mice in each group were dissected and homogenized by a Dounce homogenizer in the dissociation buffer (1%BSA, 1mM EDTA). Myelin debris was removed by 30% Percoll density gradient centrifugation, followed by the labeling of microglia with CD45-APC-Cy7 (1:200, #103116, BioLegend) and IL-33R (ST2)-PE (1:200; #12-9333-80, Invitrogen). Dead cells were excluded using DAPI (1 μg/ml; D9542, Sigma). Cells were acquired using BD FACSCantoII flow cytometer. Microglia were gated on CD45^low cells and analyzed using Flowjo (v10.8.1).

Brains were processed for transmission electron microscopy (TEM) analysis as previously described.¹⁰⁸ Briefly, animals were perfused with the 4% PFA (Electron Microscopy Sciences, #15714-S), and post-fixed in the fixative containing 2.5% glutaraldehyde/2% paraformaldehyde/0.15M cacodylate buffer with 2mM CaCl2 (Electron Microscopy Sciences, #16300, #15710) at 4°C for 24 h before further dissection. Post fixation, brain samples were cut in sagittal direction near the midline into 100 μm thick sections with a vibratome (Leica VT1200S, Vienna, Austria). Whole brain sections containing CC were rinsed in 0.15 M cacodylate buffer containing 2 mM calcium chloride 3 times for 10 min each followed by a secondary fixation in 1% osmium tetroxide/1.5% potassium ferrocyanide in 0.15 M cacodylate buffer containing 2 mM calcium chloride for 1 h in the dark. The samples were then rinsed 3 times in ultrapure water for 10 min each and en bloc stained with 2% aqueous uranyl acetate overnight at 4°C in the dark. After another 4 washes in ultrapure water, the samples were dehydrated in a graded ethanol series (30%, 50%, 70%, 90%, 100% x3) for 10 min each step. Once dehydrated, samples were infiltrated with LX112 resin (Electron Microscopy Sciences) and flat embedded and polymerized at 60°C for 48 h. Post curing, medial regions of the CC were excised and mounted on blank epoxy stubs for cross sectioning. 70 nm sections were then cut, post-stained with 2% aqueous uranyl acetate and Sato’s lead and imaged on a TEM (JEOL JEM-1400 Plus) at an operating voltage of 120 Kv.¹⁰⁹

For each TEM analysis, five to ten randomly selected CC regions in each mouse were imaged. Data analysis was performed using ImageJ 1.53a software. The myelinated ratios and g-ratios were measured with ImageJ on transverse electron micrographs at 3,000× magnification. The perimeters of each axon and the myelin sheath were measured with the freehand tool by tracing the outer surfaces of each structure. The g-ratio was calculated as the perimeter of the axon over the perimeter of its myelin.

Myelin sheaths were stained with luxol fast blue (Electron microscopy sciences, Luxol Fast Blue, 0.1% in 95% alcohol, #26056–15). For the measurement of LFB intensity, the micrograph was first converted to a binary image with the ImageJ “Make Binary” tool. The area of MCC was then selected manually on the binary image using the ImageJ “Freehand Selections” tool. The mean LFB intensity value was measured by the ImageJ “Measure” tool, using “mean gray value” as the readout. For the quantification of LFB intensity, the mean intensity values in 3 adjacent brain sections from one animal were calculated and used for further statistical analysis. At least 3 animals per genotype were analyzed and the n was specified in the figure legends.

Microglia detached from the mixed glia culture were seeded in 24-well plate at a density of 100,000 cells/well in complete RPMI 1640 medium containing 10% LCCM. After one day, microglia were stimulated with myelin debris (500 μg/ml) for 48hrs, and the microglial lysis was harvested for PROS1 ELISA (#AE25758MO, Abebio). Total protein concentration was quantified by DC protein assay (#5000112, BIO-RAD).

For fluorescent labeling of thymocyte membranes, WT thymocytes were incubated with CFSE dye (5μM; C1157, ThermoFisher Scientific) for 8 min with constant agitation at room temperature. CFSE-labeled thymocytes were then incubated with 2μM Dexamethasone (D1756, Sigma) for 6 h in 37°C incubator to induce apoptosis. Our preliminary experiment confirmed that around 70% CFSE-labeled thymocytes became apoptotic after dexamethasone incubation by Annexin V and 7AAD staining (#559763, BD Biosciences). The CFSE-labeled apoptotic thymocytes were then incubated for 10 min with or without 10nM recombinant mouse PROS1 (#9740-PS-050, R&D systems), added to the astrocyte culture at a ratio of 20:1 (apoptotic thymocytes: astrocytes), and incubated for 12 h at 37°C. Astrocytes were then briefly washed in PBS for 3 times, incubated for 3 min at 37°C in 0.25% trypsin and detached by vigorous pipetting. Phagocytic astrocytes were assessed using LSR Fortessa (BD Biosciences). Dead cells were excluded by PI staining. The post-acquisition data was analyzed by FlowJo software (v10.8.1).