Single cell atlas of spinal cord injury in mice reveals a pro-regenerative signature in spinocerebellar neurons (original) (raw)

Nat Commun. 2022; 13: 5628.

Kaya J. E. Matson,1,2 Daniel E. Russ,3 Claudia Kathe,4,5 Isabelle Hua,1 Dragan Maric,6 Yi Ding,7 Jonathan Krynitsky,8 Randall Pursley,8 Anupama Sathyamurthy,1,9 Jordan W. Squair,4,5 Boaz P. Levi,7 Gregoire Courtine,4,5 and Ariel J. Levine corresponding author 1

Kaya J. E. Matson

1Spinal Circuits and Plasticity Unit, National Institute of Neurological Disorders and Stroke, National Institutes of Health, Bethesda, MD USA

2Johns Hopkins University Department of Biology, Baltimore, MD USA

Daniel E. Russ

3Division of Cancer Epidemiology and Genetics, Data Science Research Group, National Cancer Institute, NIH, Rockville, MD USA

Claudia Kathe

4Center for Neuroprosthetics and Brain Mind Institute, Faculty of Life Sciences, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland

5NeuroRestore, Department of Clinical Neuroscience, Lausanne University Hospital (CHUV) and University of Lausanne (UNIL), Lausanne, Switzerland

Isabelle Hua

1Spinal Circuits and Plasticity Unit, National Institute of Neurological Disorders and Stroke, National Institutes of Health, Bethesda, MD USA

Dragan Maric

6National Institute of Neurological Disorders and Stroke, Bethesda, MD USA

Yi Ding

7Allen Institute for Brain Science, Seattle, WA USA

Jonathan Krynitsky

8Signal Processing and Instrumentation Section, Center for Information Technology, National Institutes of Health, Bethesda, MD USA

Randall Pursley

8Signal Processing and Instrumentation Section, Center for Information Technology, National Institutes of Health, Bethesda, MD USA

Anupama Sathyamurthy

1Spinal Circuits and Plasticity Unit, National Institute of Neurological Disorders and Stroke, National Institutes of Health, Bethesda, MD USA

9Present Address: Centre for Neuroscience, Indian Institute of Science, Bangalore, India

Jordan W. Squair

4Center for Neuroprosthetics and Brain Mind Institute, Faculty of Life Sciences, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland

5NeuroRestore, Department of Clinical Neuroscience, Lausanne University Hospital (CHUV) and University of Lausanne (UNIL), Lausanne, Switzerland

Boaz P. Levi

7Allen Institute for Brain Science, Seattle, WA USA

Gregoire Courtine

4Center for Neuroprosthetics and Brain Mind Institute, Faculty of Life Sciences, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland

5NeuroRestore, Department of Clinical Neuroscience, Lausanne University Hospital (CHUV) and University of Lausanne (UNIL), Lausanne, Switzerland

Ariel J. Levine

1Spinal Circuits and Plasticity Unit, National Institute of Neurological Disorders and Stroke, National Institutes of Health, Bethesda, MD USA

2Johns Hopkins University Department of Biology, Baltimore, MD USA

3Division of Cancer Epidemiology and Genetics, Data Science Research Group, National Cancer Institute, NIH, Rockville, MD USA

4Center for Neuroprosthetics and Brain Mind Institute, Faculty of Life Sciences, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland

5NeuroRestore, Department of Clinical Neuroscience, Lausanne University Hospital (CHUV) and University of Lausanne (UNIL), Lausanne, Switzerland

6National Institute of Neurological Disorders and Stroke, Bethesda, MD USA

7Allen Institute for Brain Science, Seattle, WA USA

8Signal Processing and Instrumentation Section, Center for Information Technology, National Institutes of Health, Bethesda, MD USA

9Present Address: Centre for Neuroscience, Indian Institute of Science, Bangalore, India

Corresponding author.

Received 2022 Jun 8; Accepted 2022 Aug 31.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Supplementary Materials

Supplementary Information

GUID: 2FBAA7D9-5DB2-4296-B142-9137358E62E4

Description of Additional Supplementary Files

GUID: 58A852EA-3022-4BCD-B83F-5113038142D9

Supplementary Dataset 1

GUID: 732EC4D9-9BF9-4767-AC00-D3FB0A93386F

Supplementary Dataset 2

GUID: C1B04402-5357-41A9-A090-B8E7445578CE

Supplementary Dataset 3

GUID: 13C146BA-6204-4AF0-AC4D-A2610D6F67C2

Supplementary Dataset 4

GUID: A0B1447F-6E33-4B46-96E9-0AC1A959BE69

Supplementary Dataset 5

GUID: CFF63044-EF1C-48D2-8F39-A92800385612

Supplementary Dataset 6

GUID: 8060126C-CEB7-4D77-A2AC-E1BDF709091B

Supplementary Dataset 7

GUID: 92759F60-83BA-4588-92D7-265B5AF8DBCE

Reporting Summary

GUID: 81FACED0-F0FF-42D3-939E-5D98BE2E119D

Data Availability Statement

A searchable version of this data is available at https://seqseek.ninds.nih.gov/spinalcordinjury. Raw sequencing data and count matrices have been deposited to the Gene Expression Omnibus (GSE172167). Source data are provided with this paper.

Custom MATLAB-based code for quantification of cell counts is available at https://github.com/ArielLevineLabNINDS/CellCounter (10.5281/zenodo.6967482). Custom python-based code for quantification of pixels on thresholded images is available at https://github.com/ijhua/pixel_counts. Code for Augur is available at https://github.com/neurorestore/Augur.

Abstract

After spinal cord injury, tissue distal to the lesion contains undamaged cells that could support or augment recovery. Targeting these cells requires a clearer understanding of their injury responses and capacity for repair. Here, we use single nucleus RNA sequencing to profile how each cell type in the lumbar spinal cord changes after a thoracic injury in mice. We present an atlas of these dynamic responses across dozens of cell types in the acute, subacute, and chronically injured spinal cord. Using this resource, we find rare spinal neurons that express a signature of regeneration in response to injury, including a major population that represent spinocerebellar projection neurons. We characterize these cells anatomically and observed axonal sparing, outgrowth, and remodeling in the spinal cord and cerebellum. Together, this work provides a key resource for studying cellular responses to injury and uncovers the spontaneous plasticity of spinocerebellar neurons, uncovering a potential candidate for targeted therapy.

Subject terms: Molecular neuroscience, Spinal cord injury

Matson et al. performed single nucleus sequencing of the “spared” spinal cord tissue distal to an injury in mice. They found that spinocerebellar neurons expressed a pro-regenerative gene signature and showed axon outgrowth after injury.

Introduction

The brain, spinal cord, and peripheral nervous system are comprised of diverse cell types that operate together as global and local communities to enable normal physiology. Following acute trauma, a complex interplay of cellular responses shapes the outcome. Whether the tissue can restrict the damage, promote structural remodeling and functional compensation, and ultimately achieve recovery depends on a myriad of dynamic molecular changes amongst neurons, astrocytes, microglia, oligodendrocytes, vascular cells, and many other cell types1,2.

Spinal cord injury (SCI) is a traumatic event that can cause long-lasting paralysis, pain, autonomic dysregulation, and body-wide physiological changes3. While understanding and developing therapeutics that target cellular changes within the lesion epicenter is undoubtedly valuable, there is an emerging focus on molecular and neural engineering approaches to augment plasticity and reorganization in anatomically incomplete injuries, which are the most common in patients4–6. Molecular approaches can induce sprouting of descending spared projections7,8 or reorganization of the neural circuits below the injury4,6. Neural engineering approaches such as epidural electrical stimulation combined with rehabilitation training can promote impressive gains in motor function and autonomic control and provide enhanced quality of life9–11, underscoring the importance of understanding intrinsic potential in the tissue below the injury site. Understanding the cellular mechanisms by which reorganization occurs in spinal cord neurons is a crucial step to promote recovery.

What may be the underlying mechanisms of plasticity in the lumbar spinal cord that enable recovery after injury? There are many forms of plasticity, from synaptic remodeling to local axonal sprouting and long-distance axon growth12,13. After injury in the peripheral nervous system, damaged axons can successfully regrow and innervate their targets to restore function14,15. However, regeneration in the CNS is generally limited, due to both cell-intrinsic and cell-extrinsic factors16. Although most CNS neurons are capable of some plasticity, these processes are limited to local changes such as synaptic remodeling. However, there are a few examples of CNS neurons that can regenerate17, suggesting that regenerative capacity is a cell type-specific feature. We hypothesized that the many neuronal subpopulations in the lumbar cord—with their diverse molecular identities18,19, physical properties, and connectivity20, – may differentially contribute to recovery after spinal cord injury21, and that specific neuronal subpopulations may display their own strategies for repair.

Here, we sought to uncover the dynamic cell type-specific responses of the lumbar spinal cord following SCI to identify the cell type-specific molecular and cellular mechanisms that promote or restrict recovery. First, we performed severe thoracic contusion spinal cord injuries in mice and tracked the progression of injury responses from acute to chronic time points. To profile the diverse cell types within the lumbar spinal cord following thoracic injury, we used single nucleus RNA Sequencing (snRNA-seq) and created an atlas of the lumbar cell types after injury (https://seqseek.ninds.nih.gov/spinalcordinjury). This resource reveals both the changes in relative composition of cell types following SCI and the changes in gene expression within each cell type. The size and scope of this dataset allowed identification of rare cellular populations that displayed molecular pathways with the potential to support recovery. Specifically, we identified neuronal populations that expressed regeneration-associated genes (RAGs). These neurons were largely excitatory, and their spatial distribution, as well as gene expression, suggested that they are ascending projection neurons that link the spinal cord and brain. We identified the RAG-expressing neurons to be Shox2-expressing V2d and spinocerebellar neurons. Using viral-labeling strategies, we showed that after thoracic injury spinocerebellar neurons increased axons and collaterals below the injury site, indicating structural remodeling. Together these findings shed light on the limited spontaneous mechanisms of repair in the tissue below the lesion and the latent potential for targeted neuro-regeneration and tissue remodeling therapies.

Results

A single-cell atlas of the lumbar cord after spinal cord injury

We profiled the cell types in the lumbar cord after spinal cord injury at translationally relevant time points to create an atlas of cell type-specific responses and uncover biological changes. A severe contusion was delivered to the thoracic (vertebral level T9) spinal cord of mice, resulting in paralysis (Supplementary Fig 1, Supplementary Dataset 1). Locomotor function was tracked over a range of time points that, in mice, corresponded to the acute injury period (1 day post injury (dpi)), the subacute and intermediate stages that are typically targeted with therapeutic interventions (1 week and 3 weeks post injury/wpi), and a chronic time point at which recovery typically plateaus (6 wpi; Supplementary Fig 1). We then dissected the lumbar spinal cord of three animals from each time point and performed single nucleus isolation and RNA sequencing (Fig. 1a, b)22. We clustered the data, integrating by time point, and removed clusters that were low-quality as well as doublet clusters, yielding a final dataset of 67,903 nuclei (see “Methods” for details).

Single nucleus RNA sequencing of the lumbar spinal cord after thoracic contusion.

a Schematics depicting the experimental design for snRNA-seq, showing the injured thoracic cord and lumbar cord (with dark red representing the lesion) as well as nuclei isolation from the intact lumbar cord followed by droplet-based barcoding for single nucleus RNA sequencing. b Top, an overview of experimental design for injury and tissue collection. The lumbar spinal cord of three animals from each time point: uninjured, 1 dpi (day post injury), 1 wpi (week post injury), 3 wpi, and 6 wpi. c Uniform manifold approximation and projection (UMAP) visualization of 67,903 nuclei from uninjured and injured lumbar spinal cords, revealing 9 classes and 39 subtypes. Colored by green (neurons), yellow, astrocytes, orange-red (microglia), purple (OPCs), blue (oligodendrocytes), light blue (Schwann), light pink (pericytes), pink (ependymal), magenta (leptomeninges), brick-red (endothelial). d Multiplex immunohistochemistry (IHC) of the lumbar spinal cord from uninjured, 1 wpi and 3 wpi. Tissue was stained for NeuN (green), GFAP (yellow), IBA1 (red), TMEM119 (dark orange), and OLIG2 (blue). Scale bars are 200 µm (main) and 50 µm (inset), respectively. e Quantification of the proportion of cell types from the snRNA-seq data and immunohistochemistry in tissue. Mean ± SEM; snRNAseq, N = 3; immunohistochemistry, N = 4. Source data are provided as a Source Data file.

To create an atlas of cell type responses in the lumbar cord after injury, we first clustered the nuclei at a coarse level to highlight large-scale changes. We identified 9 major cell classes: neurons, astrocytes, microglia/hematopoietic cells, oligodendrocyte lineage cells, Schwann cells, endothelial cells, pericytes, ependymal cells, and leptomeninges (Fig. 1c). Each of these major classes was identified using well-established markers19 (Supplementary Fig. 2, Supplementary Dataset 2). Neurons expressed Snhg11, Rbfox1, Rbfox2, Snap25. Astrocytes expressed Slc7a10, Agt, Gfap, and Vim. Microglia/hematopoietic cells expressed C1qa, Ctss, Gpnmb, Lgals3, Itgax, Ms4a4b, Cd3g, and Nkg7. Oligodendrocyte lineage cells including oligodendrocyte precursor cells (OPCs) expressed Cspg5 and Tnr; committed oligodendrocyte progenitor cells (COPs) expressed Fyn and Tcf7l2, whereas myelinating and mature oligodendrocytes expressed Plp1, Mag, and Mog. Schwann cells, which were part of the lumbar spinal cord roots that were in the dissected tissue, expressed Mpz and Pmp22. Vascular cells included endothelial cells, which express Bsg and Cldn5, and pericytes, which express Vtn and Pdgfrb. Ependymal cells expressed Nnat and Dnah12. Leptomeninges expressed Dcn and Col1a1 (Supplementary Fig. 2, Supplementary Dataset 2). With this approach, we have identified the major cell classes of the uninjured and injured mouse spinal cord.

Coarse cell types in snRNA-seq and in tissue

We next compared the proportion of the coarse cell types in the snRNA-seq dataset to the proportion of these cell types as detected by immunohistochemistry in tissue sections. In the snRNA-seq data, we calculated the proportion of each cell type within a given sample. In the uninjured tissue, oligodendrocytes represented the largest proportion of the uninjured spinal cord (34.6 ± 1.4%), followed by neurons (29.4 ± 1.7%), astrocytes (10.5 ± 1.3%), and microglia/hematopoietic cells (3.7 ± 0.3%; mean ± SEM, N = 3, Fig. 1e). We found that the overall cellular composition stayed relatively stable after injury, with few exceptions. Astrocytes appeared to decrease in proportion at 1 and 3 wpi (p = 0.013 and p = 0.013) while microglia increased in proportion at 1 wpi (p = 0.006) and neurons increased in proportion at 3 wpi (p = 0.030).

To determine whether these observations accurately reflect endogenous changes in the lumbar spinal cord, we performed immunohistochemical staining and in situ hybridization experiments in tissue sections from healthy animals, 1 week, and 3 weeks after thoracic contusion injury. We stained for neurons (NeuN), astrocytes (GFAP and SOX9), oligodendrocytes (OLIG2), as well as microglia and macrophages (TMEM119, IBA1, Fig. 1d, e). Quantitative analysis of the multiplexed immunohistochemistry confirmed the snRNA-seq cellular composition. In the uninjured tissue, the proportions of astrocytes, microglia, and oligodendrocytes were not significantly different between the two technical approaches (astrocytes 11.7 ± 0.7%; microglia/hematopoietic 4.7 ± 0.3%; oligodendrocytes 30.5 ± 0.4%; mean ± SEM; p > 0.05; Supplementary Dataset 3). However, neurons were represented at a larger fraction in the snRNA-seq dataset compared to in tissue quantification of NeuN-positive cells (neurons 20.8 ± 0.6%, p = 0.040). This might reflect a decrease in non-neural cells at 3 wpi or a bias in our snRNA-seq dataset toward neurons, which had higher genes per nucleus.

While astrocytes appeared to decrease at 1 and 3 wpi in the snRNA-seq dataset, this result was not confirmed with immunohistochemistry in tissue (Fig. 1d, e). Rather, we observed no significant change in SOX9-expressing astrocytes over time (Fig. 1e). To further characterize the proportion of astrocytes using combinatorial RNA expression, we performed multiplexed RNA in situ hybridization using Agt, Gja1, and Aqp4 markers. In this context, a modest decrease was observed from 11.7% in uninjured lumbar cords to 9.5% 1 wpi (±0.7%, 0.9%, respectively, Supplementary Fig. 3). This suggests that the proportional change in astrocytes that we observed in the single-cell atlas did not reflect endogenous cell type changes, but is likely due to the overall proportional shifts of cell types in the injured spinal cord.

We emphasize the importance of in situ validation for cell proportion changes in single-cell RNA sequencing data to distinguish authentic differences in endogenous cell compositions. While some cell types were significantly different between snRNA-seq and in tissue proportions, the composition of cell types in tissue generally reflected those observed in the sequencing analysis. Overall, single nucleus RNA sequencing provides an unbiased profiling of cell types that reflects in-tissue spinal cord biology.

Composition and changes within 39 refined cell types

Given the size of this dataset at 67,903 nuclei, we were able to cluster at a higher resolution to identify rare cell types. In this second level of hierarchical clustering, we subclustered neurons, astrocytes, microglia/hematopoietic, and oligodendrocyte lineage/Schwann cells yielding 39 cell types (Fig. 2a–d). Replicates are shown in the uninjured spinal cord by cell type, calculated by the percent within a given sample (Fig. 2e, Supplementary Dataset 3). To assess if the subpopulations increase or decrease after injury, we used scCODA, a Bayesian model for compositional single-cell data analysis (Fig. 2f)23. The scCODA framework models cell type counts while considering negative correlative bias via joint modeling of all measured cell type proportions. Here, we highlight findings from the four major cell classes (neurons, astrocytes, hematopoietic, and oligodendrocyte lineage/Schwann cells), including cell markers, composition, and changes after injury.

Cell type composition in the uninjured spinal cord and after injury.

a–d UMAPs depicting subclustering of major cell types: a neurons, b astrocytes, c microglia/hematopoietic cells, d oligodendrocyte lineage and Schwann cells. e A bar plot showing the 39 cell types in the atlas and their relative percent in each sample in the uninjured spinal cord. Individual replicates (N = 3) are shown as well as mean ± SEM. f The relative composition of the 39 cell types comparing injured samples (1 dpi, 1 wpi, 3 wpi, and 6 wpi) to uninjured, generated using scCODA showing the final parameter output from scCODA (confidence interval shown as 3–97% high-density index around the mean). Cell types with an inclusion probability > 0.85 were deemed significant. Significance depicted with a red asterisk. N = 3. For the log FC of individual replicates, see Supplementary Fig. 8. Source data are provided as a Source Data file.

Neurons

The 23,651 neurons were subclustered and annotated by a previously established atlas of the lumbar spinal cord neuronal subtypes19 using label transfer (Seurat, see “Methods”). This yielded 17 neuronal populations, including 8 excitatory, 8 inhibitory, and motoneurons (Fig. 2a, Supplementary Fig. 2, Supplementary Fig. 4). Excitatory neurons were marked by expression of Slc17a6, inhibitory neurons by Gad1, Gad2, and Slc32a1, and motoneurons by Slc5a7, Chat, and Prph. The excitatory neuron families as defined by Russ et al.19 are Cpne4, Maf, Reln, Rreb1, Sox5, Megf11, ME (mid-excitatory), and VE (ventral excitatory). Inhibitory neuron families, are Adamts5, Cdh3, Pdyn, Npy, Chat, MI (mid-inhibitory), VI (ventral inhibitory), and CSF-c (cerebral spinal fluid-contacting neurons). The proportions of these neuronal subtypes did not significantly change after injury (Fig. 2f, Supplementary Dataset 3).

Astrocytes

Astrocyte subtypes were identified by subclustering 4525 nuclei from coarse clustering (Fig. 1c, Supplementary Fig. 3). Astrocytes included five subtypes that putatively reflect two homeostatic populations, white matter astrocytes and reactive astrocytes (Fig. 2b, Supplementary Fig. 3). All astrocytes expressed Gfap, Agt, and Gja1. The two homeostatic populations, “astrocytes 1” and “astrocytes 2” expressed Gpc5 and Slc7a10, while putative white matter astrocytes did not. White matter astrocytes expressed higher levels of Gfap, as well as A2m, Cd44, C3, and Vim. Reactive astrocytes significantly increase acutely after injury (from 0.3 ± 0.01% of a sample in the uninjured cord to 2.3 ± 1.2% of a sample 1 dpi, N = 3, Fig. 2f), and expressed Lcn2, Rgs20, Hpgd, Serpina3n, Iigp1, Gbp2, and Slc10a624. Both coarse and refined clustering of astrocytes indicated a decrease in astrocyte proportions (Figs. 1e and 2f); however, these changes were not confirmed with immunohistochemical detection in tissue as discussed above (Fig. 1d, e, Supplementary Fig. 3b). The apparent decrease in astrocytes after injury in the snRNA-seq dataset may be due to overall proportional increases in cell types in the injured spinal cord or selective vulnerability of astrocytes after injury.

Microglia and hematopoietic cells

To explore microglia and related hematopoietic cell types in greater depth, we independently clustered the 5080 nuclei and observed three homeostatic and two activated microglia populations, a cluster of macrophages, and a cluster of natural killer and T cells (Fig. 2c, Supplementary Fig. 6). Homeostatic microglial clusters were defined by Cst3, C1qa, Ctss, Hexb, Trem2, P2ry12, and Tmem119 (Supplementary Fig. 2). The activated microglia populations expressed lower levels of P2ry12 and Tmem119 and induced expression of the phagocytic markers Cd68 and Lyz2. In addition, “activated microglia A” expressed Gpnmb, Apoe, Lgals3, Igf1, and Spp1, while “activated microglia B” expressed genes associated with pro-inflammatory microglia, such as Ccl2, Ccl3, Ccl4, and Lpl. Macrophages expressed the genes Mrc1, Cd74, and H2-Ab1 and did not express the microglia-specific genes Tmem119 or P2ry12. Natural killer and T cells clustered together and expressed the genes Ms4a4b, Cd52, Ptprc, Nkg7, and Cd3g. All microglia/hematopoietic cell types increased in proportion relative to other cell types at 1 wpi, particularly activated microglia A, which increased 14.8-fold, from 0.1% (±0.04) of all cells in the uninjured cord to 1.6% (±0.60) 1 wpi (Mean ± SEM; Supplementary Fig. 6) as was reflected in the coarse cell type analysis. Both activated microglia subtypes were still present at 6 wpi, suggesting that they may play an ongoing role at chronic time points (Supplementary Fig. 6), even in lumbar tissue distal to the injury.

Notably, the gene expression profile of “activated microglia A” strongly resembled a signature observed recently in postnatal myelin-phagocytosing microglia, in postnatal microglia that can promote SCI repair, and in disease-associated microglia in conditions such as Alzheimer’s disease in the brain and ALS in the spinal cord25–28. In addition, recent work examining the lesion site of SCI has identified an “injury associated microglia” cell type with a similar expression profile29,30. Pathway analysis of the genes that characterized “activated microglia A” revealed an enrichment of genes associated with (1) phagocytosis (such as Lyz2, Gpnmb, Itgax, and Cd68; GO terms: lysosome, p = 3 × 10−13, antigen presentation, p = 0.0003, phagosome, p = 0.0363), (2) lipid metabolism (such as Fabp5, Lgals3, Apoe, Soat1, and Abca1; GO term: lipoprotein, p = 0.0041), and (3) secreted proteins (such as Spp1 (OPN) and Igf1; GO term: extracellular secretion, p = 0.0282, Supplementary Fig. 6c, Supplementary Dataset 5).

To determine whether the gene expression pattern that we observed corresponded to an in vivo cell type, we performed in situ hybridization in tissue sections of spared lumbar cords, using C1qa (a general microglial marker), Spp1 (OPN) (which marked activated microglia as well as some ventral horn neurons19,31,32), and Gpnmb (which was specific to activated microglia A). These cells were observed within the white matter of the injured spinal cord, and appeared consistently in small clusters along the putative rubrospinal tract (RST) and the dorsolateral corticospinal tract (CST) in the lateral white matter at 1, 3, and 6 wpi (Supplementary Fig. 6e–h). Interestingly, this region showed loss of longitudinal axons from the descending tracts but also showed no change in the presence of residual myelin (Supplementary Fig. 6i–l). Together, these data suggest that Activated microglia A cells are similar to previously described “DAM/PAM” cells and were found in the white matter of the injured spinal cord, where they may play a role in the phagocytosis of degenerating axons or apoptotic cells33.

Oligodendrocytes and oligodendrocyte lineage cells

Oligodendrocytes and oligodendrocyte lineage cells comprised the largest proportion of the lumbar spinal cord. We subclustered 32,287 oligodendrocyte and oligodendrocyte lineage cells, as well as 315 Schwann cells. Although the Schwann cells in this study are likely from lumbar spinal cord roots that were in the dissected tissue, we included them in the downstream analysis due to the interest in Schwann cells as a source of repair and remyelination after injury34. Oligodendrocyte lineage cells, including oligodendrocyte precursor cells (OPCs), expressed Cspg5 and Tnr; committed oligodendrocyte progenitor cells (COPs) expressed Fyn and Tcf7l225,26 (Fig. 2d). Newly formed oligodendrocytes (NFOL) expressed Man1a and Synpr, whereas myelinating and mature oligodendrocytes expressed Plp1, Mag and Mog (Supplementary Figs. 2 and 7a–c). Myelinating oligodendrocytes expressed Opalin and Kirrel3, representing a transitional population between newly formed and mature oligodendrocytes. Mature oligodendrocytes expressed Klk6 and Art3. Schwann cells expressed Mpz and Pmp22. The proportion of oligodendrocyte lineage cells did not significantly change after injury with the exception of COPs, which increased fivefold from 0.3% (±0.1) of cells in the uninjured spinal cord to 1.7% (±0.2) of cells at 1 wpi (Fig. 2f, Supplementary Dataset 3). These committed oligodendrocyte precursor cells resemble previously identified “COPs” in the juvenile mouse25,26. The expansion of this population after injury indicates the emergence of new myelinating cells27,28 or potential roles for COPs themselves29.

To determine if the oligodendrocytes identified here were similar to those previously described after spinal cord injury, we compared our oligodendrocyte lineage cells to those identified in Floriddia et al.30. In this previous study, the diversity of mature oligodendrocytes was cataloged, including spatial distribution and susceptibility to spinal cord injury. Of interest, this study found MOL2 and MOL5/6 to be enriched in the spinal cord, with MOL5/6 increasing with age and enriched in the spinal cord lesion site. We analyzed how oligodendrocyte subtypes compare in both studies in the different temporal conditions (using Pearson’s correlation of the shared top 2000 variable genes). Oligodendrocytes from both studies correlated by cell type rather than injury (Supplementary Fig. 7). MOL2 from Floriddia et al. correlated most with the mature oligodendrocytes 1 and 2 (MOL-1, MOL2) from this study across injury conditions. MOL2 was enriched in distal areas of the injury site, where Wallerian degeneration took place. Both datasets similarly classified OPCs, COPs, and MFOLs. These similarities support our classification of oligodendrocyte subtypes.

In all, we provide an in-depth analysis of 39 cell types in the lumbar spinal cord and their compositional changes after injury. The relative abundance of most cell types did not significantly change, except for several glial subtypes (Fig. 2f). The increase in reactivity of cells such as reactive astrocytes at 1 dpi, and microglia at 1 and 3 wpi likely reflected a coordinated response by glia to the assault on the lesion site above. By 6 wpi, only the proportion of were astrocytes was significantly different, indicating a relative return to a native state by this chronic time point.

Cell type-specific changes in gene expression in the lumbar spinal cord after injury

While the overall composition of cell types is relatively stable after injury, we hypothesized that gene expression within cell types would change after injury, endowing particular cell types with new properties and functions. We next analyzed gene expression across cell types to understand how specific cell types change their molecular repertoire after injury. It is important to note that statistical differences in gene expression can largely be driven by the size of the cluster, with larger clusters having the power to resolve more differentially expressed genes. To determine which cell types changed significantly following injury, we implemented Augur, a method to rank responsiveness of cell types in single-cell data that is not biased by cluster size9. We applied Augur to the 39 cell types detected in this study, including 17 neuronal clusters, 4 astrocyte clusters, 7 microglial/hematopoietic clusters, 6 oligodendrocyte-related clusters, Schwann cells, endothelial cells, pericytes, ependymal cells, and leptomeninges. A cell type prioritization score was generated from Augur, indicating the responsiveness of cell types after injury. We found that the average cell type prioritization score across all cell types decreased with time after the acute response to injury, suggesting a progressive return to homeostatic conditions across multiple populations (average 0.55 AUC, Supplementary Fig. 9a, Supplementary Dataset 4). At 1 dpi we found that microglia were the most perturbation-responsive cell type (0.86 average AUC for Microglia 1, 2, and 3). These results suggest that microglia could play an important role in the immediate phase after injury, even in spared tissue distant from the site of injury.

To understand the cellular processes driving these changes, we examined differentially expressed genes within these cell types at each stage after injury and performed gene ontology (GO) and pathway analysis (Supplementary Fig. 9, 10, Supplementary Dataset 5). At 1 dpi, microglia/hematopoietic cells changed dramatically and were characterized by a burst of expression of genes related to cell metabolism, which may support the significant expansion in microglial cell numbers that we observed above. Reactive astrocytes also emerged at 1 dpi, expressing pan-reactive genes such as Lcn2 and Serpina3n, as well as markers for pro-inflammatory and neuroprotective astrocytes24 such as Gbp2 and Slc10a6 (Supplementary Fig. 9b). Similarly, endothelial cells and meninges displayed acute gene expression changes at 1 dpi, particularly in molecules related to structure, cell–cell connections, and extracellular matrix (ECM) adhesion. In addition, nearly all cell types showed increased expression of Mt1 (metallothionein-1) and Fth1 (Ferritin Heavy Chain 1) genes which are both involved in binding heavy metals, including iron. This suggests that extravascular blood may be an early environmental cue to reach the tissue of the lumbar cord (Supplementary Fig. 9b).

While non-neuronal populations changed most acutely, neurons displayed more delayed responses. One week post injury, neuronal populations showed altered expression of genes related to cell stress, including oxidation-reduction processes and protein folding stress response, and molecules related to neurotransmission and ion channel activity. Intriguingly, specific populations of excitatory neurons in the dorsal horn and inhibitory neurons within the ventral horn displayed changes in synaptic organization and plasticity-relate genes. This suggests the potential for tailored circuit remodeling. At the same time, oligodendrocytes altered their cellular metabolism molecules and genes related to cell–cell adhesion. Three weeks after injury, multiple neuronal populations continued to show signatures of cell stress and changes in cell–cell contacts, while oligodendrocytes continued to show changes related to cell metabolism and cell–cell adhesion. Finally, by 6 wpi, many cell types showed modest or no changes in gene expression compared to the uninjured state, showing an overall return to baseline expression patterns by this point (Supplementary Fig. 9b).

Together with the cellular composition analysis above, these data comprise a natural history time course of the lumbar spinal cord response to injury. Within each time window or epoch, the community of cell types responds with diverse molecular signatures which can be observed across time.

Transcriptional changes in neurons of the lumbar spinal cord after injury

Neurons in the lumbar cord are largely spared, unlike the neurons at the lesion site35. However, many of these cells undergo axotomy and abrupt changes to their descending input36,37. We sought to characterize the changes within all neurons after injury. Histological analysis of the tissue, revealed no significant change in the number of neurons (Fig. 1e), neurofilament signal, or MAP2+ dendritic signal (Fig. 3b, c) after injury, suggesting no large-scale death or loss of neuronal processes. However, overall, neurons displayed dynamic changes in gene expression, particularly at 1 and 3 wpi, with more genes enriched in uninjured neurons compared to injured neurons (Fig. 3d). Pathway analysis showed that mitochondrial, endoplasmic reticulum, and ATP synthesis pathways are enriched in the uninjured cord (Fig. 3e), likely reflecting a greater level of neuronal firing and homeostatic function. At 1 and 3 wpi, pathways associated with plasticity such as post-synaptic density, neurotransmitter receptors, and axon guidance were upregulated. Particular genes that are upregulated include neurotransmitter receptors (such as Gria1, Gria2, Gria3, Grid1, Grik1, and Chrm2), those related to synaptogenesis (Nrnx3, Nlgn1, Lrrmt4, Tenm2, Lrrc4c, and Dscam), and synaptic structure (Bdnf, Stat3, Socs3, Klf6, Gap43, Atf and Sprr1a, Fig. 3g). The upregulation of these pathways suggests a broad change in neurotransmission and synaptic remodeling amongst neurons in the first few weeks after spinal cord injury. These changes would likely render spinal neurons more responsive to low levels of neurotransmitters and their activity within local circuits.

Plasticity-related expression in neurons after injury.

a Schematic depicting lumbar spinal cord neurons and their response to injury, whether that be an ascending neuron or an interneuron. b Immunohistochemistry of the lumbar spinal cord from uninjured, 1 wpi and 3 wpi. Tissue was stained for neurofilament (a cocktail of neurofilament-light, neurofilament-medium, neurofilament-heavy; purple) and MAP2 (green). Scale bars are 200 µm. c Quantification of neurofilament and dendritic changes. Pixels were quantified from thresholded images of neurofilament and MAP2. Error bars are mean ± SEM (N = 4). d Differential gene expression analysis comparing uninjured to 1 and 3 wpi neurons, the time points of maximal neuronal gene expression changes. Black dots indicate p value adjusted < 0.001, gray indicate p ≥ 0.001. e Pathway analysis for differentially expressed genes between uninjured and injured time points. _X_-axis indicates −log(p val adj) of GO and KEGG pathway clusters. P values (adjusted) were calculated using Benjamini–Hochberg false discovery rate (FDR). Yellow indicates relatively high normalized average expression and dark blue indicates relatively low normalized average expression. f Chord plot indicating shared genes between top 5 GO terms from genes upregulated 1 and 3 wpi. g Heatmaps showing average neuronal gene expression from top GO terms, including neurotransmitter receptors, synaptogenesis, and synaptic structure. Yellow indicates relatively high normalized average expression (1) and dark blue indicates relatively low normalized average expression (−1).

Rare populations of spinal neurons induce a gene expression signature of regeneration

We observed broad changes in the expression of genes related to neural excitability, plasticity, and circuit structure that could support tissue-wide changes in function (Fig. 3). In addition to these broad effects that could alter local spinal network function, we hypothesized that specific neuronal populations might be capable of spontaneous long-range remodeling. Such changes would be challenging to observe without the resolution of single nucleus sequencing, and thus the dataset that we generated held opportunities for discovery.

As was described above, we determined the refined identities of neuronal populations in the dataset using annotations described in a recent atlas of mouse spinal cord cell types19. While we observed nearly all the cell types that we expected, one group of neurons remained challenging to categorize (Fig. 4a, b, Supplementary Fig. 5b, c). By clustering neuronal populations without label-transfer annotations, we found that these neurons were distinguished by genes associated with axon regeneration, including Atf3, Sprr1a, Tnfrsf12a (Fn14), Sox11, Klf6, Bdnf, Adcyap1, and elevated expression of Gap43 (Fig. 4c, Supplementary Fig. 4c, d). These genes all belong to a class of “regeneration-associated genes” (RAGs) that is induced robustly in the peripheral nervous system after nerve injury and can enable regeneration of axons14,15,38. The expression of individual regeneration-associated genes has been reported in the spinal cord after injury38–43. However, the full spatiotemporal profile of this rare central nervous system phenotype, the cellular identity of these neurons, and the association with axon outgrowth in this context are all unknown. This has been a major obstacle in understanding what enables, restricts, or modulates circuit reorganization after injury.

Specific neurons express genes associated with regeneration.

a UMAP showing predicted neuronal families. b Targeted view of RAG-expressing cluster over injury time points. c Featureplots showing RAGs expressed in neurons. d A dotplot showing expression of RAGs within targeted RAG-expressing cluster (cluster 23 defined by independent clustering in Supplementary Fig. 5). Average expression (avg exp) is indicated by color (gray to blue) and percent expressed (pct exp) is indicated by the size of the circle. e, f RNAscope in situ hybridization showing expression of Sprr1a and Atf3 in the uninjured cord and 1 wpi. Scale bars are 200 and 50 µm, respectively. g Quantification of _Sprr1a_+ cells and _Atf3_+ cells in the uninjured and 1 wpi injury cord (p = 0.001 shown as ***p = 0.0037 shown as **, two-sided unpaired t-test, Error bars indicate ± SEM, N = 7, 10). h Diagram of spatial location of transcription types, including ChAT (light green), Vsx2 (dark green), and Shox2 (orange). i–l Visualization and quantification of VGluT2/Slc17a6+, Shox2+, Chx10/_Vsx_2+, and Sprr1a_-expressing cells in the ventral spinal cord. Scale bars are 50 µm. Error bars indicate ± SEM (N = 5, 4, 4, 6 animals). m Diagram of spatial location of connectivity types, including ascending neurons labeled by dextran (aqua) and spinocerebellar (SCT, orange) neurons. n Visualization and quantification of lumbar spinal cord neurons labeled by dextran injected into a thoracic contusion lesion site. Aqua arrows indicate ATF3 and dextran overlap. Scale bars are 50 µm. Error bars indicate ± SEM (N = 5 animals). Visualization and quantification of the percent of RAG-expressing cells—_Sox11 (o), Sprr1a (p), and ATF3 (q) that are spinocerebellar. Spinocerebellar neurons are shown in green and RAGs are shown in red. Orange arrows indicate RAG gene and spinocerebellar dual-labeled cells. Scale bars are 50 µm. Error bars indicate ± SEM (N = 4, 4, 5 animals). Source data are provided as a Source Data file.

We examined the expression of cell type marker genes in this RAG+ cluster and found evidence of a mixed cell type origin. Cells with the highest RAG expression largely downregulated their endogenous gene expression, as previously reported in peripheral neurons, confounding the initial molecular definition of their cell type14,15. Despite this caveat, the RAG+ cluster expressed a diverse set of genes associated with ascending projection neurons from the lumbar spinal cord to the brain, including Lypd1, Tacr1, Zfhx3, Pou6f2, and Tac144–47 (Fig. 4d). Additionally, the expression of Slc17a6 (vGlut2), Zfhx3, and Pou6f2 suggested that some of these cells were likely excitatory neurons that resided in the lateral part of deep dorsal or ventral horn19.

We next probed the expression of Sprr1a and Atf3 in tissue to test whether the RAG signature in the sequencing data is reflected in vivo. Expression of Sprr1a and Atf3 were observed in lumbar spinal cord tissue beginning at 1 week following spinal cord injury and extending to chronic time points at three and 6 weeks after injury (Fig. 4e–g, Supplementary Fig. 5g, Supplementary Fig. 15), thereby confirming the transcriptional data above. Spatial analysis using in situ hybridization revealed a rostral-caudal gradient in the number of RAG-expressing cells, with a greater number of _Sprr1a_-positive cells at rostral lumbar segments closer to the lesion site (Supplementary Fig. 11a, b). The cellular distribution within the transverse plane also shifted such that RAG-expressing cells were found in the dorsal, mid, and ventral horns at L2, but were restricted to the ventral horn in L5. We compared RAG expression to that of the excitatory marker Slc17a6 and confirmed that most RAG-expressing neurons were indeed excitatory (Fig. 4i), and we next examined whether these represent any of the V0c, V2a, or V2d ventral excitatory populations48,49. We did not detect any co-expression of RAGs with the V0c marker Chat or the V2a marker Vsx2. In contrast, a small subset of RAG neurons expressed Shox2 (Fig. 4j, k), a marker of V2d excitatory, rhythm-generating central pattern generator neurons of the ventral horn50–52.

To determine if certain neurons have a molecular predisposition to express RAGs over others, we applied single nucleus ATAC-seq to the lumbar spinal cord of uninjured mice. We leveraged our snRNA-seq annotations to identify cell types, first assessing neurons compared to non-neurons (Supplementary Fig. 12a–d). We detected increased chromatin accessibility of Sprr1a, Sox11, and Gap43 RAG loci within neurons. To provide greater resolution of neurons, we next grouped these into families of cell types, including dorsal excitatory (DE), dorsal inhibitory (DI), mid-excitatory (ME), mid-inhibitory (MI), ventral excitatory (VE), ventral inhibitory (VI) and motoneurons. We found no clear differential pattern of chromatin accessibility between families of neurons (Supplementary Fig. 12e–h). Although this dataset did not provide deep neuronal subtype resolution, these findings do not support the existence of specific neuronal subpopulations that are primed to express RAG genes.

If molecularly defined populations do not explain the majority of RAG-expressing neurons, what other factors should be considered? Ascending projection neurons are important candidates for RAG expression after injury based on the expression of ascending markers in the RAG+ cluster and the passage of these neurons through the lesion area. We next hypothesized that the lumbar projection neurons with axons directly injured by the thoracic contusion respond by expressing this pro-regenerative signature.

To determine whether direct injury to axons elicited RAG expression, we performed complete transection injuries of the thoracic spinal cord plus dextran injection to the transected area (Supplementary Fig. 13a). One week later, we found that 4.6% of neurons in the lumbar cord were labeled by dextran (taken up by cut axons and transported to the cell body), with less than 1% of total neurons expressing ATF3 (Supplementary Fig. 13b, c). Of the ATF3-expressing neurons after complete transection, less than half (48%) were dextran positive (Supplementary Fig. 13b, c). We repeated this study using a lateral hemisection injury to distinguish ipsilateral and contralateral-projecting populations and found similar results (Supplementary Fig. 13d, e). Thus, ATF3 is either expressed by both injured and non-injured neurons or the dextran labeling was not complete in labeling all directly injured neurons. More importantly, we found that only a subset of lumbar projection neurons that were directly injured by transection (and took up dextran label at the injury site) induced expression of ATF3. Finally, we examined what percent of lumbar neurons are axotomized during a contusion, by administering dextran to the contusion site at the time of injury. In a contusion model, we found 3 weeks post injury that 18% of neurons were directly injured ascending neurons, as marked by dextran (18%, Supplementary Fig. 11e). In conclusion, the expression of RAGs after the injury is not found in a single molecularly defined population, nor is it a general feature of axotomized spinal neurons.

Spinocerebellar neurons express RAGs and display axon sprouting below the lesion

After the injury, the rare RAG-expressing neurons in the lumbar spinal cord displayed a heterogeneous gene expression pattern. Still, their rostral-caudal distribution, location in the deep dorsal and ventral horn, and excitatory neurotransmitter status suggested that they may be spinocerebellar neurons. Previously, Shox2 has been reported as a marker of spinocerebellar neurons in the cervical cord53. We found that this is not the case in the lumbar spinal cord (0% of spinocerebellar neurons were Shox2 +, ±0.0, N = 3 animals). To test if these neurons are spinocerebellar, we injected AAV2retro-pmSyn1-EBFP-Cre in the cerebellum and AAV1-Syn-DIO-GFP bilaterally in the lumbar spinal cord to label both ipsi- and contra-laterally projecting spinocerebellar neurons (Fig. 5a). Two weeks later, we delivered a severe thoracic (vertebral level T9) contusion and dextran to label directly injured neurons. After bilateral spinocerebellar labeling and 3 weeks post thoracic injury, we found that 65.2% (±4.3) of _Sox11-_expressing neurons, 41.8% (±1.5) of _Sprr1a_-expressing neurons and 38.0% (±5.6 SEM) of ATF3-expressing neurons were indeed spinocerebellar (Fig. 4o–q, Supplementary Fig 11d, e), while only 1.3% of the spinocerebellar neurons were dextran positive suggesting that spinocerebellar neurons are largely not axotomized by the severe thoracic contusion injury (Supplementary Fig. 11e). Ventral spinocerebellar neurons may avoid direct injury based on the ventral location of their axons. From this data, at least ~55% of ATF3-expressing neurons are ascending, as revealed by dextran labeling (18% directly injured) or spinocerebellar viral labeling (38%). Based on the inefficiencies and variability of these labeling approaches, this likely under-represents the proportion of ascending neurons amongst the RAG+ population.

Spinocerebellar neurons express RAGs and display thoracic sprouting after injury.

a Schematic of AAV- injection and injury. b–d Spinocerebellar neurons and their cell bodies, axons, and mossy fibers in the cerebellum, thoracic and lumbar spinal cord. Virus expression is shown in green. Scale bars are 500 µm in the cerebellum and 200 µm for spinal cord sections (thoracic and lumbar). e Quantification of dendritic arborizations in SCT neurons (ns, p = 0.206, p = 0.211, Mann–Whitney test). f Quantification of thoracic axons rostral (R) and caudal (C) to the injury site (p = 0.0012, p = 0.035, Mann–Whitney test). g Quantification of gray matter collaterals rostral (R) and caudal (C) to the injury site, as measured by pixels after thresholding (p = 0.0485, two-way ANOVA). h Quantification of mossy fibers, normalized by the number of SCT neurons in the same animal (p = 0.336, p = 0.039, two-sided unpaired _t_-test). Mean ± SEM, N = 4, 5, 6 animals. *p < 0.05; ****p < 0.0001; n.s. not significant. Source data are provided as a Source Data file.

We next asked whether there is there an axon outgrowth or remodeling phenotype that correlates with the RAG gene expression signature that we observed in spinocerebellar neurons? Based on the cell-filling viral label, we examined dendritic arborizations, thoracic axons, collateral axons, and mossy fiber terminals of spinocerebellar neurons after injury and performed this analysis in two independent experimental paradigms: bilateral or contralateral spinocerebellar labeling. Dendritic structure did not change after injury (Fig. 5e, p = 0.206, p = 0.211, Mann–Whitney test). In contrast, we observed an increase in the number of axons found in the white matter of the thoracic spinal cord below the injury site at both three and 6 weeks after injury (Fig. 5f, p < 0.05, p < 0.05, Mann–Whitney test) and a dramatic increase in gray matter collaterals of these axons caudal to the injury site (p < 0.005, two-way ANOVA). Specifically, these collaterals targeted lamina VII of the ventral horn (Fig. 5b–d). In contrast to these findings, there was no significant change in spinal cord axons in the tissue above the lesion (Fig. 5f, ns). In the cerebellum, there was a trend for an increase in spinocerebellar mossy fibers at 3 weeks after injury (p = 0.336, Mann–Whitney test), which then decreased significantly by 6 weeks after injury (Fig. 5g, p < 0.05, Mann–Whitney test). This indicates that the terminal synapses of spinocerebellar neurons are largely preserved after spinal cord injury and may display dynamic reorganization. Together, these results show that lumbar spinocerebellar neurons expressed regeneration-associated genes after spinal cord injury, have axons that were spared after a severe contusion, and underwent structural remodeling below the injury site including axonal outgrowth after injury. This highlights an example of spontaneous neuronal remodeling after spinal cord injury, discovered through the power and resolution of single nucleus sequencing.

Discussion

SCI disrupts neuronal connections, eliciting trauma on a wide array of cells within the tissue. While the capacity for axonal regeneration and recovery after SCI is limited54,55, there may be latent mechanisms for spontaneous recovery within the spinal cord, particularly in anatomically incomplete injury where reorganization of circuits below the lesion site can support the restoration of function54,55. We used single nucleus RNA sequencing to profile the lumbar spinal cord after a severe thoracic contusion injury to track the cell type-specific injury responses and identify spontaneous changes that could be leveraged for recovery. We present a natural history time course extended from acute through chronic time points and an accompanying interactive website as a resource for the field (https://seqseek.ninds.nih.gov/spinalcordinjury). We observed rare spinal neurons that expressed a pro-regenerative transcriptional signature, identified a major subset of these cells as spinocerebellar neurons and demonstrated axonal sparing and outgrowth of these cells after spinal cord injury.

Our findings build on prior work that used single-cell or nucleus sequencing to profile the cell type-specific responses to spinal cord injury. Most of these studies have focused on non-neuronal cells in and surrounding the lesion area or throughout the spinal tissue, such as the cells that make up the fibrotic core, glial scar, myeloid cells, vascular cells, and oligodendrocyte lineage30,56,57. In particular, an emerging finding from multiple studies, bolstered by our own observations, is the presence of damage-associated microglia in disease and injury56,58–64. Of the microglia that respond after SCI, activated microglia A cells expressed neuroprotective growth factors and were transcriptionally similar to proliferative axon tract-associated microglia in postnatal mice33,60,62, disease-associated microglia58, as well as microglia from the human spinal cord65. Future work is needed to address the implications of this disease-associated microglia for recovery after injury in adults61. While these studies have highlighted the spontaneous trauma responses amongst glia, they have left neuronal plasticity mechanisms largely unaddressed. We previously used single nucleus RNA sequencing to identify neurons that respond to epidural electrical stimulation66, but did not examine the effect of spinal cord injury itself or explore the molecular pathways that could mediate spontaneous or therapeutic circuit remodeling.

The cell type- and context-specific factors that control a neuron’s response to injury are not well understood, hampering efforts to target, expand, and modulate this cellular potential. We leveraged the resolution of single nucleus sequencing to explore this question in the context of spinal cord injury. We found that broad neuronal responses mainly included the downregulation of physiological pathways and the upregulation of genes associated with neurotransmission and synaptic structure. In contrast, two distinct subsets of spinal cord neurons—Shox2-expressing V2d neurons and spinocerebellar neurons—expressed RAGs after injury. These genes included transcription factor RAGs (such as Atf3 and Sox11) that act as indirect regulators of many growth proteins, as well as effector RAGs (such as Sprr1a and Gap43) that play a direct effector role in growth cone and axon outgrowth cytoskeletal changes40,67,68. By defining the identities of spinal neurons that express RAGs, our work opens the door to tracking the axonal remodeling that may accompany the transcriptional regeneration signature.

RAGs can be necessary and sufficient to support axon regrowth after injury14,15,39,40,69,70, prompting us to test whether spinocerebellar neurons displayed structural plasticity after injury. Indeed, we found that spinocerebellar axons are spared by severe contusion injury and even show increased numbers in the thoracic spinal cord, together with enhanced axon collaterals in the thoracic gray matter and evidence of cerebellar mossy fiber reorganization. Importantly, these structural changes confirm that spinocerebellar neurons, that express RAGs after spinal cord injury, undergo spontaneous axon outgrowth.

What may be the functional consequences of spinocerebellar circuit remodeling? Sprouting of existing axon fibers is an important component of recovery from spinal cord injury, permitting spared neurons to make new connections, serve as circuit relays, and take on compensatory roles for improved behavioral function13,71–73. We found that spinocerebellar neurons showed significant remodeling in the caudal thoracic cord just below the lesion, with a major expansion in the lamina VII of the ventral horn. Here, these neurons may contact local central pattern generator circuits or connect with spared and regenerating descending pathways37,74. Spinocerebellar neurons also showed structural plasticity at their mossy fiber terminals in the cerebellum. Given the importance of these connections in locomotion and motor learning, such anatomical changes could provide a key substrate for therapeutic interventions75,76.

There are several limitations to consider in this work. First, we extracted nuclei instead of whole cells for transcriptional profiling to avoid experimentally induced gene expression and selective cell death common in single-cell profiling18,32,77,78. However, this approach may yield fewer genes detected per cell/nucleus and may slightly bias the cellular composition79. For example, snRNA-seq on the human brain showed that using nuclei for transcriptional profiling depleted activated microglia, compared to using cells80. While it is possible that we under-represented activated microglia in our data, both single-cell and single nucleus RNA sequencing approaches have limitations81 that are important to consider when choosing a technique. Second, this study provided a global overview of all cell types in the lumbar spinal cord, and should be followed up by future studies on specific cell types after spinal cord injury enabling deeper analysis of the molecular pathways and trauma responses of each cell type. Third, the atlas component of this work examines changes at the gene expression level and does not address post-transcriptional cellular mechanisms. Despite these limitations, this work provides a powerful and temporally resolved reference atlas of cell type-specific changes after traumatic injury and allowed us to discover rare molecular changes that could provide the substrates of repair and recovery.

Here, we sought to uncover the endogenous mechanisms by which the adult nervous system can recover from a severe SCI. The single-cell atlas, discovery of RAG-expressing neurons, and the plasticity in spinocerebellar neurons following severe but incomplete thoracic contusion injury provide important insight into the natural mechanisms of recovery after SCI. Understanding of these intrinsic mechanisms will provide therapeutic targets to control or even reverse pathological changes across a wide variety of injuries and diseases.

Methods

Mice

This study including all procedures and experiments received ethical approval from both the Veterinary Office of the Canton of Geneva (Switzerland) and the National Institute of Neurological Disorders and Stroke Animal Care and Use Committee (ACUC protocol number 1384). Mice for RNA sequencing were female C57BL/6 (12–30 weeks of age). For all other experiments, balanced samples of male and female C57BL/6 mice (12–30 weeks of age) were used. Mice used in this study were housed under a 12-h light–dark cycle (06:00–18:00 light), with ad libitum access to food and water. Room temperature was between 20–26 °C and humidity was between 30–70%.

Surgical procedures

Surgical procedures were performed as previously described82. Briefly, following a mid-thoracic laminectomy (T9 vertebra), a spinal cord impactor (IH-0400 Impactor, Precision Systems and Instrumentation LLC) was used to induce a contusion injury. The applied force was set to 90 kdyn. Spinal transections were performed following a mid-thoracic laminectomy (T9 vertebra), cutting the spinal cord with spring scissors before filling the void with gel foam. Animal care, including manual bladder voiding, was performed twice daily or as needed following injury.

For dextran-labeling experiments, 1 μL dextran (Rhodamine B, 10,000 MW, ThermoFisher Scientific, Catalog Number D1824) was injected into gel foam separating the transected cord.

Two days after injury, all mice were evaluated in an open field, and all animals exhibiting any hindlimb movements were not further studied. For single nucleus RNA sequencing experiments, a larger cohort of mice was taken through kinematic analysis, and three mice representative of each time point were selected. For histology experiments, at least four mice were used for each condition. Mice with bone-hit contusions or injuries that fell a standard deviation outside of the average behavior for each time point were excluded. Due to animal care requirements injury experiments were not performed blinded.

Viruses

AAV2retro-pmSyn1-EBFP-Cre virus was produced at Addgene (Plasmid #51507)83. AAV1-Syn-DIO-GFP virus was produced by Vigene (CV17077-AV1). Viral particles were injected at a titer of 5E12−1E13 genome copies per mL.

Intraspinal injections

Intraspinal injections were performed as previously described84,85. Briefly, mice were anesthetized by intraperitoneal injection of a drug cocktail containing fentanyl (0.2 mg kg), dexmedetomidine (1 mg/kg), and midazolam (5 mg/kg) dissolved in saline. For spinal injections, a small incision was made in the skin, and the underlying musculature and adipose tissue were teased apart to reveal the vertebral column. Tissue joining the dorsal processes of consecutive vertebrae was removed, and the vertebral surfaces were cleaned with fine forceps and gently separated to reveal the dorsal surface of the spinal cord (at spinal levels L2 and L5). The dura was punctured by pinching with sharp forceps to facilitate the smooth entry of the needle. A glass pulled needle was lowered to a depth of 300 μm, halfway between the midline of the spinal cord and the lateral edge. The needle was then pulled back to 250 μm before pressure injecting 250 nL of viral particles at a rate of 75 nL per minute. Following virus injection, the needle was left in place for one minute before it was removed. Each spinal cord received three unilateral injections at L2, L3/4, and L586.

The overlying muscle was sutured after injections, and the skin incision was closed using wound clips. Anesthesia was reversed by intraperitoneal administration of buprenorphine (0.1 mg/kg), atipamezole (2.5 mg/kg), and flumenazil (0.2 mg/kg) in saline. Additionally, mice received an intradermal injection of meloxicam SR for analgesia and were returned to their home cages.

Mice were excluded when viral labeling showed less than 20 cells in the lumbar spinal cord.

Intracranial injections

Mice were anesthetized by intraperitoneal injection of a drug cocktail containing fentanyl (0.2 mg kg), dexmedetomidine (1 mg/kg), and midazolam (5 mg/kg) dissolved in saline. A small incision was made in the scalp, and a craniotomy was made at −5.8 mm AP, −4.0 mm ML, referencing from bregma (1Cb-4/5Cb)87. Virus (500 nL volume) was pressure injected through a pulled glass needle at a rate of 150 nL per minute, starting at a depth of 1.8 mm, slowly raising to a depth of 1.5 mm. Following virus injection, the needle was left in place for one minute before it was removed. The craniotomy was closed with gel foam followed by bone wax, and the scalp was closed with a wound clip. Anesthesia was reversed by intraperitoneal administration of buprenorphine (0.1 mg/kg), atipamezole (2.5 mg/kg), and flumenazil (0.2 mg/kg) in saline. Additionally, mice received an intradermal injection of meloxicam SR for analgesia and were returned to their home cages.

Behavioral assessments

All procedures used for mice in the sequencing experiment have been described in detail previously82. Limb movements were evaluated while running on a horizontal walkway. Bilateral leg kinematics were captured with the Vicon Motion Systems, UK (combining 12 infrared cameras) for tracking with reflective markers on the crest, hip, knee, ankle joints, and distal toes. The limbs were modeled as an interconnected chain of segments. Based on these, a total of 80 gait parameters were computed for each limb for each gait cycle. All gait parameters are reported in Supplementary Dataset 1. In Supplementary Fig. 1, quantifications of key gait parameters are shown: Step height (calculated from the toe), % drag (as percent of the gait cycle), whole limb oscillation (difference between maximum and minimum angle of limb axis, that is crest to toes, in XY plane within a gait cycle), whole limb oscillation velocity (velocity of the previous), ankle/knee/hip joint oscillation (difference between maximum and minimum joint angle within a gait cycle) and ankle joint oscillation velocity (velocity of previous). Differences among groups were calculated using two-tailed _t_-tests (unpaired) and were considered significant if p < 0.05. Data are represented as mean ± SEM unless otherwise indicated. Statistical analyses were performed using GraphPad prism software.

Analysis of kinematic data

A total of 78 gait parameters were computed for each limb for each gait cycle. We chose to represent the following gait parameters: step height, drag percentage, amp limb, amp speed limb, amp join 1, amp joint 2, amp joint 3, and amp speed joint 3. Differences among groups were calculated using two-tailed _t_-tests (unpaired) and were considered significant if p < 0.05. Data are represented as mean ± SEM unless otherwise indicated. Statistical analyses were performed using GraphPad Prism software.

Nuclei isolation

Nuclei were isolated from adult mouse lumbar cords with proximal dorsal and ventral roots attached using a triton-based protocol adapted from Matson et al.22. Briefly, mice were euthanized according to IACUC guidelines. The spinal cord was rapidly dissected and frozen on dry ice. Later, fresh frozen lumbar cords (spinal segment L2-S1) were placed in a Dounce homogenizer (Kontes Dounce Tissue Grinder) containing 500 μL of lysis buffer (0.32 M sucrose, 10 mM HEPES [pH 8.0], 5 mM CaCl2, 3 mM MgAc, 0.1 mM ETDA, 1 mM DTT, 0.1% Triton X-100). The cords were dounced with 5 strokes of pestle A, then 5-10 strokes of pestle B. The lysate was diluted in 3 mL of sucrose buffer (0.32 M sucrose, 10 mM HEPES [pH 8.0], 5 mM CaCl2, 3 mM MgAc, 0.1 mM ETDA, 1 mM DTT) and passed over a 40 μm strainer. The filtered lysate was centrifuged at 3200 × g for 10 min at 4 °C. After centrifugation, the pellet was resuspended in sucrose buffer and incubated for 2 min on ice. The sample was transferred to an Oak Ridge tube and homogenized for 1 min using an Ultra-Turrax Homogenizer (IKA). Then, 12.5 mL of density sucrose buffer (1 M sucrose, 10 mM HEPES [pH 8.0], 3 mM MgAc, 1 mM DTT) was layered below the sample. The tube was centrifuged at 3200 × g for 20 min and the supernatant immediately poured off. The nuclei on the side of the tube were resuspended with 100 μL of PBS with 0.04% BSA and 0.2 U/μL RNase inhibitor. Nuclei were inspected for visual appearance and cell lysis using trypan blue and quantified with a hemocytometer before being adjusted to a concentration of 1000 nuclei/μL.

Single nucleus RNA sequencing

Single nucleus RNA sequencing was carried out using single-cell gene expression 3′ v2 kit on the Chromium platform (10X Genomics) according to the manufacturer’s instructions with one modification. Following reverse-transcription, an additional PCR cycle was added to the number of cycles for cDNA amplification to compensate for decreased cDNA abundance in nuclei compared to cells. Approximately 8000–9000 nuclei were loaded in each well and 3000–7000 nuclei were recovered per sample.

Libraries were sequenced to a minimum depth of 20,000 reads per nucleus using an Illumina HiSeq 3000 (PE 26–8–98 bp). Raw sequencing reads were demultiplexed, aligned, and a count matrix was generated using CellRanger. For alignment, introns and exons were included in the reference genome (mm10).

Clustering

Analysis was performed in R (version 3.6.1), using the following packages: Seurat (version 3.2.2), RColorBrewer (version 1.1-2), ggplot2 (version 3.3.2), ggrepel (version 0.8.2), dbplyr (version 1.4.4), hrhrthemes (version 0.8.0), plyr (version 1.8.6), viridis (version 0.5.1), tibble (version 3.0.3) and tidyverse (version 1.1.0).

Seurat v3.2.2 was used to filter, normalize, anchor, and cluster the dataset88. We filtered nuclei for downstream analysis and included only those with greater than 200 genes detected per nucleus and less than 20% of reads coming from mitochondrial genes. For neurons, the minimum threshold was increased to 500 genes per nucleus as we have previously found that more genes are detected per nucleus in neurons compared with other cell types18,19. All genes analyzed were present in greater than three nuclei.

We performed integration of the five conditions (three samples each) using the uninjured samples as a reference, followed by standard log normalization and scaled the data using 2000 most highly variable genes, while regressing out percent mitochondria and nUMI. We used principal component analysis for dimensionality reduction. The number of principal components was selected based on elbow plot inflection, jackstraw plot significance, and PC heatmaps (inspecting gene loadings in each PC and their patterns) for individual principal components. Clustering was performed at two levels—first we performed coarse clustering using 25 PCs and a resolution of 3. After coarse clustering, remaining nuclei were re-normalized, scaled, and 50 principal components were used for dimensionality reduction, with a resolution of 3.

Clusters were visualized using Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP), and cluster markers were found using the “auroc” test in Seurat. Clusters with less than 3 significant markers and had low nUMI, or that were not defined by a cohesive set of genes and had low nUMI, were identified as low-quality clusters and discarded from downstream analysis (Supplementary Fig. 14). We used clustering to identify doublets, rather than defining them on a cell-by-cell basis to avoid discarding cells that may have hybrid or continuum status between two related cell types. All clusters were tested for potential doublet status by examining marker lists for defining genes, using DoubletFinder17, and by looking for co-expression of the top markers using FeatureScatter (using top ten cell type markers from a previously published spinal cord atlas19). If these analyses supported a doublet identity for the cluster, it was removed from downstream analysis. Overall, 22,435 nuclei were discarded as low-quality or doublets. We next subclustered the neurons, oligodendrocytes, astrocytes, and microglia independently. Raw data for the nuclei in each class was re-analyzed with standard log normalization and a new principal component analysis. We used the following principal components and resolution for each subclustering: neurons: 40 dimensions and resolution 1; oligodendrocytes: 16 dimensions and resolution 0.3; astrocytes: 8 dimensions and resolution 0.3; and microglia: 11 dimensions and resolution 0.8. In addition, we also used label transfer88 to analyze the neurons and oligodendrocytes, as described in the Supplemental data. Nuclei had on average 1,392 genes per nucleus in neurons and 471 genes per nucleus in non-neuronal cells (Supplementary Dataset 7).

Single nucleus ATAC sequencing

Single nucleus ATAC sequencing was carried out using Chromium Next GEM Single Cell ATAC v.1.1 kit on the Chromium platform (10X Genomics) according to manufacturer’s instructions. Libraries were sequenced to a minimum depth of 10,000 reads per nucleus using an Illumina MiSeq (PE 50–8–16–50 bp). Raw sequencing reads were demultiplexed, aligned, and a count matrix was generated using CellRanger-atac 2.0. Cell type-specific dimensionality reduction and cluster analysis for snATAC-seq was performed using ArchR (version 1.0.1). To cluster our scATAC-seq data (for both broad clustering and neuronal subclustering), we used ArchR’s addIterativeLSI function to perform iterative LSI clustering.

Cell types were determined using Seurat’s label-transfer algorithm with cell type annotations in snRNA-seq cells as the reference. Neurons were further subclustered and annotated into “families” (DE, DI, ME, MI, VE, VI, MN) using Seurat’s label-transfer algorithm.

Gene activity scores were calculated using ArchR v1.0.1 with default parameters by using a distance-weighted accessibility model that aggregates signal inside the gene body and in the local genomic region89. The resulting gene activity scores were additionally imputed using MAGIC90 to reduce sparsity noise in the scATAC-seq data. For peak calling and sequencing tracks, we used the addReproduciblePeakSet function from ArchR (v.1.0.1) with default parameters to call accessible chromatin peaks using MACS2 (v.2.2.7.1) in each cell type subcluster. Marker peaks were identified using ArchR’s getMarkerFeatures function. Sequencing tracks were created using ArchR’s plotBrowserTrack function. All tracks show data that have been normalized by ‘reads-in-TSS’ to account for differences in signal-to-background ratios across samples, unless otherwise stated. For all tracks, genes on the plus strand are shown in red and genes on the minus strand are shown in blue.

Cell type prioritization by AUGUR and GO analysis

Augur was implemented as previously described using default parameters to rank which cell types changed the most after injury66. This approach uses a random forest classifier on subsampled matrices and reports the mean cross-validation AUC across many small subsamples (code is available on GitHub, see below). (The AUC is a measure of the performance of a classifier, with 1 being a perfect classification, 0 being random and negative values indicating poor performance.) For pathway analysis, differential gene expression across conditions was generated using FindMarkers using the Wilcox test. GO Analysis was done using all differentially expressed genes with p_adj < 0.05 using medium stringency and default parameters at https://david.ncifcrf.gov/. GO biological process, cell compartment, and molecular function were analyzed and the clustering annotation tool was selected. Only clusters with an enrichment score (−log of p value) greater than 1.3 were considered. In cases in which multiple clusters had the same genes and similar terms, only the most significant is shown.

Immunohistochemistry and in situ hybridization

Animals were euthanized with avertin and perfused with PBS and then 4% paraformaldehyde. Spinal cords were dissected, fixed in 4% paraformaldehyde overnight, washed in PBS for one hour, then dehydrated in 30% sucrose an additional night before being embedded in OCT.

Immunohistochemistry was performed as previously described18. Briefly, spinal cords were cut at 50 µm, placed in blocking buffer (1% IgG-free BSA, 10% normal donkey serum, 0.1% Triton X-100 in PBS) for one hour prior to incubation in blocking buffer and primary antibody for 48 h at 4 °C. Primary antibody was washed off three times in PBS before a 2-h incubation in secondary antibody at room temperature. The secondary antibody was washed off three times in PBS before adding a coverslip.

Multiplex immunohistochemistry was performed as previously described on 10-µm-thick tissue sections91. In situ hybridization was performed according to the manufacturer’s instructions for fixed frozen tissue RNAscope (ACD Bio).

Immunohistochemistry antibodies

Primary antibodies

IBA1 (dilution 1:100, Cedarlane Labs, 234006(SY)), TMEM119 (dilution 1:100, Cedarlane Labs, 400004(SY)), CD11c (dilution 1:100, GeneTex, GTX74935), Myelin-MBP (dilution 1:100, BioLegend, 808402), NF-L (dilution 1:200, Cell Signaling, 2835S), NF-M (dilution 1:200, Cell Signaling, 2838 S), NF-H (dilution 1:200, Cell Signaling, 2836 S), NeuN (dilution 1:500, Millipore Sigma, ABN90P), CD68 (dilution 1:100, Abcam, ab125212), CNPase (dilution 1:200, Millipore Sigma, MAB326), GFAP (dilution 1:500, Agilent/Dako, Z033429-2), Cleaved Caspase 3 (dilution 1:100, Cell Signalling Tech, 9661L), and Phospho-IGF1R (dilution 1:100, Invitrogen, PA5-104773).

Secondary antibodies

DAPI (dilution 1:1000, Thermo Fisher Scientific, 62248), Goat anti-Hamster IgG Alexa Fluor 647 (dilution 1:200, Jackson ImmunoResearch, 127-605-160), Goat anti-Rabbit IgG Alexa Fluor 430 (dilution 1:200, Thermo Fisher Scientific, A11064), Goat anti-Rabbit IgG Alexa Fluor 430 (dilution 1:200, Thermo Fisher Scientific, A11064), Goat anti-Rabbit IgG Alexa Fluor 594 (dilution 1:200, Thermo Fisher Scientific, A11037), Goat anti-Mouse IgG1 Alexa Fluor 488 (dilution 1:200, Thermo Fisher Scientific, A21121), Goat anti-Mouse IgG1 Alexa Fluor 647 (dilution 1:200, Thermo Fisher Scientific, A21240), Goat anti-Mouse IgG2a Alexa Fluor 546 (dilution 1:200, Thermo Fisher Scientific, A21133), Goat anti-Mouse IgG2b Alexa Fluor 647 (dilution 1:200, Thermo Fisher Scientific, A21242), Goat anti-Guinea Pig IgG Alexa Fluor 546 (dilution 1:200, Thermo Fisher Scientific, A11074), Goat anti-Guinea Pig IgG Alexa Fluor 594 (dilution 1:200, Thermo Fisher Scientific, A11076), Donkey anti-Chicken IgG IRDye 680LT (dilution 1:200, Li-Cor Biosciences, 926-68028), Donkey anti-Rabbit IgG IRDye 800CW (dilution 1:200, Li-Cor Biosciences, 926-32213), Donkey anti-Chicken IgG Alexa Fluor 488 (dilution 1:500, Jackson ImmuoResearch, 703-545-155), and Donkey anti-Guinea Pig IgG Alexa Fluor 647 (dilution 1:500, Jackson ImmuoResearch, 706-605-148).

RNAscope in situ hybridization probes

ACDbio RNAscope probes: Spp1 (435191), Vsx2 (438341), Gap43 (318621), Chat (408731), Itgax (311501), Mdga1 (546411), Sprr1a (426871-C2), Vgf (517421-C2), Igf1 (443901-C2), Sprr1a (426871-C2), C1qa (441221-C2), Megf11 (504281-C2), GFP (409011-C2), Apoe (313271-C3), Tnfrsf12a (429311-C3), Atf3 (426891-C3), Gpr83 (317431-C3), Pdgfra (480661-C3), Shox2 (554291-C3), Gpnmb (489511-C3).

Imaging

Images of immunohistochemistry and in situ hybridization samples were imaged using a Zeiss 800 LSM confocal microscope. For quantification, a tile scan image spanning the section was generated for ≥ 3 sections from ≥ 3 mice. Brightness and contrast were adjusted in Photoshop (Adobe), standardized across images.

Quantification of cell counts from images of immunohistochemistry

Quantification of NeuN, Olig2, and DAPI immunohistochemistry for Fig. 1 was done through a custom MATLAB-based image analysis program (code is available on Github, see below). The software automatically identifies and counts cells based on a criterion that constrains size at a user-selectable intensity threshold. A manual selection tool is also available to identify additional cells that are more difficult to detect. A second channel of the same image shows can be used to count cells that are labeled with both stains using the results from the first channel and a second set of user-selectable thresholds.

Quantification of all other immunohistochemistry was done using FIJI (ImageJ) and photoshop counting tools.

Pixel quantification

Pixels were quantified using a custom python script (code is available on Github, see below), after a standardized thresholding of images using FIJI (ImageJ).

Fluorescence intensity

Mean fluorescence intensity was quantified using FIJI (ImageJ) after manually drawing borders based on DAPI and NeuN expression.

Histological quantification

Two-tailed Mann–Whitney _t_-tests (unpaired) were used for quantification of immunohistochemistry and in situ hybridization. Differences among groups were considered significant if p < 0.05. Data are represented as mean ± SEM unless otherwise indicated. Statistical analyses were performed using GraphPad Prism software (version 9.0.0). In all plots, dots represent individual mice. Imaging and most measurements were not done blinded. Only the quantification of pixels for gray matter collaterals was blinded.

Statistics and reproducibility

Sample sizes are estimated based on previous studies, with at least four animals per condition for anatomical and histological work, and three animals per condition for snRNAseq. For anatomical and histological work, four animals per condition enabled statistical analysis despite a small amount of variability expected after a spinal cord contusion injury. For snRNAseq, three animals were chosen as the minimum number to allow statistical analysis of each time point, while balancing sequencing costs.

AAV injections were replicated in two completely independent experimental cohorts. Mice were randomly assigned to injury groups and AAV-injection groups in each cohort, with AAV injections occurring 2 weeks prior to the injury surgery. Both cohorts were included in the final analysis.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Supplementary information

Acknowledgements

The authors gratefully acknowledge technical support from Li Li, Aaron Bickert, and Stefan Stoica, and advice and comments from Claire Le Pichon and Michael O’Donovan, as well as the veterinary staff of the Porter Research Building, including Heather Narver. This research was supported in part by the Intramural Research Program of the NIH, NINDS, and CIT. A.S., I.H., K.J.E.M., and A.J.L. were supported by NIH Intramural Funds through 1ZIA NS003153; and J.K. and R.P. were supported by Z99 CT999999. G.C., C.K., and J.W.S. were supported by the European Research Council [ERC-2015-CoG HOW2WALKAGAIN 682999] and the Swiss National Science Foundation (subside 310030_192558).

Author contributions

This project was conceptualized by A.J.L. and K.J.E.M. with input from C.K., J.W.S., and G.C. Experiments were performed by K.J.E.M., D.M., and C.K. Data was analyzed by K.J.E.M., D.E.R., J.W.S., J.K., R.P., I.H., Y.D., B.P.L., and A.S. The manuscript was written by K.J.E.M. and A.J.L with input from all authors. Supervision and funding acquisition were by A.J.L. and G.C.

Peer review

Peer review information

Nature Communications thanks the anonymous reviewer(s) for their contribution to the peer review of this work.

Funding

Open Access funding provided by the National Institutes of Health (NIH).

Competing interests

The authors declare the following competing financial interests: G.C. is a consultant with ONWARD Medical. G.C. is a shareholder of ONWARD Medical, a company developing an EES-based therapy to restore movement after spinal cord injury. The remaining authors declare no competing interests.

Footnotes

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

The online version contains supplementary material available at 10.1038/s41467-022-33184-1.

References

1. Kigerl KA, et al. Identification of two distinct macrophage subsets with divergent effects causing either neurotoxicity or regeneration in the injured mouse spinal cord. J. Neurosci. 2009;29:13435–13444. doi: 10.1523/JNEUROSCI.3257-09.2009. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

2. Sofroniew MV, Vinters HV. Astrocytes: biology and pathology. Acta Neuropathol. 2010;119:7–35. doi: 10.1007/s00401-009-0619-8. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

3. McDonald JW, Sadowsky C. Spinal-cord injury. Lancet. 2002;359:417–425. doi: 10.1016/S0140-6736(02)07603-1. [PubMed] [CrossRef] [Google Scholar]

4. Courtine G, et al. Recovery of supraspinal control of stepping via indirect propriospinal relay connections after spinal cord injury. Nat. Med. 2008;14:69–74. doi: 10.1038/nm1682. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

5. Ballermann M, Fouad K. Spontaneous locomotor recovery in spinal cord injured rats is accompanied by anatomical plasticity of reticulospinal fibers. Eur. J. Neurosci. 2006;23:1988–1996. doi: 10.1111/j.1460-9568.2006.04726.x. [PubMed] [CrossRef] [Google Scholar]

6. Bareyre FM, et al. The injured spinal cord spontaneously forms a new intraspinal circuit in adult rats. Nat. Neurosci. 2004;7:269–277. doi: 10.1038/nn1195. [PubMed] [CrossRef] [Google Scholar]

7. Liu Y, et al. A Sensitized IGF1 treatment restores corticospinal axon-dependent functions. Neuron. 2017;95:817–833 e814. doi: 10.1016/j.neuron.2017.07.037. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

8. Jin D, et al. Restoration of skilled locomotion by sprouting corticospinal axons induced by co-deletion of PTEN and SOCS3. Nat. Commun. 2015;6:8074. doi: 10.1038/ncomms9074. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

9. Angeli CA, et al. Recovery of over-ground walking after chronic motor complete spinal cord injury. N. Engl. J. Med. 2018;379:1244–1250. doi: 10.1056/NEJMoa1803588. [PubMed] [CrossRef] [Google Scholar]

10. Gill ML, et al. Neuromodulation of lumbosacral spinal networks enables independent stepping after complete paraplegia. Nat. Med. 2018;24:1677–1682. doi: 10.1038/s41591-018-0175-7. [PubMed] [CrossRef] [Google Scholar]

11. Wagner FB, et al. Targeted neurotechnology restores walking in humans with spinal cord injury. Nature. 2018;563:65–71. doi: 10.1038/s41586-018-0649-2. [PubMed] [CrossRef] [Google Scholar]

12. Courtine G, Sofroniew MV. Spinal cord repair: advances in biology and technology. Nat. Med. 2019;25:898–908. doi: 10.1038/s41591-019-0475-6. [PubMed] [CrossRef] [Google Scholar]

13. Sofroniew MV. Dissecting spinal cord regeneration. Nature. 2018;557:343–350. doi: 10.1038/s41586-018-0068-4. [PubMed] [CrossRef] [Google Scholar]

14. Nguyen, M. Q., Le Pichon, C. E. & Ryba, N. Stereotyped transcriptomic transformation of somatosensory neurons in response to injury. Elife8, 10.7554/eLife.49679 (2019). [PMC free article] [PubMed]

15. Renthal W, et al. Transcriptional reprogramming of distinct peripheral sensory neuron subtypes after axonal injury. Neuron. 2020;108:128–144 e129. doi: 10.1016/j.neuron.2020.07.026. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

16. Kaplan A, Ong Tone S, Fournier AE. Extrinsic and intrinsic regulation of axon regeneration at a crossroads. Front. Mol. Neurosci. 2015;8:27. doi: 10.3389/fnmol.2015.00027. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

17. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 2019;8:329–337 e324. doi: 10.1016/j.cels.2019.03.003. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

18. Sathyamurthy A, et al. Massively parallel single nucleus transcriptional profiling defines spinal cord neurons and their activity during behavior. Cell Rep. 2018;22:2216–2225. doi: 10.1016/j.celrep.2018.02.003. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

19. Russ DE, et al. A harmonized atlas of mouse spinal cord cell types and their spatial organization. Nat. Commun. 2021;12:5722. doi: 10.1038/s41467-021-25125-1. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

20. Dobrott CI, Sathyamurthy A, Levine AJ. Decoding cell type diversity within the spinal cord. Curr. Opin. Physiol. 2019;8:1–6. doi: 10.1016/j.cophys.2018.11.006. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

21. Squair JW, Gautier M, Sofroniew MV, Courtine G, Anderson MA. Engineering spinal cord repair. Curr. Opin. Biotechnol. 2021;72:48–53. doi: 10.1016/j.copbio.2021.10.006. [PubMed] [CrossRef] [Google Scholar]

22. Matson KJE, et al. Isolation of adult spinal cord nuclei for massively parallel single-nucleus RNA sequencing. J. Vis. Exp. 2018 doi: 10.3791/58413. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

23. Buttner M, Ostner J, Muller CL, Theis FJ, Schubert B. scCODA is a Bayesian model for compositional single-cell data analysis. Nat. Commun. 2021;12:6876. doi: 10.1038/s41467-021-27150-6. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

24. Liddelow SA, et al. Neurotoxic reactive astrocytes are induced by activated microglia. Nature. 2017;541:481–487. doi: 10.1038/nature21029. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

25. Marques S, et al. Oligodendrocyte heterogeneity in the mouse juvenile and adult central nervous system. Science. 2016;352:1326–1329. doi: 10.1126/science.aaf6463. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

26. Zeisel A, et al. Molecular architecture of the mouse nervous system. Cell. 2018;174:999–1014 e1022. doi: 10.1016/j.cell.2018.06.021. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

27. Assinck P, et al. Myelinogenic plasticity of oligodendrocyte precursor cells following spinal cord contusion injury. J. Neurosci. 2017;37:8635–8654. doi: 10.1523/JNEUROSCI.2409-16.2017. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

28. Hesp ZC, Goldstein EZ, Miranda CJ, Kaspar BK, McTigue DM. Chronic oligodendrogenesis and remyelination after spinal cord injury in mice and rats. J. Neurosci. 2015;35:1274–1290. doi: 10.1523/JNEUROSCI.2568-14.2015. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

29. Marisca R, et al. Functionally distinct subgroups of oligodendrocyte precursor cells integrate neural activity and execute myelin formation. Nat. Neurosci. 2020;23:363–374. doi: 10.1038/s41593-019-0581-2. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

30. Floriddia EM, et al. Distinct oligodendrocyte populations have spatial preference and different responses to spinal cord injury. Nat. Commun. 2020;11:5860. doi: 10.1038/s41467-020-19453-x. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

31. Yamamoto T, et al. SPP1 expression in spinal motor neurons of the macaque monkey. Neurosci. Res. 2011;69:81–86. doi: 10.1016/j.neures.2010.09.010. [PubMed] [CrossRef] [Google Scholar]

32. Alkaslasi MR, et al. Single nucleus RNA-sequencing defines unexpected diversity of cholinergic neuron types in the adult mouse spinal cord. Nat. Commun. 2021;12:2471. doi: 10.1038/s41467-021-22691-2. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

33. Wlodarczyk A, et al. A novel microglial subset plays a key role in myelinogenesis in developing brain. EMBO J. 2017;36:3292–3308. doi: 10.15252/embj.201696056. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

34. Jessen KR, Mirsky R. The success and failure of the Schwann cell response to nerve injury. Front. Cell. Neurosci. 2019;13:33. doi: 10.3389/fncel.2019.00033. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

35. Liu XZ, et al. Neuronal and glial apoptosis after traumatic spinal cord injury. J. Neurosci. 1997;17:5395–5406. doi: 10.1523/JNEUROSCI.17-14-05395.1997. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

36. D’Amico JM, Condliffe EG, Martins KJ, Bennett DJ, Gorassini MA. Recovery of neuronal and network excitability after spinal cord injury and implications for spasticity. Front Integr. Neurosci. 2014;8:36. [PMC free article] [PubMed] [Google Scholar]

37. Wang Z, Maunze B, Wang Y, Tsoulfas P, Blackmore MG. Global connectivity and function of descending spinal input revealed by 3D microscopy and retrograde transduction. J. Neurosci. 2018;38:10566–10581. doi: 10.1523/JNEUROSCI.1196-18.2018. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

38. Holland, S. D., Ramer, L. M., McMahon, S. B., Denk, F. & Ramer, M. S. An ATF3-CreERT2 knock-in mouse for axotomy-induced genetic editing: proof of principle. eNeuro6, 10.1523/ENEURO.0025-19.2019 (2019). [PMC free article] [PubMed]

39. Starkey ML, et al. Expression of the regeneration-associated protein SPRR1A in primary sensory neurons and spinal cord of the adult mouse following peripheral and central injury. J. Comp. Neurol. 2009;513:51–68. doi: 10.1002/cne.21944. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

40. Bonilla IE, Tanabe K, Strittmatter SM. Small proline-rich repeat protein 1A is expressed by axotomized neurons and promotes axonal outgrowth. J. Neurosci. 2002;22:1303–1315. doi: 10.1523/JNEUROSCI.22-04-01303.2002. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

41. Kim JE, Liu BP, Park JH, Strittmatter SM. Nogo-66 receptor prevents raphespinal and rubrospinal axon regeneration and limits functional recovery from spinal cord injury. Neuron. 2004;44:439–451. doi: 10.1016/j.neuron.2004.10.015. [PubMed] [CrossRef] [Google Scholar]

42. Forstner P, et al. Neuroinflammation after traumatic brain injury is enhanced in activating transcription factor 3 mutant mice. J. Neurotrauma. 2018;35:2317–2329. doi: 10.1089/neu.2017.5593. [PubMed] [CrossRef] [Google Scholar]

43. Dhara SP, et al. Cellular reprogramming for successful CNS axon regeneration is driven by a temporally changing cast of transcription factors. Sci. Rep. 2019;9:14198. doi: 10.1038/s41598-019-50485-6. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

44. Haring M, et al. Neuronal atlas of the dorsal horn defines its architecture and links sensory input to transcriptional cell types. Nat. Neurosci. 2018;21:869–880. doi: 10.1038/s41593-018-0141-1. [PubMed] [CrossRef] [Google Scholar]

45. Barik, A. et al. A spinoparabrachial circuit defined by Tacr1 expression drives pain. Elife10, 10.7554/eLife.61135 (2021). [PMC free article] [PubMed]

46. Roome RB, et al. Phox2a defines a developmental origin of the anterolateral system in mice and humans. Cell Rep. 2020;33:108425. doi: 10.1016/j.celrep.2020.108425. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

47. Chiang MC, et al. Divergent neural pathways emanating from the lateral parabrachial nucleus mediate distinct components of the pain response. Neuron. 2020;106:927–939 e925. doi: 10.1016/j.neuron.2020.03.014. [PubMed] [CrossRef] [Google Scholar]

48. Arber S. Motor circuits in action: specification, connectivity, and function. Neuron. 2012;74:975–989. doi: 10.1016/j.neuron.2012.05.011. [PubMed] [CrossRef] [Google Scholar]

49. Gosgnach S, et al. Delineating the diversity of spinal interneurons in locomotor circuits. J. Neurosci. 2017;37:10835–10841. doi: 10.1523/JNEUROSCI.1829-17.2017. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

50. Ha, N. T. & Dougherty, K. J. Spinal Shox2 interneuron interconnectivity related to function and development. Elife7, 10.7554/eLife.42519 (2018). [PMC free article] [PubMed]

51. Dougherty KJ, et al. Locomotor rhythm generation linked to the output of spinal shox2 excitatory interneurons. Neuron. 2013;80:920–933. doi: 10.1016/j.neuron.2013.08.015. [PubMed] [CrossRef] [Google Scholar]

52. Garcia-Ramirez DL, Ha NTB, Bibu S, Stachowski NJ, Dougherty KJ. Spinal cord injury alters spinal Shox2 interneurons by enhancing excitatory synaptic input and serotonergic modulation while maintaining intrinsic properties in mouse. J. Neurosci. 2021 doi: 10.1523/JNEUROSCI.1576-20.2021. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

53. Baek M, Menon V, Jessell TM, Hantman AW, Dasen JS. Molecular logic of spinocerebellar tract neuron diversity and connectivity. Cell Rep. 2019;27:2620–2635 e2624. doi: 10.1016/j.celrep.2019.04.113. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

54. O’Shea TM, Burda JE, Sofroniew MV. Cell biology of spinal cord injury and repair. J. Clin. Investig. 2017;127:3259–3270. doi: 10.1172/JCI90608. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

55. Tran AP, Warren PM, Silver J. The biology of regeneration failure and success after spinal cord injury. Physiol. Rev. 2018;98:881–917. doi: 10.1152/physrev.00017.2017. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

56. Milich LM, et al. Single-cell analysis of the cellular heterogeneity and interactions in the injured mouse spinal cord. J. Exp. Med. 2021;218:e20210040. doi: 10.1084/jem.20210040. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

57. Wahane, S. et al. Diversified transcriptional responses of myeloid and glial cells in spinal cord injury shaped by HDAC3 activity. Sci. Adv.7, 10.1126/sciadv.abd8811 (2021). [PMC free article] [PubMed]

58. Keren-Shaul H, et al. A unique microglia type associated with restricting development of Alzheimer’s disease. Cell. 2017;169:1276–1290 e1217. doi: 10.1016/j.cell.2017.05.018. [PubMed] [CrossRef] [Google Scholar]

59. Hakim R, et al. Spinal cord injury induces permanent reprogramming of microglia into a disease-associated state which contributes to functional recovery. J. Neurosci. 2021;41:8441–8459. doi: 10.1523/JNEUROSCI.0860-21.2021. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

60. Li Q, et al. Developmental heterogeneity of microglia and brain myeloid cells revealed by deep single-cell RNA sequencing. Neuron. 2019;101:207–223 e210. doi: 10.1016/j.neuron.2018.12.006. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

61. Li Y, et al. Microglia-organized scar-free spinal cord repair in neonatal mice. Nature. 2020;587:613–618. doi: 10.1038/s41586-020-2795-6. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

62. Hammond TR, et al. Single-Cell RNA sequencing of microglia throughout the mouse lifespan and in the injured brain reveals complex cell-state changes. Immunity. 2019;50:253–271 e256. doi: 10.1016/j.immuni.2018.11.004. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

63. Benmamar-Badel A, Owens T, Wlodarczyk A. Protective microglial subset in development, aging, and disease: lessons from transcriptomic studies. Front. Immunol. 2020;11:430. doi: 10.3389/fimmu.2020.00430. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

64. Kamphuis W, Kooijman L, Schetters S, Orre M, Hol EM. Transcriptional profiling of CD11c-positive microglia accumulating around amyloid plaques in a mouse model for Alzheimer’s disease. Biochim. Biophys. Acta. 2016;1862:1847–1860. doi: 10.1016/j.bbadis.2016.07.007. [PubMed] [CrossRef] [Google Scholar]

65. Tansley S, et al. Single-cell RNA sequencing reveals time- and sex-specific responses of mouse spinal cord microglia to peripheral nerve injury and links ApoE to chronic pain. Nat. Commun. 2022;13:843. doi: 10.1038/s41467-022-28473-8. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

66. Skinnider MA, et al. Cell type prioritization in single-cell data. Nat. Biotechnol. 2021;39:30–34. doi: 10.1038/s41587-020-0605-1. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

67. Ma TC, Willis DE. What makes a RAG regeneration associated? Front Mol. Neurosci. 2015;8:43. doi: 10.3389/fnmol.2015.00043. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

68. Bomze HM, Bulsara KR, Iskandar BJ, Caroni P, Skene JH. Spinal axon regeneration evoked by replacing two growth cone proteins in adult neurons. Nat. Neurosci. 2001;4:38–43. doi: 10.1038/82881. [PubMed] [CrossRef] [Google Scholar]

69. Fagoe ND, Attwell CL, Kouwenhoven D, Verhaagen J, Mason MR. Overexpression of ATF3 or the combination of ATF3, c-Jun, STAT3 and Smad1 promotes regeneration of the central axon branch of sensory neurons but without synergistic effects. Hum. Mol. Genet. 2015;24:6788–6800. doi: 10.1093/hmg/ddv383. [PubMed] [CrossRef] [Google Scholar]

70. Jing X, Wang T, Huang S, Glorioso JC, Albers KM. The transcription factor Sox11 promotes nerve regeneration through activation of the regeneration-associated gene Sprr1a. Exp. Neurol. 2012;233:221–232. doi: 10.1016/j.expneurol.2011.10.005. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

71. Collyer E, et al. Sprouting of axonal collaterals after spinal cord injury is prevented by delayed axonal degeneration. Exp. Neurol. 2014;261:451–461. doi: 10.1016/j.expneurol.2014.07.014. [PubMed] [CrossRef] [Google Scholar]

72. Cafferty WB, McGee AW, Strittmatter SM. Axonal growth therapeutics: regeneration or sprouting or plasticity? Trends Neurosci. 2008;31:215–220. doi: 10.1016/j.tins.2008.02.004. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

73. Jain N, Florence SL, Kaas JH. Limits on plasticity in somatosensory cortex of adult rats: hindlimb cortex is not reactivated after dorsal column section. J. Neurophysiol. 1995;73:1537–1546. doi: 10.1152/jn.1995.73.4.1537. [PubMed] [CrossRef] [Google Scholar]

74. Barriere G, Leblond H, Provencher J, Rossignol S. Prominent role of the spinal central pattern generator in the recovery of locomotion after partial spinal cord injuries. J. Neurosci. 2008;28:3976–3987. doi: 10.1523/JNEUROSCI.5692-07.2008. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

75. Chalif JI, Martinez-Silva ML, Pagiazitis JG, Murray AJ, Mentis GZ. Control of mammalian locomotion by ventral spinocerebellar tract neurons. Cell. 2022;185:328–344 e326. doi: 10.1016/j.cell.2021.12.014. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

76. Lei Y, Perez MA. Cerebellar contribution to sensorimotor adaptation deficits in humans with spinal cord injury. Sci. Rep. 2021;11:2507. doi: 10.1038/s41598-020-77543-8. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

77. Lacar B, et al. Nuclear RNA-seq of single neurons reveals molecular signatures of activation. Nat. Commun. 2016;7:11022. doi: 10.1038/ncomms11022. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

78. Wu YE, Pan L, Zuo Y, Li X, Hong W. Detecting activated cell populations using single-cell RNA-Seq. Neuron. 2017;96:313–329 e316. doi: 10.1016/j.neuron.2017.09.026. [PubMed] [CrossRef] [Google Scholar]

79. Bakken TE, et al. Single-nucleus and single-cell transcriptomes compared in matched cortical cell types. PLoS ONE. 2018;13:e0209648. doi: 10.1371/journal.pone.0209648. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

80. Thrupp N, et al. Single-nucleus RNA-Seq is not suitable for detection of microglial activation genes in humans. Cell Rep. 2020;32:108189. doi: 10.1016/j.celrep.2020.108189. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

81. Marsh SE, et al. Dissection of artifactual and confounding glial signatures by single-cell sequencing of mouse and human brain. Nat. Neurosci. 2022;25:306–316. doi: 10.1038/s41593-022-01022-8. [PubMed] [CrossRef] [Google Scholar]

82. Asboth L, et al. Cortico-reticulo-spinal circuit reorganization enables functional recovery after severe spinal cord contusion. Nat. Neurosci. 2018;21:576–588. doi: 10.1038/s41593-018-0093-5. [PubMed] [CrossRef] [Google Scholar]

83. Madisen L, et al. Transgenic mice for intersectional targeting of neural sensors and effectors with high specificity and performance. Neuron. 2015;85:942–958. doi: 10.1016/j.neuron.2015.02.022. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

84. Sathyamurthy A, et al. Cerebellospinal neurons regulate motor performance and motor learning. Cell Rep. 2020;31:107595. doi: 10.1016/j.celrep.2020.107595. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

85. Chaterji S, Barik A, Sathyamurthy A. Intraspinal injection of adeno-associated viruses into the adult mouse spinal cord. STAR Protoc. 2021;2:100786. doi: 10.1016/j.xpro.2021.100786. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

86. Harrison M, et al. Vertebral landmarks for the identification of spinal cord segments in the mouse. Neuroimage. 2013;68:22–29. doi: 10.1016/j.neuroimage.2012.11.048. [PubMed] [CrossRef] [Google Scholar]

87. Sengul G, Fu Y, Yu Y, Paxinos G. Spinal cord projections to the cerebellum in the mouse. Brain Struct. Funct. 2015;220:2997–3009. doi: 10.1007/s00429-014-0840-7. [PubMed] [CrossRef] [Google Scholar]

88. Stuart T, et al. Comprehensive integration of single-cell data. Cell. 2019;177:1888–1902 e1821. doi: 10.1016/j.cell.2019.05.031. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

89. Granja JM, et al. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat. Genet. 2021;53:403–411. doi: 10.1038/s41588-021-00790-6. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

90. Roopra A. MAGIC: a tool for predicting transcription factors and cofactors driving gene sets using ENCODE data. PLoS Comput. Biol. 2020;16:e1007800. doi: 10.1371/journal.pcbi.1007800. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

91. Maric D, et al. Whole-brain tissue mapping toolkit using large-scale highly multiplexed immunofluorescence imaging and deep neural networks. Nat. Commun. 2021;12:1550. doi: 10.1038/s41467-021-21735-x. [PMC free article] [PubMed] [CrossRef] [Google Scholar]

Articles from Nature Communications are provided here courtesy of Nature Publishing Group