The neurons that restore walking after paralysis

Abstract

A spinal cord injury interrupts pathways from the brain and brainstem that project to the lumbar spinal cord, leading to paralysis. Here we show that spatiotemporal epidural electrical stimulation (EES) of the lumbar spinal cord1-3 applied during neurorehabilitation4,5 (EESREHAB) restored walking in nine individuals with chronic spinal cord injury. This recovery involved a reduction in neuronal activity in the lumbar spinal cord of humans during walking. We hypothesized that this unexpected reduction reflects activity-dependent selection of specific neuronal subpopulations that become essential for a patient to walk after spinal cord injury. To identify these putative neurons, we modelled the technological and therapeutic features underlying EESREHAB in mice. We applied single-nucleus RNA sequencing6-9 and spatial transcriptomics10,11 to the spinal cords of these mice to chart a spatially resolved molecular atlas of recovery from paralysis. We then employed cell type12,13 and spatial prioritization to identify the neurons involved in the recovery of walking. A single population of excitatory interneurons nested within intermediate laminae emerged. Although these neurons are not required for walking before spinal cord injury, we demonstrate that they are essential for the recovery of walking with EES following spinal cord injury. Augmenting the activity of these neurons phenocopied the recovery of walking enabled by EESREHAB, whereas ablating them prevented the recovery of walking that occurs spontaneously after moderate spinal cord injury. We thus identified a recovery-organizing neuronal subpopulation that is necessary and sufficient to regain walking after paralysis. Moreover, our methodology establishes a framework for using molecular cartography to identify the neurons that produce complex behaviours.

Conflict of interest statement

The authors declare competing interests: G.C., J.B., F.W., L.A., R.D., S.K., S.P.L. and J.W.S. hold various patents in relation to the present work. G.C. is a consultant for ONWARD medical. G.C., J.B. and S.P.L. are minority shareholders of ONWARD, a company related to the presented work. The other authors declare no competing interests.

© 2022. The Author(s).

Figures

Fig. 1. EESREHAB remodels the spinal cord of humans and mice.

a, Body weight support system enabling overground walking and wireless implantable pulse generator operating in closed loop, connected to a paddle lead targeting the dorsal roots that innervate lumbosacral segments. b, Chronophotography showing the transitioning from sitting to walking in a representative participant. c, 18FDG-PET projected onto a personalized model of the spinal cord elaborated from high-resolution MRI (participant ID DM002), showing the metabolic activity of the spinal cord—expressed as standardized uptake value (SUVbw)—in response to walking before and after EESREHAB. d, Bar plots reporting the relative change in normalized FDG-PET metabolic activity during walking before and after EESREHAB, the lower limb motor scores, and the distance covered during the 6-min walk test (n = 9; metabolic activity mixed-effects model: t = −3.2, P = 0.002; lower limb motor scores, paired samples two-tailed t-test: t = 3.7, P = 0.0063; distance covered, paired samples two-tailed t-test: t = 3.5; P = 0.0076). e, Left, body weight support system enabling overground walking in mice, with implantable electrodes to deliver EES. Right, spinal cord visualization of projections from neurons in the motor cortex and glutamatergic (vGluT2ON) neurons in the reticular formation, traced with AAV5-CAG-COMET-GFP and AAV5-CAG-DIO-COMET-tdTomato, respectively. Scale bars, 1 mm. f, Chronophotography of representative mice with SCI only (SCI, EESOFF) or SCI with EESREHAB (EESREHAB, EESOFF). g, Lumbar spinal cord expression of cFos following walking with EESON following SCI or SCI with EESREHAB. Scale bars, 500 μm. h, Walking performance of uninjured mice (n = 3), mice with SCI (n = 10), and mice with SCI and EESREHAB tested with EESOFF (n = 10) or EESON (n = 10) (one-way ANOVA; Tukey’s honest significant difference for SCI versus EESREHAB→EESOFF: P = 3.3 × 10–11). i, The number of neurons expressing cFos(cFosON) (mice with SCI with EESON, n = 4; mice with EESREHAB and EESON, n = 4; independent samples two-tailed t-test: t = –5.7; P = 0.001). h,i, Bars show mean ± s.e.m. with individual points overlaid. *P < 0.05, **P < 0.01, ***P < 0.001. Source data

Fig. 2. Molecular cartography of EESREHAB.

a, Overview of the eight experimental conditions capturing the key therapeutic features of EESREHAB. A detailed description is provided in  Methods, ‘Experimental conditions’. b, Uniform manifold approximation and projection (UMAP) visualization of 20,990 nuclei revealing 36 neuron subpopulations. Five dorsal and ventral populations are highlighted on the basis of their marker genes. In each corner, an UMAP visualization coloured by the expression of classical marker genes reveals the cardinal organization of neuronal subpopulations across dorsal–ventral and excitatory–inhibitory axes. MN, motor neuron; VI, ventral-inhibitory; VE, ventral-excitatory; CSF-N, cerebrospinal-fluid contacting neurons; Ia-IN, Ia inhibitory interneurons; Rora-I, inhibitory neurons expressing Rora; Rorb-I, inhibitory neurons expressing Rorb. c, Left, visualization of 22,127 barcodes registered to a common coordinate framework highlighting the expression of classical excitatory–inhibitory and ventral–dorsal marker genes. Second from left, the localization of all 36 neuron subpopulations, with each spatial barcode coloured according to cellular identity. Five classical dorsal and ventral populations are highlighted, with the spatial expression of their marker genes shown below the image. Right, RNAscope analysis, confirming the spatial location of these five neuronal subpopulations.

Fig. 3. Molecular compass to identify recovery-organizing neurons.

a, UMAP visualization of 20,990 neurons, coloured by Augur cell type prioritization, identifying perturbation-responsive subpopulations in a representative experimental comparison. b, Identification of perturbation-responsive subpopulations across experimental comparisons with Augur. Top, clustering tree of neuronal subpopulations. Middle, heat map showing scaled areas under the curve (AUCs) for individual comparisons: (1) SCI versus SCI→EES::walking; (2) SCI versus EESREHAB; (3) EESREHAB versus EESREHAB→cortex; (4) EESREHAB versus EESREHAB→EES; (5) EESREHAB→EES::walking versus SCI→EES::walking; (6) EESREHAB versus EESREHAB→EES::walking. Bottom, distribution of AUCs across all comparisons. c, Spatial visualization of 22,127 barcodes in the common coordinate space coloured by Magellan spatial prioritization for the same representative experimental comparison as in a. d, Synaptic inputs and outputs from SCVsx2::Hoxa10 neurons. Inputs: SCVsx2::Hoxa10 neurons and their projections. Projections from vGi neurons (AAV5-CAG-COMET-GFP) and large-diameter afferent neurons (PVcre::Advillin:tdTomato) onto SCVsx2::Hoxa10 neurons. Outputs: synaptic appositions from SCVsx2::Hoxa10 neurons (AAV-Dj-hSyn-flex-mGFP-2A-synaptophysin-mRuby) onto glutamatergic, GABAergic and motor neurons. C, caudal; R, rostral; PV, parvalbumin; syn, synaptophysin. e, Single-unit recordings of optotagged SCVsx2::Hoxa10 neurons. Responses of an optotagged single unit to optogenetic stimulations. Right, the number of units responding to each type of stimulation. f, Number of vGluT1 synapses from large-diameter afferents apposing SCVsx2::Hoxa10 neurons (n = 6 out of 10 mice per group; independent samples two-tailed t-test: t = 4.9, P = 0.002). g, Number of reticulospinal neurons projecting to SCVsx2::Hoxa10 neurons (n = 4 mice per group; independent samples two-tailed t-test: t = 4.8, P = 0.0029). h, Percentage of motor neurons, glutamatergic and GABAergic neurons receiving at least one synapse from SCVsx2::Hoxa10 neurons (n = 6; one-way analysis of variance (ANOVA); choline acetyltransferase (ChAT): F = 0.45, P = 0.39; glutamatergic: F = 0.52, P = 0.76; GABAergic: statistics indicate Tukey’s honest significant difference tests: uninjured versus SCI: P = 5.2 × 10−6; SCI versus EESREHAB: P = 1.8 × 10−5). i, Fluorescent intensity of cFos in Vsx2cre neurons following walking with EESON (n = 5 out of 3 mice per group; independent samples two-tailed t-test: t = 5.7, P = 0.0013). f–i, Bars show mean ± s.e.m. with individual points overlaid. AU, arbitrary units. Source data

Fig. 4. SCVsx2::Hoxa10 neurons are required for the recovery of walking after paralysis.

a, Implantable optoelectronic device to deliver EES and photostimulation. After EESREHAB, chronophotography of walking during EESON in a representative mouse. Red-shifted light is delivered for a few seconds to silence SCVsx2::Hoxa10 neurons (AAV5-Syn-flex-ChrimsonR-tdTomato; n = 4; Tukey’s honest significant difference tests following repeated measures one-way analysis of variance (ANOVA): P = 0.0023). b, Chronophotography of walking after EESREHAB in a representative mouse with EESON before and after chemogenetic silencing of SCVsx2::Hoxa10 neurons (AAV5-hSyn-DIO-hm4D-(Gi)-mCherry; n = 4; paired samples two-tailed t-test; t = −21.3, P = 0.0002). c, Chronophotography of walking in representative mice with no EESREHAB, with EESOFF before and after chemogenetic activation of SCVsx2::Hoxa10 neurons (AAV5-hSyn-DIO-hm3D-(Gq)-mCherry; n = 4; paired samples two-tailed t-test; t = 5.3, P = 0.0013). d, After EESREHAB, chronophotography of walking in representative mice with EESON. SCVsx2::Hoxa10 neurons were silenced during the entire period of EESREHAB (n = 5; independent samples two-tailed t-test: t = −3.5, P = 0.008). e, Chronophotography of representative mice walking without EES after recovery from a lateral hemisection SCI. Kinematic limb reconstruction is overlaid. Right, same condition, but SCVsx2::Hoxa10 neurons located in lumbar segments were ablated before the SCI (AAV5-CAG-flex-DTR; n = 5 per group; independent samples two-tailed t-test: t = 5.9, P = 0.0004). a–e, Bars show mean ± s.e.m. with individual points overlaid. Source data

Extended Data Fig. 1. EESREHAB restores walking in nine humans with chronic SCI.

a, Chronophotography of nine human participants with chronic SCI who are walking independently with EESON and the help of crutches or a front-wheel walker at the end of the neurorehabilitation period. b, Bar plots reporting the volitional modulation of the step for each participant during EESON (paired samples two-tailed t-tests: P9: n = 13 small and 12 large steps, t = 10.6, P = 2.3 × 10–10; P1: n = 13 small and 9 large steps, t = 6.4, P = 3.1 × 10–6; P6: n = 13 small and 10 large steps, t = 4.8, P = 9.9 × 10–6; P3: n = 32 small and 24 large steps, t = 23.0, P < 10–15; P4: n = 10 small and 7 large steps, t = 12.1, P < 10–15; P2: n = 15 small and 15 large steps, t = 3.0, P = 0.006; P8: n = 15 small and 10 large steps, t = 5.3, P = 2.0 × 10–5; P5: n = 12 small and 8 large steps, t = 3.0, P = 0.0073). c, Body schemes report the lower extremity motor score (LEMS) measured before EESREHAB and after EESREHAB. Relative changes in lower limb motor and sensory scores after EESREHAB are also reported. d, WISCIII score measured before EESREHAB and after EESREHAB with EESOFF and EESON. e, Bar plots reporting the outcome of the clinically standardized six-minute walk test and 10-meters walk test before EESREHAB, and after EESREHAB with EESOFF and EESON. Participants used a front-wheel walker for stability, but did not receive any external assistance. Two participants could not walk without assistance at the hip or at least 10 to 30% of body weight support for stability. Consequently, they did not perform these evaluations. Line plots report the evolution of the required amount of body weight support (% of bodyweight) necessary to enable each participant to walk independently over the period of EESREHAB. Bar plots show the mean with individual data points overlaid. Error bars show the standard error of the mean. *, P < 0.05; **, P < 0.01; ***, P < 0.001. Source data

Extended Data Fig. 2. Mouse model replicating the key technological and therapeutic features of EESREHAB in humans.

a, Three-dimensional visualization and quantification of the contusion SCI. Photographs show multiple 3D views of a contused spinal cord that has been cleared using uDISCO. The extent of spinal cord damage was quantified from coronal sections immunolabeled against glial fibrillary acidic protein (GFAP). The relative amount of spared tissues was quantified at the lesion epicentre and at regular intervals from the lesion epicentre in the rostral and caudal directions, as reported in the line plots. Sparing was consistent across groups (mean = 10.6%, n = 5 mice per group; independent samples two-tailed t-test on lesion epicenter, t = 0.5, P = 0.61). The impact of the contusion SCI on descending pathways compared to an uninjured spinal cord was visualized in whole brain-spinal cord preparations cleared with CLARITY. Projections from neurons located in the leg region of the motor cortex (AAV5-CAG-COMET-GFP) and from GluT2ON neurons located in the reticular formation (AAV5-CAGDIO-COMET-tdTomato) were labeled with virus infusions, as indicated in the photograph. Insets show the complete interruption of corticospinal fibers and partial preservation of reticulospinal fibers. b, Adaptations of the technological features of EESREHAB in mice, including a new robotic body weight support interface with 1g accuracy and a pair of electrodes attached to L2 and S1 spinal segments to deliver EES. Kinematic recordings of whole-body movements were mapped onto a 3D model of the mouse including bones and skin contours to generate realistic visualization of leg movements. The line shows the trajectory of the lower limb endpoint (toe) while circles indicate the maximum step height. Representative leg movements are shown concomitantly to the electromyographic (EMG) activity of the tibialis anterior, recorded from chronically implanted bipolar electrodes into the muscle. At one week, EES was insufficient to reactivate the spinal cord to a level that enabled walking in mice. Consequently, a small bolus of agonists to 5H1A/7 (8-OH-DPAT, 0.05-0.2 mg/kg) and 5HT2A receptors (quipazine 0.2-0.3 mg/kg) was administered to augment the response of the spinal cord. The combination of EES and 5-HT agonists enabled all the mice to walk overground as early as one week after the contusion SCI. 5-HT agonists were used to enable sustained walking during EESREHAB, but the amount of drug was progressively decreased over the course of the recovery, and eventually suppressed. Final testing was performed without 5-HT agonists. c, We noticed that compared to rats and humans, EES was not as effective to enable walking in mice with SCI. Inspection of muscle responses to EES revealed that the stimulation recruited motor nerves (efferents) at relatively low amplitudes. This off-target stimulation of the ventral roots is due to the relatively small size of the mouse spinal cord. We previously showed that high-frequency bursts are effective to maximize the activation of large-diameter afferent fibers targeted by EES. To adapt this stimulation protocol to mice, we elaborated a computational model of the mouse spinal cord. The model includes a spiking neural network model of muscle spindle feedback circuits for a pair of antagonistic muscles. Response in motor neurons are shown for two amplitudes of EES that correspond to the recruitment of 20% and 60% of the entire afferent population. Simulations suggested that, compared to conventional protocols (40 Hz), high-frequency burst stimulation (carrier frequency 600 Hz, modulating frequency 30 Hz) leads to a summation of synaptic events into motor neurons, which augment the probability to activate muscles via the recruitment of afferents compared to efferents. To confirm these findings experimentally, we compared conventional stimulation protocols to high-frequency burst stimulation during stepping on a treadmill in mice with contusion SCI. Reconstructed leg movements including foot trajectories are shown for each experimental condition at two amplitudes of EES. High frequency burst stimulation increased the therapeutic window through which stimulation could be delivered. Concretely, step height scaled linearly up to 1.5× motor threshold using burst stimulation protocols, and stimulation could be delivered effectively up to 2.5× motor threshold. In contrast, conventional stimulation at 40 Hz facilitated stepping around motor threshold, but additional increase in amplitude that would be necessary to promote robust stepping led to tonic activation of leg muscle activation and thus cessation of stepping. d, To model the volitional control of stepping observed in humans during EESON, we manipulated cortical activity with optogenetics. We expressed channelrhodopsin in the neurons of the leg motor cortex using targeted injections of AAV5-hSyn-ChR2 (n = 4 mice). Mice with contusion SCI were recorded during stepping on a treadmill with EESON (no 5-HT agonists). Photostimulation of the motor cortex induced a significant increase of the step height that scaled up with laser intensity, as shown in the bar plots (n = 4 mice per group; statistics indicate Tukey HSD tests following one-way ANOVA, P = 0.0003, P = 0.025, respectively). e, Leg movements and muscle activity during overground walking without any intervention or support recorded at four weeks after a contusion SCI in a mouse that did not undergo EESREHAB and in a mouse that underwent EESREHAB. Locomotor performance was quantified using principal component analysis applied to 80 gait parameters calculated from kinematic recordings (Supplementary Table 2). In this denoised space, each dot represents a gait cycle (n = 20 per mouse, n = 9 mice per group). Larger dots represent the mean of each experimental group. The first principal component (PC1) distinguished gaits from mice with SCI that did not undergo EESREHAB from mice that underwent EESREHAB. Locomotor performances were thus quantified as the scores on PC1 (reported in Fig. 1g). Analysis of factor loadings on PC1 revealed that the percentage of paw dragging, the extent of whole-limb oscillation (virtual limb connecting the hip to the toe) and step height were the parameters that showed the highest correlation with PC1. Bars report the mean values of these gait parameters (n = 9 mice per group; statistics indicate Tukey HSD tests following oneway ANOVA, P = 2.1 × 10–9, P = 7.4 × 10–6, P = 0.0055, respectively). Data from mice that underwent EESREHAB are shown during EESOFF and EESON. f, Photographs show whole mouse spinal cords before and after processing with iDISCO+,, during which the spinal cords underwent immunolabeling of cFos followed by clearing. The spinal cords were imaged with the mesoSPIM lightsheet microscope. Representative microscopy images show a raw coronal optical slice of the cFos labelling in the spinal cord and after application of automated 3D nuclear spot detection. Images were then reconstructed to visualize the entire lumbar spinal cord. g, 3D cFos quantifications were confirmed using immunohistochemistry and labelling for cFos on coronal sections of the lumbar spinal cord, as illustrated in the representative photographs of spinal cord sections from mice with SCI and mice that followed EESREHAB. The bar plot reports the mean number of cFos labeled cells per section in the spinal grey matter across the whole section, in Lamina I–III (dorsal), in Lamina IV–VI (intermediate), in Lamina VII–IX (ventral) and in Lamina X (central canal) (n = 4 mice per group; independent samples two-tailed t-test, t = –2.7, P = 0.042; t = 0.60, P = 0.57; t = 3.7, P = 0.010; t = 1.2, P = 0.27; t = 1.0, P = 0.35, respectively). Bar plots show the mean with individual data points overlaid. Error bars show the standard error of the mean. *, P < 0.05; **, P < 0.01; ***, P < 0.001. Source data

Extended Data Fig. 3. Single-nucleus RNA-sequencing of the mouse lumbar spinal cord.

a–f, Quality control statistics for 82,093 nuclei from the mouse lumbar spinal cord. a, Number of unique molecular identifiers (UMIs) per nucleus. Inset text shows the median number of genes detected. b, Number of genes detected per nucleus. Inset text shows the median number of genes detected. c, Proportion of mitochondrial counts per nucleus. Inset text shows the median proportion of mitochondrial counts. d, As in a, but shown for individual libraries from each experimental condition. e, As in b, but shown for individual libraries from each experimental condition. f, As in c, but shown for individual libraries from each experimental condition. d–f, n = 3 mice per condition. g, UMAP visualization of 82,093 nuclei colored by experimental condition, revealing consistent recovery of all major transcriptional states in each experimental condition. h, Proportions of each major cell type of the lumbar spinal cord recovered in individual libraries from each experimental condition. i, UMAP visualizations of nuclei from each experimental condition (n = 3 mice per condition). j, Number of UMIs quantified per nucleus in each major cell type of the mouse spinal cord. k, Number of genes detected per nucleus in each major cell type of the mouse spinal cord. l, Proportion of mitochondrial counts per nucleus in each major cell type of the mouse spinal cord. j–l, n = 3 mice per condition. m, Expression of key marker genes for the six major cell types of the lumbar spinal cord. Box plots show the median (horizontal line), interquartile range (hinges) and smallest and largest values no more than 1.5 times the interquartile range (whiskers), and the error bars show the standard deviation.

Extended Data Fig. 4. A single-nucleus atlas of neuronal subpopulations in the mouse lumbar spinal cord.

a, Clustering tree 52 of neurons in the mouse lumbar spinal cord, revealing a hierarchical taxonomy of transcriptionally defined neuron subpopulations over six clustering resolutions. b, Dot plot showing expression of a single marker gene per cell type for the 37 transcriptionally defined neuron subpopulations of the mouse lumbar spinal cord. c, Expression of key neurotransmitters for excitatory (Slc17a6) and inhibitory (Slc31a1) genes demonstrating a logical topography of neurons in the spinal cord. d, UMAP visualization of 20,990 neuronal nuclei colored by experimental condition, revealing consistent recovery of all major neuron transcriptional states in each experimental condition. e, Number of nuclei from each transcriptionally defined neuron subpopulation identified across all experimental conditions. f, Proportions of each neuronal subpopulation of the lumbar spinal cord recovered in each of eight experimental conditions. g, UMAP visualizations of neurons from each experimental condition (n = 3 mice per condition).

Extended Data Fig. 5. A spatial transcriptomic atlas of the mouse lumbar spinal cord.

a–c, Quality control statistics for 64 sections from the lumbar spinal cord profiled by spatial transcriptomics (a, mean number of UMIs per barcode; b, mean number of genes detected per barcode; c, mean proportion of mitochondrial counts per barcode). Shaded areas show the standard deviation. Sections are colored by the slide on which each section was sequenced. Sections highlighted in red failed quality control and were removed prior to further analysis. d–f, Overall distributions of quality control statistics aggregated over all sections (d, number of UMIs per barcode; e, number of genes detected per barcode; f, proportion of mitochondrial counts per barcode), shown separately for barcodes that passed or failed quality control, respectively. g, Relationship between the number of UMIs and number of genes detected per spatial barcode, with barcodes that failed quality control highlighted. h, Schematic overview of the registration procedure used to establish a cohesive atlas of spatial transcriptomes, with all 61 spinal cords registered to a common coordinate system. i, j, Quality control statistics visualized within the common coordinate system of the mouse lumbar spinal cord (i, number of UMIs per barcode; j, number of genes detected per barcode). k–l, Quality control statistics shown separately for each region of the lumbar spinal cord identified by registration to the Allen Brain Atlas (k, number of UMIs per barcode; l, number of genes detected per barcode; n = 61 sections from 12 mice). m, Proportion of barcodes from each spinal cord region recovered in individual replicates from each of four experimental conditions. n, Visualization of barcodes from each of four experimental conditions within the common coordinate system of the mouse lumbar spinal cord. Box plots show the median (horizontal line), interquartile range (hinges) and smallest and largest values no more than 1.5 times the interquartile range (whiskers), and the error bars show the standard deviation.

Extended Data Fig. 6. Integration of single-nucleus and spatial transcriptomes.

a, Spatial expression of key marker genes for canonical neuronal subpopulations of the mouse lumbar spinal cord. b, Deconvolution of spatial barcodes using snRNA-seq data to establish a spatial atlas of neuronal subpopulations. Panels show the score assigned by robust cell type decomposition (RCTD) to each of the 36 neuronal subpopulations for each barcode in the spatial transcriptomics dataset. c, Mean RCTD scores assigned for each of the 36 neuronal subpopulations to barcodes within five regions of the spinal cord.

Extended Data Fig. 7. RNAscope ISH confirms the spatial location and neurotransmitter phenotype of five neuron subtypes.

Spatial atlases of 5 neuronal subpopulations (Vsx2, Maf, Tac2, Chat and Pkd1l2) after deconvolution of spatial barcodes using snRNA-seq data and corresponding RNAScope ISH labelling for neuron subtype marker and neurotransmitter phenotypes Slc17a6 (Glut) and Slc32a1 (GABA).

Extended Data Fig. 8. Cell type prioritization to identify recovery-organizing neurons.

a, Expression of the immediate early gene cFos in neuronal nuclei from each of eight experimental conditions. b–g, Robustness of cell type prioritizations across six comparisons capturing the key therapeutic features of EESREHAB. Top, cell type prioritizations are shown within a clustering tree of spinal cord neurons declined in five different clustering resolutions, demonstrating the robustness of these prioritizations to the resolution at which transcriptionally defined neuronal subpopulations are defined. The five cell types with the highest AUCs in each comparison are annotated. Bottom, cell type prioritizations over increasingly granular clustering resolutions are visualized on a progression of UMAPs, with neuronal subpopulations colored by the strength of the perturbation response, as inferred by Augur. h, Volcano plots showing the strength of the perturbation response (x-axis) and its statistical significance (permutation test, y-axis) for six direct comparisons between two experimental conditions each capturing the key therapeutic features of EESREHAB.

Extended Data Fig. 9. Differential expression analysis of SCVsx2::Hoxa10 neurons.

a, Enrichment of immediate early genes among acutely perturbed SCVsx2::Hoxa10 neurons exposed to immediate perturbations (e.g., in response to walking with EES), as compared to chronically perturbed SCVsx2::Hoxa10 neurons (e.g., following EESREHAB). b, Top-ranked genes upregulated in chronically perturbed SCVsx2::Hoxa10 neurons, as compared to immediately perturbed neurons. c, Top-ranked GO pathways enriched in chronically perturbed SCVsx2::Hoxa10 neurons, as compared to immediately perturbed neurons.

Extended Data Fig. 10. Spatial prioritization of perturbation responses with Magellan.

a, Schematic overview of Magellan. b–g, Validation of Magellan in simulated data. Six different patterns of perturbation were simulated. Perturbation patterns were selected in order to illustrate key desiderata for a robust spatial prioritization method, including the ability to detect sharp borders of perturbation-responsiveness, smooth gradients of perturbation-responsiveness, and multifaceted perturbation responses of differing intensities. Top left, the ground-truth pattern of perturbation responsiveness in simulated spatial transcriptomics data. Top right, spatial prioritizations assigned by Magellan, with default parameters. Bottom left, robustness of spatial prioritizations assigned by Magellan when varying the number of spatial nearest-neighbors, k, used to assign a perturbation score to each barcode (default, k = 20). The parameter k controls the spatial resolution at which Magellan circumscribes the perturbation response. Bottom right, reproducibility of perturbation scores when varying the number of times three-fold cross-validation is repeated for each barcode (default, 50 times). h, Robustness of spatial prioritization to the number of spatial nearest-neighbors, k, used to assign a perturbation score to each barcode, shown for four direct comparisons between two experimental conditions each capturing key therapeutic features of EESREHAB.

Extended Data Fig. 11. Spatial embedding of single-nucleus transcriptomes.

a, Comparison of Magellan spatial prioritizations, left, with cell type prioritization scores assigned by Augur within the snRNA-seq dataset and embedded within the common coordinate system of the lumbar spinal cord by Tangram, right. b, Correlation between Magellan spatial prioritizations and Augur cell type prioritization scores assigned by Augur within the snRNA-seq dataset and assigned to matching spatial coordinates by Tangram. c, Comparison of Augur cell type prioritizations to Magellan spatial prioritizations for barcodes aligned with single-nucleus transcriptomes of each neuronal subtype. Left, Augur cell type prioritizations, as shown in Fig. 3b. Right, mean AUCs assigned by Magellan to spatial barcodes aligned with single-nucleus transcriptomes of each neuronal subtype in the snRNA-seq dataset. Source data

Extended Data Fig. 12. The connectome and projectome of SCVsx2::Hoxa10 neurons.

a, Schematic detailing the experimental set-up to study synaptic inputs to the SCVsx2::Hoxa10 neurons. Monosynaptically restricted EnvA pseudotyped and G-protein deleted rabies expressing mCherry revealed reticulospinal neurons in the vGi and PVON neurons in the DRG give direct inputs to SCVsx2::Hoxa10 neurons. b, Schematic detailing the experimental set-up to study synaptic inputs to the SCVsx2::Hoxa10 neurons. PVON afferent fibres are labeled in the transgenic PVCre::Advillin:tdTomato mice. Injections targeted to the vGi labeled descending reticulospinal projections. Vsx2ON neurons were labeled with RNAScope ISH. c, Schematic detailing the experimental set-up to study outputs from SCVsx2::Hoxa10 neurons. SCVsx2::Hoxa10 neurons were traced with AAVDj-hSyn-Flex-mGFPfp-2A-Synaptophysin-mRuby that labels the axons and pre-synaptic terminals. Synaptic appositions from SCVsx2::Hoxa10 neurons on to glutamatergic, GABAergic and motor neurons, labeled using RNAScope ISH were quantified. Bar graphs represent the mean with individual data points overlaid. Error bars reflect the standard error of the mean. Source data

Extended Data Fig. 13. Anatomical and functional features of SCVsx2::Hoxa10 neurons.

a, Optimization of the timing to identify supraspinal neurons connected to SCVsx2::Hoxa10 neurons using rabies viruses. Tracing was conducted using infusions of Cre-dependent G-deleted EnvA rabies or Cre-dependent pseudorabies (PRV) into the lumbar spinal cord of Vsx2Cre mice. For pseudorabies experiments, tissue was harvested at 2, 2.5, 3, 3.5 and 4 days after the infusion. Sagittal brain sections were used for 3D reconstructions of brains in Neurolucida, shown here for key timepoints. Quantifications report the number of neurons identified in each region of the brain and brainstem after rabies infusion, and at different timepoints after pseudorabies infusions. At 3.5 days days after pseudorabies injections, neurons were labeled in new brain regions compared to regions labeled with monosynaptically-restricted rabies, suggesting that secondorder neurons became transfected by the pseudorabies. We thus used a window of 3 days to study the connectome of SCVsx2::Hoxa10 neurons after SCI. b, Transverse spinal cord section from a Vsx2Cre mouse showing Cre-dependent expression of ChrimsonR in SCVsx2::Hoxa10 and the tract resulting from the insertion of one electrode shank. Inset shows the cell bodies of SCVsx2::Hoxa10 neurons in the vicinity of the tract. The schematic displays the 4-shank, 64-channel silicon probe used for single-unit recordings. The histogram reports the timing of spikes with respect to the onset of photostimulation for the SCVsx2::Hoxa10 optotagged single unit shown in Fig. 3d. The waveforms display spontaneous (gray) vs. optogenetically-evoked (red) spikes from an SCVsx2::Hoxa10 optotagged single unit. Data are shown as mean traces with standard error of the mean ribbons. The traces on the right were obtained from the same recording site during EES (top) and stimulation of the reticular formation (bottom). In these trials, two single-units responded to EES, whereas only one of these two units also responded to stimulation of the reticular formation (C56). The plot on the right reports the latencies of spikes from SCVsx2::Hoxa10 optotagged single units following EES across all trials (150 trials, top) and stimulation of the reticular stimulation (100 trials, bottom). Neurons that consistently responded to stimulations with short latencies (putatively monosynaptic) are highlighted in red. c, To evaluate the response of the spinal cord to EES, we measured muscle responses in the tibialis anterior when delivering a single pulse of EES. We tested uninjured mice, mice with chronic SCI, and mice that had undergone EESREHAB. Mice with chronic SCI exhibited abnormal long-latency responses (range: 10 to 20 ms), which are highlighted within the grey area and quantified in the bar plots as integral of root mean square (RMS) (n = 6 mice per group; statistics indicate Tukey HSD tests of the key comparison following one-way ANOVA, P = 0.0001). d, The evaluations reported in c for mice were also conducted in the cohort of human participants. Three of the participants (DM002, GO004, HT008) showed abnormal long-latency responses (range: 50 to 100 ms) to single-pulse of EES before EESREHAB. EESREHAB nearly completely abolished these responses. The bar plot reports the normalised amplitude of these responses before and after EESREHAB (n = 3, trials > 5/session; mixed-effects model, t = –6.40, P = 2.9 × 10–8). e, f, The role of SCVsx2::Hoxa10 neurons in the reorganization of muscle responses to EES was evaluated using Cre-dependent expression of Gi, e, or Gq DREADDs in SCVsx2::Hoxa10 neurons, f. The timelines summarize the timing of the various experimental procedures. The photographs illustrate the expression of DREADD receptors in SCVsx2::Hoxa10 neurons. Long-latency muscle responses to single-pulse of EES were quantified before and 30 min after CNO injections that either silenced (Gi) or activated (Gq) SCVsx2::Hoxa10 neurons. Bar plots report the integral of the RMS of long-latency muscle responses for each each experimental condition (n = 5 mice per group (Gi); paired samples two-tailed t-test, t = 2.4, P = 0.046; n = 5 mice per group (Gq); paired samples two-tailed t-test, t = 2.8, P = 0.047). g, Three complementary strategies were used to evaluate the role of SCVsx2::Hoxa10 neurons during basic unskilled walking in uninjured mice: targeted injections of AAV5-flex-hm4di (Gi) or AAV5-flex-Jaws in the lumbar spinal cord of Vsx2Cre mice to evaluate the short-term or immediate impact of silencing SCVsx2::Hoxa10 neurons, and targeted injection of AAV-flex-DTR to evaluate the long-term impact of the complete ablation of SCVsx2::Hoxa10 neurons. Locomotor performance was evaluated during overground walking using the procedures detailed in Extended Data Fig. 2e (n > 20 gait cycles per mouse, n = 5 mice per group). The bar plot reports locomotor performance, quantified as average scores on PC1 (n = 5 mice per group except for DTR with n = 8 mice per group; paired samples two-tailed t-test, Gi: t = 1.4; P = 0.23; LED: t = 1.5, P = 0.21; DTR: t = 1.1, P = 0.33), and the number of SCVsx2::Hoxa10 neurons per tissue section (SCVsx2::Hoxa10 ablation, n = 7 mice; no ablation, n = 4 mice; independent samples two-tailed t-test, t = 12.3, P = 6.2 × 10–7). Bar plots show the mean with individual data points overlaid. Error bars show the standard error of the mean. *, P < 0.05; **, P < 0.01; ***, P < 0.001. Source data

Extended Data Fig. 14. Gait analysis of uninjured mice after diphtheria toxin ablation of SCVsx2::Hoxa10 neurons.

a, Workflow of kinematic gait analysis. Step 1: Locomotor performances of uninjured mice were evaluated during quadrupedal walking on a treadmill. Bilateral leg kinematics were captured with twelve infrared cameras that tracked reflective markers attached to the crest, hip, knee, ankle joints and distal toes. Step 2: The limbs were modelled as an interconnected chain of segments and a total of 78 gait parameters were calculated from the recordings. All gait parameters are reported in Supplementary Table 2. Step 3–6: To evaluate differences between experimental conditions, as well as to identify the most relevant parameters to account for these differences, we implemented a multistep multifactorial analysis based on principal component analysis, which we described in detail previously,. b, Multifactorial analysis including the experimental groups of uninjured mice before and after diphtheria toxin ablation. Step 3a: In the principal component analysis plot, each dot represents a gait cycle (n = 20 per mouse, n = 8 mice per group). Larger dots represent the mean of each experimental group. The first principal component (PC1) usually distinguishes the largest differences in gait. In this case, only minimal non-significant differences were detected. The barplots quantify locomotor performances as the scores on PC1 (n = 8; paired samples two-tailed t-test; t = 0.8874; P = 0.4043). Step 4a: Analysis of factor loadings on PC1 reveal potential differences in 3 functional gait clusters (Step 5a). c, As in b, but including an additional group of mice that underwent hemisection injuries. Only the gait data from ipsilesional hindlimb are included in the analysis. Please note the obvious differences on PC1 and the higher correlation of gait parameters onto PC1 (n = 8; statistics indicate Tukey HSD tests following repeated measures one-way analysis of variance (ANOVA); P = 1.7 x 10−8). d, The barplots quantify key gait parameters found in the gait analysis in b-c. Bars report n = 8 mice per group, statistics indicate Tukey HSD tests following repeated measures ANOVA, P = 8.9 x 10−11, P = 2.7 x 10−6, P = 1.2 x 10−8, P = 2.5 x 10−13, P = 2.4 x 10−8, P = 6.4 x 10−4, P = 9.8 x 10−8 and P = 6.9 x 10−7, respectively. Bar graphs represent the mean with individual data points overlaid. Error bars reflect the standard error of the mean.  *, **, *** indicate P < 0.05, 0.01, and 0.001, respectively. Source data

Extended Data Fig. 15. SCVsx2::Hoxa10 neurons are necessary and sufficient for the recovery of walking after paralysis.

a, Exploded diagram of a new wireless optoelectronic system that integrates on the same implant red-shifted microLEDs to deliver deeply-penetrating photostimulation and electrodes to deliver EES. The zoom shows the microLEDs and electrodes for EES. b, Optogenetic silencing of SCVsx2::Hoxa10 neurons in mice that underwent four weeks of EESREHAB. Experiments and data are shown using the same conventions as in Extended Data Figs. 2 and 139 (n > 20 per mouse, n = 4 mice per group). Bars report n = 4 mice per group; statistics indicate Tukey HSD tests following repeated measures ANOVA, P < 10–15, P = 0.003 and P = 0.002, respectively. c, As in b, but during acute silencing of SCVsx2::Hoxa10 neurons with Gi in mice that underwent four weeks of EESREHAB (n = 4 mice per group; statistics indicate Tukey HSD tests following one-way ANOVA, P = 7.2 × 10–10, P = 2.1 × 10–5, P = 0.0006 respectively). d, As in c, but during acute activation of SCVsx2::Hoxa10 neurons with Gq in mice tested at four weeks post-SCI (no EESREHAB) with EESOFF and EESON (n = 4 mice per group; statistics indicate Tukey HSD tests following one-way ANOVA, P = 3.1 × 10–5, P = 4.4 × 10–7, P = 3.1 × 10–5, P = 0.0069 respectively). e, As in d, but following the chronic silencing of SCVsx2::Hoxa10 neurons (Gi, CNO in drinking water) in mice that underwent four weeks of EESREHAB, chronic activation of SCVsx2::Hoxa10 neurons (Gq, CNO in drinking water) in mice after SCI, and chronic activation of SCVsx2::Hoxa10 neurons (Gq, CNO in drinking water) in mice that underwent four weeks of rehabilitation without EES (n = 5 mice per group; statistics indicate Tukey HSD tests following one-way ANOVA, P = 0.0007, P = 1.6 × 10–5, P = 0.004, P = 0.011, respectively). f, Chronic silencing of SCVsx2::Hoxa10 neurons during EESREHAB altered the connectome and projectome of SCVsx2::Hoxa10 neurons compared to mice that underwent EESREHAB. 3D reconstructions show the labeled neurons in the reticular formation 3 days following the infusion of Cre-dependent PRV Ba2017 in the lumbar spinal cord of Vsx2Cre mice. The density of vGluT1 synapses apposing SCVsx2::Hoxa10 neurons was quantified to evaluate the synaptic projection from large-diameter afferent fibers onto SCVsx2::Hoxa10 neurons. The projectome of SCVsx2::Hoxa10 neurons was vizualized using infusions of AAV-flex-tdTomato in the lumbar spinal cord. The bar plots show the average number of vGluT1 synapses apposing SCVsx2::Hoxa10 neurons (n = 6 mice per group; independent samples two-tailed t-test, t = –3.1, P = 0.018), the number of neurons in the reticular formation (n = 4 and 3 mice per group, respectively; independent samples two-tailed t-test, t = –4.3, P = 0.0076), and the density of projections from SCVsx2::Hoxa10 neurons in the ventral horn (n = 4 mice per group; independent samples two-tailed t-test, t = –9.2, P = 0.0003). g, The impact of the chronic ablation of SCVsx2::Hoxa10 neurons on the natural recovery of mice with thoracic lateral hemisection SCI was evaluated as in Extended Data Fig. 13g, and summarized in the timeline of experiments. Photographs show coronal sections of the hemisected spinal cord at the lesion epicenter, while CLARITY-optimized light sheet microscopy illustrates the interruption of the reticulospinal tract on the hemisected side. Locomotor performances were evaluated as detailed in the other panels. Bar plots, n = 5 mice per group; independent samples two-tailed t-test, t = 6.6, P = 0.0002; t = –3.4, P = 0.09; t = –1.9, P = 0.093, respectively). h, Photographs show coronal sections of the spinal cord with Vsx2 RNAScope, confirming the reduction in the number of SCVsx2::Hoxa10 neurons in the lumbar spinal cord of Vsx2Cre mice after diphtheria toxin injections. Bar plots report the mean number of Vsx2 labeled neurons per section in the spinal grey matter (SCI n = 3, SCI with SCVsx2::Hoxa10 ablation n = 4 mice; independent samples two-tailed t-test, t = 9.5, P = 0.0001). Photographs below show examples of projections from reticulospinal neurons in the lumbar spinal cord in both groups of mice. The plot reports the mean density of reticulospinal projections across the dorsoventral extent of the spinal cord for mice with SCI (n = 3) and mice with SCI and ablation of SCVsx2::Hoxa10 neurons (n = 4). Ribbons show the standard deviation (two-tailed Wilcoxon rank-sum test, W = 2008248, P < 10–15). Bar plots show the mean with individual data points overlaid. Error bars show the standard error of the mean. *, P < 0.05; **, P < 0.01; ***, P < 0.001. Source data