Spinal Cord Injury Induces Permanent Reprogramming of Microglia into a Disease-Associated State Which Contributes to Functional Recovery

Microglia are resident myeloid cells of the CNS. Recently, single-cell RNA sequencing (scRNAseq) has enabled description of a disease-associated microglia (DAM) with a role in neurodegeneration and demyelination. In this study, we use scRNAseq to investigate the temporal dynamics of immune cells harvested from the epicenter of traumatic spinal cord injury (SCI) induced in female mice. We find that as a consequence of SCI, baseline microglia undergo permanent transcriptional reprogramming into a previously uncharacterized subtype of microglia with striking similarities to previously reported DAM as well as a distinct microglial state found during development. Using a microglia depletion model we showed that DAM in SCI are derived from baseline microglia and strongly enhance recovery of hindlimb locomotor function following injury. SIGNIFICANCE STATEMENT Although disease-associated microglia (DAM) have been the subject of strong research interest during recent years (Keren-Shaul, 2017; Jordão, 2019), their cellular origin and their role in “normal” acute injury processes is not well understood. Our work directly addresses the origin and the role of DAM in traumatic injury response. Further, we use a microglia depletion model to prove that DAM in spinal cord injury (SCI) are indeed derived from homeostatic microglia, and that they strongly enhance recovery. Thus, in this work we significantly expand the knowledge of immune response to traumatic injury, demonstrate the applicability to human injury via our unique access to injured human spinal cord tissue, and provide the community with a comprehensive dataset for further exploration.

Although DAM have been the subject of strong research interest during recent years, their cellular origin and their role in acute injury processes is not well understood. In this study, we investigated the response of healthy tissue to injury, in the context of traumatic spinal cord injury (SCI). To obtain a detailed view of the temporal dynamics of immune cells after SCI, we applied scRNAseq to analyze 3687 single cells taken from healthy tissue and from injured spinal cord at a wide range of time points ranging from 0.5 h to 90 d post-SCI. Thus, we were able to trace an unbroken path between the acute and chronic tissue state of the injury. We determined that following SCI baseline microglia undergo permanent transcriptional reprogramming, resulting in a subtype of microglia with striking similarities to DAM, cells which we could also observe in human SCI. Using a microglial depletion model we showed that DAM in SCI are derived from baseline microglia and do in fact enhance recovery of hindlimb locomotor function in mice subjected to SCI.

Experimental design
Main experiment C57BL/6J mice were subjected to severe contusion SCI. At 0.5 h, 2 h, 6 h, 24 h, 36 h, 3 d, 7 d, 21 d, and 90 d post-SCI the injury epicenters were harvested. Healthy tissue was included as reference. Following tissue dissociation all CD45 1 cells were isolated and sequenced individually. Three to five animals were used per time point. Cell suspensions from the individual animals were pooled into one sample for each time point.SCI epicenter tissue sections were stained for APOE, CST7, IGF1, IBA1, and TMEM119 at each time point. Healthy tissue was included as reference. Three to five animals were investigated for each time point and stain.Proliferation of CD45 1 cells was investigated at 7 and 21 d post-SCI cells using flow cytometry. Three animals were investigated at each time point.
Origin of DAM Cx3cr1 1 cells were depleted by administration of tamoxifen to C57BL/ 6J-CX3CR1 CreER/1 Rosa26 DTA/1 mice at 9, 8, and 7 d prior to SCI. Severe contusion SCI was induced and CD45 1 cells were isolated and sequenced individually at 7 and 21 d post-SCI. Hind limb function was assessed during 90 d post-SCI. Wild-type animals were used as control. Each experimental group contained three or four animals.

Human tissue
Sections of injured human spinal cord were stained for APOE, CST7, IGF1, IBA1.
Bulk sequencing of spinal cord CX3CR11CD19-LY6G-CD3-NK1.1-SIGLEC-F-cells were isolated and subjected to sequencing in bulk at 7 and 21 d post-SCI. Healthy tissue was used as reference.

Mice
Wild-type female C57BL/6J mice (Scanbur), CX3CR1-GFP (005582, The Jackson Laboratory) and CX3CR1 CreER/1 Rosa26 DTA/1 (in-house breeding of 021160 and 009669 from The Jackson Laboratory by the R. A.H. lab) were used. Animal care and experiments were conducted in line with ethical permits (N196/15, N12317-2017, N138/14) approved by the local ethical committee (Stockholm, Sweden) and in line with local guidelines at Karolinska Institutet. Adult (10-15 weeks) females were used in all experiments.

Genotyping
Cx3cr1-gfp Protocols from The Jackson Laboratory were implemented. Briefly, ear biopsies were incubated for 1 h at 95°C in 100 ml solution [1:200 5 M NaOH (in-house), 1:2500 0.5 M EDTA (Invitrogen, 15575-038) in Milli-Q]. Samples were placed on ice and 100-ml solution (1:25 Tris-HCl; Merck Millipore, T5941) was added. 4 ml of the suspension was added to 16 ml of master mix (62.5% REDExtract-N-Amp PCR (Sigma, R4775), 28% nuclease free-water, 3% primer mutant, 3% primer wild type and 3% primer common). Working stock of primers (Eurofins Genomics; Table  1) was added at 20 mM. Suspension was added to a PCR plate (Microamp Optical 96-well Reaction Plate, Life Technologies N8010560). After centrifugation at 200 Â g for 60 s, the PCR plate was loaded onto a PCR instrument. A thermal cycling protocol was initiated with polymerase activation and denaturation at 94°C for 2 min, followed by 10 cycles of denaturation at 94°C for 20 s and annealing/extension at 65°C for 15 s (decreased 0.5°C/cycle) and 68°C for 10 s. This was followed by 28 cycles of denaturation at 94°C for 15 s, annealing/extension at 60°C for 15 s and 72°C for 10 s followed by 2 min at 72°C. Following thermal cycling the samples were added into an agarose gel (E-Gel Agarose Gel with SYBR Safe DNA Gel stain 2%, Invitrogen G521802) in combination with 20 ml DNA ladder (GeneRuler 100 bp, Invitrogen SM0241) and connected to a power-supply (Bio-Rad, Power-Pac 100) by which 90 V was delivered for 15-20 min. Heterozygotes were identified as expressing fluorescence at 337 and 485 bp. Mutants (485 bp only) and wild type (337 bp only) were killed. Cx3cr1 creer/1 Rosa26 DTA/1 DNA was prepared from ear biopsies. Using specific primers (Table 1) and PCR amplification, the genotype for each mouse was determined. The PCR reaction was conducted using a thermal cycler (Applied Biosystems, 2720) at 94°C for 5 min; 35 cycles at 94°C for 30 s, 60°C for 60 s, and 72°C for 60 s followed by 72°C for 2 min. The DNA fragments were measured on an agarose gel (4%) using a GelRed fluorescent nucleic acid gel stain.

Induction of SCI in mice
Intraperitoneal injection anesthesia consisting of 0.5 mg/kg medetomidin intraperitoneal (Domitor vet., Orion Pharma Animal Health 1 mg/ ml) and 75 mg/kg ketamine intraperitoneal (Ketador vet., Salfarm Scandinavia, 100 mg/ml) was administrated. Analgesics were administrated subcutaneously at induction using 0.05 mg/kg buprenorfin subcutaneously (temgesic, 0.3 mg/ml) and 5 mg/kg karprofen subcutaneously (Rimadyl, 50 mg/ml). The skin above the thoracic vertebrae was incised, the paravertebral muscles separated and the vertebral column stabilized using bilateral fixators in a stereotaxic frame (Model 900 & 900-c, Kopf). The dorsal part of the vertebrae was removed using a surgical drill (Anspach, EMAX 2). A 75 kdyn contusion SCI was induced at Th10 (Infinite Horizon, IH-0400). Bupivakain (Marcain, Aspen Nordic 2.5 mg/ml) was injected into the muscles surrounding the injury and the  moIMR0003  GGG CCA GCT CAT TCC TCC CAC TCA T  oIMR4281  CCA CAG GAT TTC AGC CTG AAC TTT G  oIMR4282  CGT GCA CTA TGC TCA GAT ATC TGT  skin was sutured (Ethicon, Vicryl 4.0). Post-operative analgesic treatment and care was given twice daily for a total of 3 d after surgery. Bladders were compressed manually until reflexive bladder emptying returned. The humane endpoint was set at 25% weight loss and was evaluated weekly.
Harvesting mouse spinal cords Animals were euthanized using a lethal dose of pentobarbitalnatrium (Apotek Produktion & Laboratorier, 60 mg/ml) and transcardially perfused with ice-cold 1Â DPBS using a peristaltic pump (Watson Marlow, 120S). Animals for histologic evaluation were further perfused with 4% paraformaldehyde (PFA; Thermo Fisher Scientific 28908) and postfixed in 4% PFA at 4°C O/N. Tissue for single-cell analysis was harvested from the SCI epicenter, defined as the site of impact 63 mm along the cranial-caudal axis in the mouse. Even in our highly controlled setting, the composition of the injury (size of cystic cavity, size of scar, degree of hematoma, force transmission and injury severity) varied slightly between the animals. Thus, a pool of three to five mice were used for each time point. This strategy allowed us to capture the heterogeneity in the immune cell response more accurately at the different time points and at the same time being able to pursue a cost-efficient and time-efficient approach.
Harvesting human spinal cords Human SCI samples were retrieved from an ongoing clinical trial [Safety and Efficacy of SC0806 (fibroblast growth factor 1 and a device) in traumatic SCI subjects] approved by the Swedish Medical Products Agency Study Protocol SC0806-A101. The study is currently ongoing at Karolinska University Hospital (Stockholm, Sweden) using ethical permits 2013/2257-31/1, 2015/1436-31, and 2019-03671. All research subjects gave their signed written informed consent. Briefly, the injured segment of the spinal cord was surgically removed in patients with SCI in the chronic stage. The rostral 1/3 of the dissected tissue was used for histologic evaluation of immune cells.

Dissociation of spinal cords
The SCI epicenters from three to five spinal cords/mice from the same treatment groups were pooled. This approach was implemented to limit variability within each time point. Spinal cords were dissociated using enzymatic dissociation with Papain (10 U/ml; Worthington, L5003126) in L-15 medium (Thermo Fisher Scientific, 11415064). Single-stranded and double-stranded DNA were digested using 200 U/ml DNase (Sigma, D7291). Papain was inactivated using 1% BSA (Invitrogen, 15 260-037) followed by myelin removal using 30% Percoll gradient (Sigma, P1644) at 750 Â g for 10 min at 4°C with slow brake. Cells were resuspended in fluorescence-activated cell sorting (FACS) buffer (1% BSA, 2 mM EDTA, 25 mM HEPES; Sigma, H0887) and passed through 100-mm (Corning, 431752) and 40-mm (Corning, 431750) cell strainers. Papain dissociation was used in favor of approaches such as mechanical dissociation and/or collagenase since we observed a superior cell survival when using Papain. The extent to which Papain affects the transcriptome of the cells is difficult to assess but Papain was used for all treatment groups. Dead cells and debris were removed using MACS Dead Cell Removal kit (Miltenyi Biotec, 130-090-101) using manufacturer's protocol. Cells were collected and resuspended in 100-ml FACS buffer and incubated in Fc-lock (CD16/32; BD, 553141) 1 mg/10 6 cells for 5 min on ice. Preconjugated antibodies (Table 2) were added onto the cell suspension at a concentration of 1 mg/10 6 cells and the cells were stained for 30 min on ice while protected from light.

Sorting single cells from spinal cord
Following staining the cells were washed in FACS buffer and eventually resuspended in 500 ml FACS buffer; 0.5-ml live/dead dye (Thermo Fisher Scientific, L34959) was added to the suspension. FACS was conducted using a Sony SH800S cell sorter. A mixture of unstained cells from healthy and injured spinal cord in equal volumes was used as negative control when setting compensations. CD45 1 cells from the spinal cords were sorted into Hard-Shell 384-well PCR plates (Bio-Rad, hsp3805) containing 2-ml lysis buffer in each well (Tables 3, 4). Plates were covered with a seal (Bio-Rad, MSF1001), centrifuged at 4000 rpm (Hettich, Universal 320R) for 2 min, then snap-frozen on dry-ice and stored at À80°C.  Library preparation and sequencing of single cells Full-leng cDNA and sequencing libraries were prepared using the Smart-Seq2 protocol (Picelli et al., 2014). Briefly, primers were annealed at 72°C for 3 min (Bio-Rad, 100°C Touch Thermal Cycler) followed by immediate reverse transcription and template switching by adding 3-ml (Eppendorf, epMotion 5073) reverse transcription mix (Tables 3, 4) to each well and kept at 42°C for 90 min and 70°C for 5 min. DNA was preamplified by adding 7.5 ml of preamplification mix (Tables 3, 4) into each well and kept at 37°C for 30 min, 95°C for 3 min, 22 cycles of: 98°C for 20 s; 67°C for 15 s; 72°C for 4 min followed by 72°C for 5 min. DNA was purified by adding solid phase reversible immobilization (SPRI) beads at a ratio of 0.75:1, incubated for 5 min at room temperature (RT), washed twice with 30 ml 80% EtOH (Solveco, 1000), air dried for 15 min and eluted in 15 ml EB buffer (QIAGEN, 19086). A magnetic plate (custom made) was used for separation of the beads during processing. DNA concentration was measured in 16 random wells from each 384-well plate using a step-wise pattern with QUBIT dsDNA HS Assay kit (Q33231); 3-ml sample was added to 97-ml buffer in a 96-well plate (Invitrogen, M33089) and concentrations were measured and compared with standard (Tecan spark 10 M, ex/em 485/528 nm). DNA was then diluted to a concentration of 0.15 mg/ml using RNase-free H 2 O. Tagmentation was performed by adding 1.8 ml Tn5 mix to each well in a 384-well plate (Table 3); 0.7 ml cDNA was then added into each well containing the Tn5 mix. The plate was vortexed (Thermo Fisher Scientific, 88880017TS), centrifuged and kept at 55°C for 10 min. Tn5 was inactivated by adding 1 ml 0.1% SDS (Sigma, 71725) to each well. Plate was vortexed and incubated at RT for 7 min. Libraries were barcoded by adding 19.5 ml PCR/bar code mix into each well now containing 3.5-ml sample (Table 3). A 2-ml primer (3.75 mM/primer) was finally added to each well. PCR/barcoding was then conducted at 72°C for 3 min, 95°C for 3 min followed by 12 cycles of: 95°C for 15 s; 55°C for 30 s; 72°C for 45 s. DNA was pooled (5 ml from each well) into an Eppendorf tube. Cleanup, as described above, was repeated and DNA was eventually eluted in 30-ml EB buffer. Concentration was again measured using QUBIT (3.0 Fluorometer) before 37 bp paired-end sequencing with dual index 8 bp using 2.6 pm on a 309 bp peak on Illumina NextSeq.

scRNAseq data analysis
Demultiplexing and alignment Demultiplexed FASTQ files were trimmed, adaptor sequences removed and aligned to the Genome Reference Consortium Mouse Build 38 (GRCm38, build mm10) using STAR (Dobin et al., 2013). Following removal of duplicate reads HT-seq was used for obtaining raw transcript counts. Eventually, the count matrix was exported into R (statistical computing software, v3.6.0; Team} RC, 2016) for downstream analysis (Table 5).
Data preprocessing and quality control Quality control wells (bulk and empty wells), ERCC, non-genes and genes present in zero genes were removed. Normalized counts per million (CPM) was estimated for each gene and cell by dividing the counts for each gene and cell with the sum of all counts for that particular cell, multiplying by 10 6 and adding 1 (pseudocount) before log2-transformation. Cells with Actb expression in the first percentile were removed (Extended Data Fig. 1-1A). Log2 of features was plotted against log2 of counts to determine the cutoff for cells with low number of counts and features. The minimum counts per cell was set at 2 16 (Extended Data Fig. 1-1B). The filtered data had a median sequencing depth of 606,166 reads/cell (IQR: 643,913) and 2191 genes/cell (IQR: 1228). Out of 4734 sequenced cells 3687 cells and 20,469 genes with a mean of 307 cells/condition remained for downstream analysis. No significant batch effects between samples sorted on different occasions were observed.
Clustering and dimensionality reduction R (Team} RC, 2016) package Seurat (v3.0.0) was implemented for downstream analysis of the preprocessed data. Normalization was conducted with scale.factor = 10 4 using NormalizeData(). Highly variable genes were defined as having a mean log2(cpm).0.3 (mean.Cutoff) and normalized dispersion (dispersion.Cutoff).0.5 keeping other parameters in FindVariableFeatures () at default (Extended Data Fig. 1-1C). Data were scaled using ScaleData() while RunPCA() was used for determining and detecting the first 50 principal components (PCs). The standard deviation for the 50 PCs were plotted to determine the highly variable PCs (Extended Data Fig. 1-1D). The first 17 PCs were selected for construction of the shared nearest neighbor (SNN) graph using FindNeighbors(). The SNN graph was used for identification of clusters using FindClusters() implementing the Louvain algorithm using resolution = 1.3. Resolution was set at 1.3 following evaluation of clustering using resolutions between 0.5 and 4. Hallmark genes for immune cells were used for evaluation of the clusters. A proper resolution was deemed to be one which did not overcluster or undercluster the data. Dimensionality reduction was conducted using RunTSNE(). The first 15 dimensions were estimated with seed.use = 132 and perplexity = 70. When including immune cells sorted following microglial depletion (plate IDs VZC01001, VZC01101) the key settings mentioned above were adjusted to: FindNeighbors (dims = 1:20), FindClusters (resolution = 0.75), RunTSNE (dims = 1:16, seed.use = 246), while the other settings were as above.
Pseudotime ordering of microglial cells Pseudotime analysis was only conducted for microglial cells (bMicroglia, aMicroglia, pmMicroglia, maMicroglia, DAM). Centroids (mean tSNE-1 and mean tSNE-2) were estimated for each condition. The Euclidian distance was estimated from each cell to each centroid, each cell being assigned to the closest centroid. Cells were ordered along pseudotime by estimating a principal curve using the centroids as anchors. A linear model was estimated for all ;20,000 genes. The 10 genes with the most positive slopes and the 10 genes with the most negative slopes were included in A heat map. The genes were ordered vertically using hierarchical clustering. log2(cpm) expression was reported for each microglial cell over pseudotime for a selected set of genes.
Cluster marker genes and differential gene expression analysis Marker genes for each cluster were identified using FindAllMarkers() using logfc.threshold=log(2) and other settings at default. Differentially expressed genes between selected clusters were identified using the edgeR (Robinson et al., 2010) and limma (Ritchie et al., 2015) packages in R (Team} RC, 2016) using the trimmed mean of M values (TMM) method for calculation of normalization.

Integration with peer datasets
The Seurat R package (v3.0.0) was implemented for integrating the current data with the three published datasets (Table 6). Each dataset was normalized using NormalizeData() and variable features were identified  using FindVariableFeatures() using selection.method=vst and nfeatures = 2000. Anchors were identified for the joint dataset using FindIntegrationAnchors() and the data were integrated using IntegrateData () both functions using the 30 first dimensions. Data were scaled using ScaleData() and RunPCA() was used to determine and detect the first 50 PCs with seed.use = 42. The first 10 PCs were selected for construction of the SNN graph using FindNeighbors(). The SNN graph was used for identification of clusters using FindClusters() implementing the Louvain Figure 1. A, Severe contusion SCI was induced in the thoracic spinal cord of adult female wild-type mice. CD45 1 immune cells from the SCI epicenter were isolated using FACS and sequenced individually following library preparation using the Smart-Seq2 protocol. Thin sections from mouse and human spinal cord were subjected to histologic analysis. B, tSNE plot of 3069 individually sequenced CD45 1 immune cells from SCI. Healthy spinal cord tissue was used as reference. Color of the cell corresponds to the time point following SCI at which the spinal cord was harvested and evaluated. C, tSNE plot of 3069 CD45 1 immune cells from mouse SCI. Eighteen clusters were detected and annotated using hallmark genes. bMicroglia: baseline microglia; aMicroglia: activated microglia; pmMicroglia: proliferation-mediating microglia; maMicroglia: monocyte-activating microglia; DAM: disease-associated microglia in SCI; elMicroglia: embryonic-like microglia, pvMacrophages: perivascular macrophages; mdMacrophages: monocyte-derived macrophages; DCs: dendritic cells; imB-cells: algorithm using resolution = 1. Dimensionality reduction was conducted using RunTSNE() using the first 10 dimensions with seed. use = 246 and perplexity = 70. In the first part of the analysis, the current data were integrated with data published by Li et al. (2019), while all four datasets, adding data from Masuda et al. (2019) and Mathys et al. (2017), were integrated at once in the end. The same settings were used during both integrations except in FindNeighbors() which was set to use the first 13, rather than 10, dimensions during the integration including four datasets.

Correlogram of transcription factors
Transcription factors defining each microglial cell cluster were identified using FindAllMarkers() using a minimum of 10% difference in fraction between clusters (min.diff.pct), logFC threshold at 0.4 (logfc.threshold) and genes that were detected in at least 30% of the cells in each cluster (min.pct), respectively. Top transcription factors for each cluster were selected based on logFC [fulfilling false discovery rate (FDR) , 0.05], and a correlogram was constructed indicating the correlation and cell type for each transcription factor.

Index sorting and FACS trace-back
Surface expression data were exported from the FACS apparatus. Because of slight variations in voltage settings between sorting sessions the expression data were separately normalized for each marker and plate. This allowed us to plot the expression for all cells at once in the same tSNE. FACS trace-back is defined as plotting all sorted and sequenced cells in a reconstructed FACS-plot (CD45 vs BSC-A) containing all analyzed cells to indicate granularity and CD45 expression of different cell types.
Depletion of Cx3cr1 1 cells Tamoxifen (Sigma, T5648) was resuspended in corn oil (Sigma, C8267) and heated to 75°C in a water bath and kept at 75°C for 60 min. The Cre-recombinase in CX3CR1 CreER/1 Rosa26 DTA/1 mice was activated by injecting 5 mg (200 ml) tamoxifen subcutaneously once daily for a total of 3 d. The depletion was evaluated and confirmed at various time points. Therefore, animals treated with tamoxifen were killed and the spinal cords harvested as described above. The spinal cords were stained and surface marker expression analyzed using flow cytometry (Beckman Coulter, Gallios). Data were analyzed using Kaluza software (Beckman Coulter). A piece of each spinal cord was preserved for histologic evaluation and was evaluated as described above.

Evaluation of hind limb locomotor function
The recovery of hindlimb locomotor function following SCI was evaluated using the Basso mouse scale (BMS; Basso et al., 2006) before, immediately following, and then weekly post-SCI until 90 d following SCI. Two researchers discussed and assessed the score. The mice were evaluated in a custom-made open-field environment constructed using Plexiglas.  Bulk RNAseq of microglia/macrophages FACS Microglia/macrophages were defined as CX3CR1 1 CD19 -Ly6G -CD3 -NK1.1 -SIGLECF -CD11b 1 CD45 high/low cells and sorted into FACS buffer using a BD Influx cell sorter.

Isolation of total RNA
Sorted microglia/macrophages in FACS buffer were collected by centrifugation at 300 Â g for 5 min and resuspended in 1-ml TRIzol reagent (Thermo Fisher Scientific, 15596026) and stored at À70°C until downstream analysis. RNA was isolated using the manufacturer's (TRIzol) protocol. Isopropanol precipitation was conducted using 75-mg glycogen (GlycoBlue Coprecipitant, Thermo Fisher Scientific AM9516) at À20°C O/N. RNA clean-up was performed using RNeasy micro kit (QIAGEN, 74004) according to the manufacturer's protocol. Contaminating genomic DNA was removed during the isolation by on-column digestion with DNase (DNase I QIAGEN, 79254). RNA was stored at À70°C.
Library preparation, sequencing, and differential gene expression analysis Total RNA from microglia/macrophages was sequenced. Sequencing libraries were prepared using SMARTer Stranded Total RNA-Seq kit-Pico Input Mammalian (Clontech) because of the low amount of RNA. Libraries were sequenced 2 Â 125 bp on three lanes using HiSeq2500 (Illumina Inc.) on high output mode. For each replicate a minimum of 15 Â 10 6 read-pairs were used. TrimGalore (Babraham Bioinformatics) was used for removal of low quality regions and adapter sequences. Read-pairs were aligned to the mouse genome (GRCm38, build mm10) using STAR (Table 5). featureCounts (Liao et al., 2014) and Ensembl annotation release 81 was used for summary of read counts over genes. Differential gene expression analysis was conducted in R (v3.6.0; Team} RC, 2016) using limma (Ritchie et al., 2015) and edgeR (Robinson et al., 2010) with annotations from Mus. musculus (Team BC, 2015).

Proliferation analysis
Spinal cords subjected to SCI were dissociated and stained with antibody against CD45 as described above. Following wash, the cells were incubated in 1:2000 Hoechst 33342 for 15 min, washed twice in 1Â DPBS and resuspended in 0.5 ml 1Â DPBS. Live/dead marker was added to the suspension in the same manner as described above (Table 2). Surface expression data were collected and analyzed offline. The CD45 1 cells were divided into CD45 low and CD45 high . For each subset (low/high) the expression of Hoechst 33342 (Thermo Fisher Scientific, 62249) was plotted in a histogram. The cells with the highest Hoechst 33342 expression were defined as cells in S-phase of the cell cycle and therefore proliferative cells.

Statistical analysis
Data for hindlimb function, proliferation rate and number of microglia following depletion was presented for each biological replicate in combination with a mean surrounded by a 95% confidence interval (CI). Hind limb function was further investigated using a mixed ANOVA followed by Wilcoxon post hoc test between groups within each time point. Correlation between current data and peer data were conducted using Pearson's product moment correlation coefficient.

Results
Microglia undergo rapid and profound changes on SCI, ultimately arresting in a chronic DAM-like state We characterized the immune cell population over time following severe traumatic contusion SCI in mice. CD45 1 cells were isolated from the injury epicenter (site of impact 63 mm along the cranial-caudal axis), and libraries for each individual cell were prepared using the Smart-Seq2 protocol (Picelli et al., 2014). To ensure that the time-points were comparable and to minimize the impact of random variation, a pooling strategy was implemented in which a suspension of cells collected from spinal cords of three to five mice were used for each time point. The immune cell response was investigated at nine time points following SCI, spanning both the acute and chronic phase, and using healthy spinal cord as reference (Fig. 1A). A total of 3069 individually sequenced immune cells were included following quality control (Extended Data Fig. 1-1A-G; Table 7). The cells organized into distinct clusters, with various degrees of temporal influence (Fig. 1B), whereas the transcriptional profile of B and T cells remained largely constant over time, microglia underwent a dramatic transformation following SCI, and never reestablished their original transcriptional profile at any time point following SCI (Extended Data Fig. 1-2A,B). Five subtypes of microglia were identified ( Fig. 1C; Extended Data Fig. 1-2B). Baseline microglia (bMicroglia) defined by expression of P2ry12, Tmem119 and Siglech underwent a rapid transformation into an activated form of microglia (aMicroglia), defined by expression of inflammation-associated markers such as Il1a and Cd83, already at 0.5 h post-SCI (Fig. 1C-E). aMicroglia had fully transformed into a distinct population of proliferation-mediating microglia (pmMicroglia) defined by expression of the growth factor Hbegf as well as Edn1 at 2 h post-SCI. The transformation halted between 6 and 36 h in a state  Cells were ordered along pseudotime by estimating a principal curve using centroids for each condition as anchors. Centroids were estimated by calculating the mean tSNE-1/2 value for all cells within each condition. Each cell was then assigned to a centroid based on a minimum Euclidean distance approach. B, Density of 1474 microglial cells in each centroid along pseudotime. C, Heat map reporting log2(cpm) expression for selected genes for 1474 microglial cells ordered along pseudotime. A linear model was estimated for all ;20,000 genes. The 10 genes with the most positive slopes and the 10 genes with the most negative slopes were included in the heat map. The genes were ordered vertically using hierarchical clustering. Rows represent genes while columns represent cells. D, log2(cpm) expression for microglial cells for selected genes arranged along pseudotime. A locally estimated scatterplot smoothing (loess) was estimated for each gene and reported using a red dotted line in combination with a 95% CI. Each dot represents one cell. E, tSNE plots of 3069 CD45 1 immune cells reporting the log2(cpm) expression with a continuous color for selected genes. F, Protein expression of APOE on mouse spinal cord sections. Scale bars: 50 mm.  Table 6). C, D, Each dataset represented by one tSNE. Cells are annotated using color representing their cell type. bMicroglia: baseline microglia; aMicroglia: activated microglia; pmMicroglia: proliferation-mediating microglia; maMicroglia: monocyte-activating microglia; DAM: disease-associated microglia in SCI; elMicroglia: embryonic-like microglia; pvMacrophages: perivascular macrophages; mdMacrophages: monocyte-derived macrophages; DCs: dendritic cells; imB-cells: immature B-cells; mB-cells: mature B-cells; pimMicroglia: postnatal immature microglia; IEG1Microglia: immediate early genes microglia; sMicroglia: sphase microglia; pPAM: postnatal proliferative-region-associated microglia; eMicroglia: embryonic microglia; pcpMacrophages: postnatal choroid plexus macrophages; acpMacrophages: adult defined by expression of Mt2 and Msr1 resulting in monocyteactivating microglia (maMicroglia). DAM-like cells appeared at 3 d post-SCI, and increased in number from 3 to 21 d post-SCI when it reached a steady state which persisted during the chronic phase of SCI (90 d post-SCI; Fig. 1C-E; Extended Data Fig. 1-2A). These cells bear all the hallmarks of DAM described in mouse models of neurodegeneration, demyelination and development (Wang et al., 2015;Keren-Shaul et al., 2017;Krasemann et al., 2017;Mathys et al., 2017;Friedman et al., 2018;Mrdjen et al., 2018;Olah et al., 2018;Sousa et al., 2018;Tay et al., 2018;Hammond et al., 2019;Jordão et al., 2019;Li et al., 2019;Masuda et al., 2019;O'Koren et al., 2019;Figs. 1D, 2A). A distinct set of transcription factors could be determined for each microglial subtype, with Egr3 and Skil associated with aMicroglia, pmMicroglia mainly associated with Bach1 and Etv3, while maMicroglia were associated with transcription factors Id2 and Myc. Glmp and Nfe2l2 were transcription factors unique for DAM in SCI, while Fosb was expressed in both DAM and bMicroglia. Similarly, Skil was expressed in both aMicroglia and bMicroglia (Fig. 2B). Protein expression of CST7, IGF1, and APOE was confirmed using immunohistochemistry on thin sections of injured spinal cord from both mice and humans (Figs. 2C, 3-6). Microglia and DAM could not be distinguished using surface expression of CD45 alone (Extended Data Fig. 1-3), but differential expression of the genes defining the two cell types could be detected in bulk RNAseq (Extended Data Fig. 1-4).
Pseudotime analysis of the microglial populations confirmed our initial classification (Fig. 7A), with bMicroglia undergoing a rapid transformation immediately following SCI (0-2 h), adopting an intermediate state between 6 and 36 h and gradually falling into the steady-state DAM profile which persisted in the chronic phase of the injury (90 d post-SCI; Fig. 7B). Genes defining bMicroglia were downregulated over pseudotime, while the opposite was true for genes defining DAM in SCI (Figs. 2A, 7C-E). Protein expression of CST7, IGF1, and APOE concorded with the gene expression, both in mouse and human tissue sections (Figs. 2C, 3-6, 7F).
Taken together, our results show a rapid transformation of microglial cells following injury, in combination with a loss of baseline markers and gradual adoption of a permanent disease-associated cell state. DAM in SCI appear to be low-proliferative microglial cells defined by genes related to lipoprotein and cholesterol metabolism, lysosomal proteolysis as well as cell and organelle adhesion. In combination with the activation of distinct transcription factor cascades, this suggests transcriptional reprogramming of bMicroglia into DAM, via a previously unrecognized intermediate state during the first week post-SCI. Interestingly, this intermediate state was more distinct from bMicroglia than chronic DAM. DAM in trauma, demyelination, degeneration, and development of the CNS are transcriptionally similar Based on marker gene expression, we found that DAM in SCI bore a striking similarity to DAM found during neurodegeneration, demyelination and development. In order to more coherently compare our cells to those found in other disease processes, we used molecular anchor-based integration of our dataset with published datasets. Li et al. (2019) sequenced 2123 microglia and myeloid cells from mouse brain during three different stages of development (Table 6), resulting in a dataset of 4659 immune cells prepared using similar methodology to our work (Fig. 8A,B). The authors detected postnatal proliferative-region-associated microglia (pPAM), a subtype of microglia with striking similarities to DAM, embedded in the white matter of developing mouse brains. pPAM were metabolically active with an amoeboid morphology and phagocytized newly formed oligodendrocytes. As expected, DAM in SCI co-localized with pPAM in the aggregated tSNE (Fig. 8C,D). aMicroglia in SCI co-localized with immediate early genes (IEG)1 Microglia described by Li et al., while cells expressing high levels of C1qa-c overlapped with embryonic microglia (eMicroglia) and were therefore termed embryonic-like microglia (elMicroglia; Fig. 8E). The proliferative cluster (G2/ M-microglia) found in embryonal tissue, was however missing from our data, which is consistent with our observed low rate of proliferation (;1.0-1.5% of all CD45 low cells were in Sphase). No significant difference in proliferation was detected at 7 and 21 d post-SCI (Fig. 8F). Integration of datasets thus confirmed the similarities between microglial subtypes in traumatic injury to a subtype of microglia found in developing brain.
We similarly analyzed two published datasets collected from models of neurodegenerative disease in the CNS (10,014 cells in total; Table 6). Masuda et al. (2019) used scRNAseq to investigate microglia in different regions of the brain when using a model of toxic demyelination (cuprizone) and neurodegeneration (unilateral n.VII axotomy). We integrated 3721 microglia harvested from healthy tissue as well as from brain and spinal cord tissue subjected to n.VII axotomy or demyelination into our dataset (Extended Data Fig. 8-1A). In the aggregated tSNE microglia during demyelination co-localized somewhat with DAM in SCI, while microglia following n.VII axotomy (3/14 d after injury) did not co-localize (Extended Data Fig. 8-1B,C). Microglia harvested from healthy 3-and 16-week-old mice were located in proximity to, but did not overlap with bMicroglia. However, embryonal microglia (healthy) formed a distinct cluster in proximity to DAM and pPAM. Thus, microglia during injury seem to have a transcriptional profile similar to embryonal and postnatal microglia. The differences between bMicroglia and microglia at 3 and 16 weeks might be explained by Masuda et al. (2019): (1) using CD-1 mice as compared with C57BL/6J mice in this study; (2) using both brain and spinal cord tissue as compared with spinal cord only in this study; (3) gating, sorting and sequencing selected microglia (CD45 int CD11b 1 Ly6C -Ly6G -CD206 -) as compared with the current study which investigated all CD45 1 cells/microglia and d) evaluating microglia during de/remyelinzation in the corpus callosum only. Mathys et al. (2017) isolated and sequenced 1634 microglial cells from the hippocampus of mice with neurodegeneration resulting in symptoms of AD (CK-p25 inducible mice). The microglial cells from these CK-p25 mice broadly matched our DAM, but no clear subtype co-localization was detected (Extended Data Fig. 8-1A-C). Correlation analysis using Pearson's product moment correlation coefficient between differential gene expression for DAM in SCI and pPAM and microglia during AD (Mathys et al., 2017), respectively, confirmed the downregulation of baseline genes and upregulation of disease-associated genes (Extended Data Fig. 8-1D). Taken together, our analysis strongly supports that DAM in trauma,  Microglia depletion refers to CD45 1 immune cells isolated from mice subjected to microglial depletion before SCI. C, tSNE plot of 3687 CD45 1 immune cells clustered and annotated according to cell type. bMicroglia: baseline microglia; aMicroglia: activated microglia; pmMicroglia: proliferation-mediating microglia; maMicroglia: monocyte-activating microglia; DAM: disease-associated microglia in SCI; elMicroglia: embryonic-like microglia, pvMacrophages: perivascular macrophages; mdMacrophages: monocyte-derived macrophages; DCs: dendritic cells; imB-cells: immature B-cells; mB-cells: mature B-cells. D, Bar plot reporting the percentage portion which each cell type makes up out of all cells at each time point of evaluation. E, tSNE plot of 3687 CD45 1 immune cells reporting log2(cpm) expression for a set of selected genes using a continuous color for each cell separately. F, Volcano plot reporting differentially expressed genes between mdMacrophages and bMicroglia. G, BMS score over time. Each dot represents one animal. Mean value for each time point and study group is surrounded by a 95% CI; *p , 0.05; ns, not significant. demyelination and development are similar both with respect to expression of effector genes, whole transcriptional profile, and morphology (where investigated).
DAM in SCI are permanently reprogrammed baseline microglia which enhance functional recovery Having established the identity of DAMs in the chronic phase SCI response, two important questions remain. First, do DAM cells originate from resident baseline microglia or peripheral infiltrating cells? Our pseudotime analysis strongly suggest that bMicroglia gradually change states to transition into DAM cells, but does not formally exclude a peripheral origin from tissue or bone marrow reservoirs. Second, do DAM cells impact functional recovery following SCI? Previous studies have used disease models in which the evaluation of the function is difficult, in contrast to gross motor evaluation in a SCI model.
To address these questions, we used a transgenic mouse model (CX3CR1 CreER/1 Rosa26 DTA/1 ), which allows tamoxifendependent conditional depletion of Cx3cr1 1 cells. Cx3cr1 is mainly expressed by CNS microglia but is also expressed by monocytes, dendritic cells and natural killer cells (The Jackson Laboratory). The lack of microglia and monocytes in the CNS at the moment of SCI allowed us to investigate the origin of DAM as well as their contribution to function. Hence, a lack of DAM following SCI in Cx3cr1-depleted mice would suggest that DAM are indeed reprogrammed bMicroglia. However, if DAM would appear this would suggest that DAM are derived from peripheral monocytes entering the CNS following SCI (Fig. 9A). As determined by flow cytometry and immunohistochemistry, depletion was almost complete, with the percentage of CD11b 1 CX3CR1 1 or CD11b 1 CD45 low at a nadir (;0.5% of cells) 7 d after tamoxifen administration (Fig. 10A-C). We isolated and sequenced 604 CD45 1 immune cells at 7 and 21 d post-SCI (14 and 28 d after full depletion, respectively). In the microglia-depleted mice, we observed a total absence of the previously observed bMicrogliaderived cells, including DAM (Fig. 9B,C). We instead observed a large increase in monocytes and monocyte-derived macrophages (mdMacrophages) expressing Ms4a7 (Fig. 9B-E), instead of genes characteristic for DAM in SCI (Fig. 9F). Thus, when no bMicroglia were present in the spinal cord at the moment of injury DAM did not appear, in line with bMicroglia being the cell of origin for DAM.
To determine whether microglia-derived cells had an effect on recovery from injury, we concomitantly measured hindlimb locomotor function during 90 d post-SCI using the BMS (Fig.  9A). Over time the BMS score was significantly lower (p , 0.05, mixed ANOVA combined with Wilcoxon post hoc test) among animals subjected to depletion of Cx3cr1 1 cells before SCI (1.83, CI: 1.51-2.16) as compared with wild-type mice (3.25, CI: 2.60-3.90) also subjected to SCI (p , 0.05; Fig. 9G). Cx3cr1 1 cell-depleted animals which were not subjected to SCI did not deteriorate in terms of BMS score and remained at a full healthy score (9) over time (data not included). The majority of the difference in BMS score between depleted and wild-type mice was observed during the first week post-SCI, concurrent to DAM transformation. Taken together, our microglia depletion experiment shows that DAM in SCI are formed through transcriptional reprogramming of bMicroglia and have a strong beneficial effect on the recovery of hindlimb locomotor function post-SCI.

Discussion
The present study provides a detailed description of the temporal dynamics of immune cells at the injury epicenter of traumatic SCI, at single-cell resolution. We find that following trauma to the CNS, baseline microglia undergo a distinct temporal transformation through cellular reprogramming to assume a diseaseassociated phenotype similar to cell states found in neurodegeneration, demyelination and development. This subtype persists during the chronic phase of disease in both rodents and humans, and contributes to functional recovery.
Baseline microglia (bMicroglia), found in healthy animals, transformed into DAM in two distinct stages with different kinetics, first transitioning within hours to an intermediary activated form (aMicroglia) and reaching a DAM-like state beginning 3 d post-SCI. The DAM state appears permanent since bMicroglia did not reappear within the 90-d period, indicating that they have undergone stable cellular reprogramming. Such stable reprogramming has previously been reported as a diseaseassociated transformation of microglia in both chronic and acute neurodegenerative disease conditions such as AD, demyelination, or inflammation/cell death (Keren-Shaul et al., 2017;Krasemann et al., 2017;Mathys et al., 2017;Sousa et al., 2018;Tay et al., 2018;Jordão et al., 2019;Masuda et al., 2019). In addition, a similar DAM-like population in microglia also exist in healthy developing mouse brain (Li et al., 2019). This suggests a canonical microglial activation pathway with cell death, by controlled apoptosis or injury, as the underlying trigger rather than disease-specific processes.
The origin of DAM, which appear following pathology and/ or injury to the CNS, is not known. DAM could either be reprogrammed bMicroglia or the progeny of peripheral immune cells entering the CNS through the blood brain/spinal cord barrier. Our results clearly indicate that DAM are permanently reprogrammed bMicroglia, the bMicroglia cells undergo a rapid and gradual transformation into DAM, consistent with cell state transition, and depletion of bMicroglia completely abolished DAM formation. In addition, the rapid transformation of bMicroglia into aMicroglia are in line with transcriptional reprogramming rather than the slower process of proliferation of infiltrating cells. These findings are in line with Plemel et al. (2020) who found that microglia, rather than infiltrating macrophages, dominate the spinal cord following demyelination. They also found that absence of microglia during demyelination resulted in a significantly higher infiltration of peripheral macrophages. Charting the molecular mechanisms of microglial reprogramming may reveal ways of modulating their function for therapeutic purposes, promising to aid recovery in traumatic or neurodegenerative disease progression.
Our microglia depletion experiment shows that microglial activation has a strong positive effect on recovery from SCI. This is in line with previous studies which have reported that pharmacological depletion of microglia (using a CSF1R inhibitor) following SCI results in impaired functional recovery (Bellver-Landete et al., 2019). However, the fully formed DAM might not be the main effector of recovery. Whereas microglia rapidly transform during the first 2 h, they halt in an intermediate state / singlets in the specified gates are reported in the plots. B, Protein expression of TMEM119 and CX3CR1 on representative thin histologic sections of spinal cord at each time point. C, Level of microglia in the spinal cord over time following injection of TAM presented as percentage of singlets in the CD11b 1 CD45 low and CD11b 1 CX3CR1 1 gates, respectively. Each dot represents one animal. Mean value at each time point of evaluation is surrounded by a 95% CI. at 6-36 h (Fig. 7B). During this time the cells express low levels of baseline microglial markers and slowly upregulate expression of genes associated with DAM. Considering that the majority of the functional recovery is observed already at 7 d post-SCI (i.e., before full maturation of DAM), it is plausible that this intermediate population mediates the majority of the beneficial effects resulting in enhanced functional recovery. We hypothesize that DAM (and/or intermediate populations) contribute to functional recovery by beneficially modulating the immune response following SCI. Removal of debris and dead cells as well as an appropriate inflammatory response modifying immune cell and astrocyte response is essential to minimize the glial scar and thereby promote axonal regeneration. Previous studies have mainly investigated neurodegenerative disease models (Multiple sclerosis, AD) in rodents, in which functional recovery (e.g., sensation, autonomic functions, memory) is more difficult to evaluate (Keren-Shaul et al., 2017;Krasemann et al., 2017;Mathys et al., 2017;Friedman et al., 2018;Mrdjen et al., 2018;Olah et al., 2018;Sousa et al., 2018;Tay et al., 2018;Hammond et al., 2019;Jordão et al., 2019;Masuda et al., 2019). The traumatic SCI model allows for exact and detailed evaluation of motor function, which has not previously been evaluated in a context of SCI/neurotrauma. The results presented in this study are consistent with, that DAM being beneficial for functional recovery. However, since Cx3cr1 is expressed by other immune cells in addition to microglia, we cannot exclude the possibility that for example dendritic cells or natural killer cells are the main factor beneficial for motor recovery. Irrefutable proof of the role of DAM in recovery will have to await the development of a fully specific depletion model.
Taken together, we have identified a disease-associated microglial cell type which originates from permanently reprogrammed baseline microglia via an activated cell state. These DAM cells appear within days following traumatic injury, and persist years later in human spinal cord at the site of injury. These cells seem to contribute to healing, and show a striking similarity to DAM found in neuronal degeneration or demyelination, as well as a specific microglial cell type found during development.