Abstract
Propriospinal neurons play crucial roles in recovery of function after spinal cord injury. The V2a class in particular contributes significantly to recovery of locomotor and respiratory function after injury in animal models. However, V2a are a diverse class with different properties and gene expression patterns, suggesting that V2a subtypes might have disparate responses to injury. Here, we used single-nucleus RNA sequencing of enriched cervical V2a neurons from adult female mice to identify specific markers of 10 V2a subtypes. We found differences in medial-lateral location as well as rostral-caudal distribution of V2a belonging to specific clusters in the adult cord. Although adult V2a can be identified as belonging to N- and Z-group divisions similarly to neonatal V2a, there is a not a 1:1 correlation between most neonatal and adult clusters. We assessed changes in gene expression in each V2a subtype 5 d following a C2 hemisection injury. Remarkably, three subtypes of V2a neurons were found to be nearly absent ipsilateral and contralateral to injury, demonstrating that some V2a neurons are more vulnerable to changes in identity or cell loss than others. Further, the remaining V2a subtypes exhibited distinctive transcriptional alterations indicating a potential rewiring of cellular interactions within the spinal cord microenvironment. Gene expression changes common to most V2a subtypes included an upregulation of RNA splicing factors. This study reveals the diverse responses of propriospinal neurons to injury and provides a foundation for understanding changes in propriospinal neurons that may lead to adaptive (or maladaptive) changes in circuit function following injury.
Significance Statement
This study is the first to describe markers of adult V2a propriospinal neuron subtypes and their relationship to developmental V2a neuron subtypes. In addition, we characterize the diverse responses of different V2a neurons to cervical spinal cord injury, including the loss of entire V2a subtypes. These analyses nominate new cellular targets for altering propriospinal neuron function, plasticity, and regeneration following traumatic injury.
Introduction
Spinal cord injuries disrupt communication between the brain and spinal cord leading to long-term disability, including muscle paralysis, loss of sensation, impaired breathing, immune dysfunction, and impaired regulation of autonomic functions (Berlowitz et al., 2016; Eldahan and Rabchevsky, 2018; Badhiwala et al., 2019; Rodgers et al., 2022). A better understanding of which cells are important for recovery of function and how they respond to injury is necessary to develop improved therapies to treat spinal cord injury. Propriospinal neurons play a critical role in recovery of motor function after spinal cord injury (Bareyre et al., 2004; Courtine et al., 2008; Filli and Schwab, 2015; Jensen et al., 2020; Zavvarian et al., 2020). Recent studies have demonstrated that the V2a class of excitatory propriospinal neuron is particularly important for recovery of function after spinal cord injury (Huang et al., 2022; Kathe et al., 2022; Van Steenbergen et al., 2023; Jensen et al., 2024). V2a are excitatory, ipsilaterally projecting neurons in the intermediate laminae of the spinal cord that play important roles in many motor behaviors including locomotion, respiration, and goal directed movements (Crone et al., 2008, 2009, 2012; Azim et al., 2014; Romer et al., 2017; Jensen et al., 2019). V2a neurons undergo changes in connectivity in both locomotor and respiratory circuits that likely contribute to recovery of function (Zholudeva et al., 2017; Van Steenbergen et al., 2023). Moreover, activating V2a neurons can restore activity to a previously paralyzed diaphragm in a mouse model of cervical spinal cord injury (Jensen et al., 2024). V2a neurons are also a target of epidural stimulation and are critical for recovery of locomotor ability after injury (Kathe et al., 2022; Squair et al., 2023). Thus, a better understanding of the changes that occur in V2a neurons after injury could help identify molecules and pathways important for driving recovery.
V2a are a heterogeneous group of neurons with diverse gene expression profiles as well as electrophysiological and functional properties (Dougherty and Kiehn, 2010; Zhong et al., 2010; Hayashi et al., 2018). V2a neurons from neonatal spinal cord can be divided into at least 11 subtypes based on gene expression profiling that are derived from 2 broad groups during development: Z-group neurons that are early born, laterally located, and mostly long-range projecting and N-group neurons that are late-born, medially located, and more likely to project locally (Delile et al., 2019; Osseward et al., 2021; Sagner et al., 2021). However, it is currently unknown whether V2a subtypes (or even the N- and Z- groups) can be identified in adult spinal cords, as many markers expressed during development are downregulated in adult cord.
The C2 hemisection injury (C2Hx) has long been used as an experimental model to study neuroplasticity in rodents and other species (Porter, 1895; Lewis and Brookhart, 1951; Goshgarian, 2003, 2009; Hoh et al., 2013; Warren and Alilain, 2014). This model causes selective damage to the ascending and descending pathways on one side of the spinal cord, resulting in paralysis of the hemidiaphragm ipsilateral to the injury, but sparing the animals ability to walk, breath, and regulate autonomic functions. Moreover, since some animals spontaneously recover diaphragm function between 4 and 14 d after C2Hx injury (Dougherty et al., 2012; Mantilla et al., 2014; Komnenov et al., 2016; Jensen et al., 2024), this model can be used to identify genes and pathways whose expression is altered during this time frame that potentially contribute to recovery.
We hypothesized that V2a neurons in the adult spinal cord comprised multiple distinct subtypes that differentially respond and adapt to spinal cord injury. To test this hypothesis, we profiled V2a neurons using single-cell genomics, prior to and following a C2 hemisection spinal cord injury. We found diverse transcriptional responses of different V2a neurons to cervical spinal cord injury, including the loss of entire V2a subtypes.
Materials and Methods
Mouse lines
All animal experiments were performed according to guidelines provided by the National Institutes of Health and approved by the Cincinnati Children’s Hospital Medical Center’s Animal Care and Use Committee. The following strains of mice were used: Chx10Cre (Azim et al., 2014; Bouvier et al., 2015; Romer et al., 2017; Jensen et al., 2019), Sun1/sfGFP reporter mice (ROSA26CAG-LSL-Sun1-sfGFP-Myc/CAG-LSL-Sun1-sfGFP-Myc; RRID:IMSR_JAX:021039), Lgr5tm1(cre/ERT2)Cle (RRID:IMSR_JAX:008875), and ROSA-tdTomato (ROSA26Sortm14(CAG-tdTomato)Hze, RRID:IMSR_JAX:007914). Mouse lines were maintained on a C57BL/6 background.
To fluorescently label nuclei of V2a neurons, V2a-Sun1GFP mice (Chx10Cre/+; ROSA26CAG-LSL-Sun1-sfGFP-Myc/+) were generated by crossing Chx10Cre/+ or Chx10Cre/Cre mice to homozygous Sun1/sfGFP reporter mice. The expression of Cre in Chx10 expressing cells removes a loxP flanked stop sequence which allows expression fusion protein containing Sun1 (a nuclear envelope protein), two copies of super folder GFP (GFP), and six copies of a Myc epitope tag (Mo et al., 2015). A total of 14 female mice were used for single-nuclei RNA sequencing: 6 uninjured (No Injury 1 and No Injury 2) and 8 following a C2 hemisection injury (Post Injury 1 and Post Injury 2) at 6–12 weeks of age. Although rare, we occasionally observe broad expression of reporter genes in Chx10Cre/+ animals, likely due to sporadic expression of Cre recombinase in the germline or early embryo, as reported for other Cre lines (Luo et al., 2020). We used PCR to detect abnormal recombination in tail DNA from V2a-Sun1GFP mice and excluded these mice from all experiments. Tissue used to count Necab1-positive V2a neurons came from Lgr5-eGFP-IRES-Cre-ERT2; ROSA-tdTomato mice; however, no tamoxifen was given. We did not count the eGFP+ (Lgr5) neurons.
Isolation of V2a nuclei from spinal cords of adult mice
Adult female V2a-Sun1GFP mice were anesthetized with pentobarbital (0.1 mg/g, i.p.) and decapitated. The cervical spinal cord was dissected and washed with 0.1 M RNase-free phosphate-buffered saline (PBS). The spinal cord was then snap-frozen in dry ice and stored at −80°C. The isolation of nuclei from the spinal cord was performed as described previously (Chamessian et al., 2018; Serafin et al., 2019, 2021). On experiment day, tissue was placed in an ice-cold homogenization buffer (HB). HB contains 250 mM sucrose; 25 mM KCL; 20 mM Tris-HCl, pH 8.0; 5 mM MgCl2; 1 tablet per 10 ml of Roche cOmplete Mini EDTA-free Protease Inhibitor Cocktail (Millipore Sigma, #4693159001), 40 U/ml RNasin Plus (Promega, Fisher Scientific, #N2611), and 1 μM DTT (Sigma). All the subsequent steps were carried out at 4°C. The cervical spinal cord was homogenized using 40 strokes of loose pestle, and 30 strokes of tight pestle. Then, 60 μl of 10% NP-40 was added to the homogenate followed by 16 more stokes of the tight pestle. Filtration and centrifugation were carried out as described previously (Chamessian et al., 2018; Serafin et al., 2019, 2021). The liberated nuclei were then suspended in 1 ml of HB.
As the cervical spinal cord has a lot of myelin, an additional myelin removal step was performed using myelin removal beads. The nuclei suspension was centrifuged for 5 min at 300 × g at 4°C. The supernatant was discarded, and the nuclei were resuspended in 1 ml of ice-cold blocking buffer (BB). BB contains 2.5 ml of 10× PBS, 2.5 ml of 10% bovine serum albumin, 125 μl of RNasin 40 U/µl RNasin Plus, and 20 ml of nuclease-free water. Then, 10 μl of Myelin Removal Beads II (Miltenyi Biotec, #130-096-731) was added and mixed by gently triturating. The solution was left for 15 min. Then, 1,000 μl of BB was then added and triturated. The nuclei/bead suspension was centrifuged for 5 min at 300 × g. The supernatant was gently removed, and the pellet was resuspended in 1,000 μl BB. The tube was then placed on a magnet for 15–20 min. Carefully the supernatant was removed without disturbing the myelin debris on the side of the tube. The supernatant was filtered through a 20 μm Sysmex filter into a FACS tube.
Two uninjured cords (No Injury 1) that were pooled together did not go through myelin removal and the rest of the spinal cords did undergo myelin removal. Both sides of the cord (ipsilateral and contralateral to injury) were collected and analyzed together.
Animal surgery
C2 hemisection was performed in adult mice (6–12 weeks of age). A laminectomy was performed using microscissors to expose the cervical spinal cord levels. Subsequently, a 30 G needle was used to lesion the left side of the spinal cord, extending from the midline to the lateral side (repeated five times to ensure a complete lesion). Muscles were sutured, and adhesive was applied to seal the skin incision. Immediately postsurgery, carprofen (5 mg/kg*BW) was administered subcutaneously as an analgesic. The mice were housed in 29°C incubator overnight and given nutritional gel and a water bottle. Subcutaneous injections of 1.0 ml saline were given twice daily for the initial 2 postoperative days. Additional, analgesia was provided if necessary. Animals were killed 5 d after the injury for tissue harvest and nuclei isolation.
Fluorescence-assisted nuclei sorting to enrich for V2a nuclei
The nuclei suspension was stained with Hoechst 33342 (Thermo Fisher Scintific; H3570). Fluorescence-assisted nuclei sorting (FANS) was performed at the Research Flow Cytometry Facility (RRID:SCR_022635) in the Division of Rheumatology at Cincinnati Children’s Hospital Medical Center. A BD FACS Aria II was used with a 70 µm nozzle to sort intact GFP+ nuclei into a well of round-bottom 96-well plate filled with 0.1 M RNase-free PBS containing 2% non-acetylated BSA (Sigma) and 0.5 U/μl RNasin Plus (Promega). GFP+/Hoechst−, GFP−/Hoechst+, and GFP−/Hoechst− control samples were used for determining gating parameters. The sorted nuclei were transferred to a 1.5 ml Eppendorf tube to be sent for sequencing.
10x Genomics library preparation and sequencing
The Single Cell Genomics Facility (RRID:SCR_022653) at Cincinnati Children’s Hospital Medical Center assessed the approximate concentration and viability of the nuclei. Two cohorts of uninjured (first containing two cords and then second containing four cords) and two replicates of injured (four cords each) were run as separate samples on the 10x Chromium instrument and constructed as separate libraries. The libraries were assessed for quality using a 2100 Bioanalyzer High Sensitivity DNA Assay (Agilent) and then sequenced on a NovaSeq 6000 sequencer (Illumina) using an SP flow cell. RNA sequencing was performed using Single Cell 3ʹ Gene Expression Assay v3 (10x Genomics) by Genomics Sequencing Facility (RRID:SCR_022630).
Bioinformatics analyses
Cell Ranger (10x Genomics) version 6.0.1 was used to demultiplex Illumina base call files (BCL) using mkfastq, align the reads to mouse reference genome mm10 (with intron), and quantification of gene expression using Cell Ranger count. Ambient RNA was removed using 0.5 contaminant cutoff in SOUPX. Further analysis was performed using Seurat v5.0.0 in R v4.4.0. A total of 2,968 nuclei were captured from the uninjured and 1,033 nuclei were captured from the injured samples. The number of nuclei captured, mean reads per count, and median genes per cell are given in detail in Table 1.
Summary of cervical spinal cord cohorts used for RNA sequencing
Uninjured cohorts were first analyzed to identify different clusters of V2a neurons. Initial quality control was done such that nuclei that expressed at least 200 genes and genes that were expressed in at least three nuclei were selected. The two uninjured cohorts raw counts were combined together using merge() function. The merged uninjured dataset was then normalized using the NormalizeData() function in Seurat where feature counts for each nucleus are divided by the total counts for that nucleus and multiplied by a scale factor of 10,000. The result is then transferred to the natural log using log1p. Merged dataset was subjected to FindVariableFeatures(), ScaleData(), and RunPCA() before integrating them using IntegrateLayers() [Reciprocal Principal Component Analysis (“RPCA”) method] and JoinLayers() functions from Seurat. We used the RPCA method as we expected there to be slight differences in cell types because of the differing methods of sorting the nuclei. This integrated dataset was used for principal component analysis, variable gene identification, Shared Nearest Neighbor (SNN) clustering analysis, and Uniform Manifold Approximation and Projection (UMAP). We used a resolution of 1 to identify 23 clusters using FindClusters(). Rbfox3, Snap25, Syp, and Snhg11 genes were used to identify neuronal nuclei; Aqp4, Atp1a2, Gja1, Slc1a2, Ctss, Itgam, Pptrc, Mbp, Mobp, Mog, and Plp1 genes were used to identify glial nuclei; Slc17a6 genes were used to identify excitatory nuclei; Gad1, Gad2, and Slca32a1 genes were used to identify inhibitory nuclei; and Vsx2, Lhx3, Lhx4, and Shox2 genes were used to identify V2a nuclei. Lmx1b, Calb1, Calb2, Tac1, Tac2, Npy1r, Grp, Grpr, Chat, and Nts were used to mark non-V2a nuclei. From a total of 23 transcriptionally distinct clusters, we removed 9 that contained glial markers, 1 that contained inhibitory neuron markers, and 3 that contained non-V2a neuronal markers. The remaining 10 clusters were subset from the original dataset and reclustered using a resolution of 2. At this stage, we removed low-quality nuclei that expressed <1,000 genes. Using res of 2, we identified and removed 2/21 clusters that did not express V2a marker genes. These two clusters were probably like V2a neurons in gene expression and thus initially clustered together with V2a nuclei. The remaining 19 clusters were labeled as V2a neurons and contained nuclei from each uninjured cohort of spinal cords. A more conservative resolution of 0.5 was used to identify V2a clusters consisting of 1,251 nuclei for the uninjured dataset.
For the injured cohorts, a similar method of analysis was performed as for the uninjured cohorts, except as noted below. Like the uninjured cohorts, we selected 1,020 nuclei that expressed at least 200 genes and genes that were expressed in at least three nuclei. Further, we selected high-quality nuclei by removing 220 nuclei that expressed <1,000 genes. Contrary to the uninjured cohort analysis, instead of RPCA, we used the CCA method for integration as we did not expect a difference in cell types between our injured samples because they were sorted similarly. A total of 800 injured V2a nuclei were retained after removing 233 nuclei.
We next used Seurat [IntegrateLayers (CCA method) and JoinLayers functions] to integrate the V2a neurons from each cohort of cords (Injured 1, Injured 2, Uninjured 1, Uninjured 2) into a single dataset containing 10 clusters of V2a. Vsx2 was detected in 7% of cells, with all clusters containing cells with V2a markers (Vsx2, Lhx3, Lhx4, Shox2). The integrated dataset (Uninjured/Injured cohorts) was used for principal component analysis, variable gene identification, SNN clustering analysis, and UMAP.
For comparing the uninjured versus injured cohorts within the integrated dataset, the FindTransferAnchors(), TransferData(), AddMetaData(), RunUMAP(), and MapQuery() functions were used to transfer the labels from the reference (Uninjured Dataset) to the query (Integrated Injured/Uninjured Dataset).
cellHarmony from AltAnalyze version 2.1.4 was used to compare the gene expression changes in cluster-cluster V2a population after injury using software defaults. The log fold cutoff threshold was set to 1.2 and the Pearson’s correlation was set to 0.3.
Cell prioritization to identify most perturbed cluster was performed using Augur (1.0.3). Clusters 7, 8, and 9 were not included for the analysis due to the low number of cells in these clusters after injury. We used the calculate_auc() function with min.cells = 20 to calculate the area under the curve (AUC) for each cluster.
Analysis of regeneration-associated genes (RAGs) was performed as follows. We calculated a regeneration module score for each cell based on the expression of RAGs (list from Saraswathy et al., 2024) using R package: UCell (Andreatta and Carmona, 2021). To assess the differences between conditions and clusters, we fitted a linear mixed-effects model with cluster and condition (injured, uninjured) as fixed effects, their interaction, and sample (different batches of sequencing) as a random effect. Type III ANOVA tested overall effects, and estimated marginal means (emmeans) with pairwise contrasts were used for post hoc comparisons within each cluster.
Analysis of neonatal V2a dataset
To analyze the published Neonatal V2a dataset (Hayashi et al., 2018), we followed the standard protocol of Seurat. We first subset cells to include features >100 and <7,000. We then normalized the data using NormalizeData(). Using the top 2,000 variable features, we scaled our data and ran principal component analysis and used FindNeighbors() and FindClusters() with resolution of 1.2 to cluster the dataset, yielding 11 clusters.
Immunohistochemistry
Cervical spinal cord sections of V2a-Sun1GFP mice were used for immunohistochemistry to validate our RNA sequencing findings. The process of tissue preparation has been previously described (Romer et al., 2017; Jensen et al., 2019). Briefly, adult mice were anesthetized using pentobarbital (0.1 mg/g, i.p.) and transcardially perfused using cold phosphate buffer (PB) followed by cold 4% paraformaldehyde (PFA) dissolved in PB. Cervical spinal cord tissue (C1-C8) was harvested, rinsed overnight in PBS, and cryoprotected overnight using 30% sucrose at 4°C. The tissues were mounted in OCT Tissue Freezing Medium and stored at −80°C. Then, 14 μm transverse sections of the cervical spinal cord were cut on a cryostat microtome and mounted onto SuperFrost Plus slides (Fisher). Every 16th section from C1-C8 (at least 45 hemisections/animal) was analyzed for each immunostaining experiment. Slides were washed with 1× PBS with 1 ml of 20% Tween-20 (PBT). A blocking step was performed using antibody dilution buffer (PBT with 10% heat-inactivated donkey serum and 0.3% Triton X-100). Primary antibodies were incubated for 2 h at room temperature followed by 10 min washes with PBT four times. This was followed by 1 h incubation of secondary antibodies and DAPI (1:60,000 in PBT) for nuclear stain. Finally, the slides were dipped in CuSO4 with 0.02% Tween 20 for 15–20 min to reduce lipofuscin autofluorescence and rinsed in PBT before mounting with ProLong Gold Antifade Reagent (Life Tech, #P36930; Schnell et al., 1999). The primary antibodies used were Chx10 (Sheep at 1:2,000, RRID:AB_302278), Nfib (Rabbit at 1:1,600, RRID:AB_1854424), Zfhx4 (Rabbit at 1:100, RRID:AB_1859013), Necab1 (Rabbit at 1:2,000, RRID:AB_1848014), GFP (Chicken at 1:2,000, RRID:AB_2307313), GFP (Rabbit at 1:5,000, RRID:AB_221570), and Cleaved Caspase-3 (Rabbit mAb at 1:200, RRID:AB_2070042). The secondaries used were anti-Rabbit Cy3 (1:500, RRID:AB_2307443), anti-Chicken Alexa Fluor 488 (1:500, RRID:AB_2340375), and anti-Rabbit Alexa Fluor 750 (1:500, RRID:AB_2924801), all made in donkey.
RNAScope for in situ RNA detection
In situ RNA detection was performed in adult female mice using the RNAScope Multiplex Fluorescent Kit v2 [Advanced Cell Diagnostics (ACD), #323100] using the manufacturer’s directions for fresh-frozen tissue. The probes used were Nfib (ACD, #586511-C2), Zfhx4 (ACD, #531081), Pde3a (ACD, #1157721-C1), Sema5a (ACD, #508091-C3), and Syt10 (ACD, #52257-C1) in conjunction with Opal 620 and Opal 690 dyes (Akoya Biosciences). Positive control probe (ACD, #320881) and negative control probe (ACD, #320871) were used as controls. Antibodies to GFP were used to label Sun1-GFP V2a nuclei due to weak endogenous GFP fluorescence following protease digestion of tissue. After the in situ hybridization steps, PBT was used to wash the slides (three times for 5 min per wash) followed by incubation with antibody dilution buffer for 1 h. Overnight incubation with primary antibodies (Rabbit anti-GFP at 1:500, RRID:AB_221570) was performed overnight at 4°C. The secondary antibody (anti-Rabbit 488 at 1:500, RRID:AB_2307443) was incubated for 1 h and rinsed with PBS (three times for 5 min per wash). Finally, DAPI was applied for 30 s. The slides were mounted with ProLong Gold Antifade Reagent and coverslipped.
Imaging
All confocal images were acquired using equipment maintained by Bio-Imaging and Analysis Facility (RRID:SCR_022628). For immunohistochemistry experiments, images of cervical (C1-C8) spinal cord segments were collected with a single optical plane tile scanning protocol to capture the whole cord using a 20× objective and the Galvano scanner of a Nikon A1R inverted LUNV confocal microscope. In addition, 60× z-stack images of individual or groups of cells were taken using the same microscope.
For in situ hybridization, 60× z-stack images of cervical spinal cord capturing most of the GFP (V2a) nuclei were taken using the Confocal Spinning Disk Microscope Yokogawa SoRa W1. The positive and negative probes were used to determine the optimal settings to take images across the range of expected expression levels.
Quantification
Counts/analyses were performed using Imaris software (version 10.0.0). For immunohistochemistry, the number of GFP+, Chx10+/GFP+, Nfib+/GFP+, Zfhx4+/GFP+, and Necab1+/Chx10+ cells containing detectable DAPI were counted using the spot feature from 20× tiled single plane images. The middle of the central canal was used as the origin for the reference frame with the positive y-axis toward the dorsal horn. The location of all the spots with the central canal as a reference was exported. To merge data from the left and right sides of the cord, all the x-axis values were converted to positive. The distance from the central canal to the lateral, dorsal, and ventral edges (at the level of the central canal) was measured to account for size differences between sections. We mapped the cells from all C1-C8 sections onto a single reference section by normalizing the coordinates to a reference section with distances between the central canal and the lateral, dorsal, and ventral edges of 650, 400, and 400 µm, respectively. These values were then imported into Matlab. The kde function was used to generate contour plots and surface plots showing the distribution of the neurons in the spinal cord.
For in situ hybridization, 60× images were used to count the number of spots of each probe touching or inside the SUN1GFP+ nuclei. Firstly, using NIS elements software (version 5.41.02), the 60× images were preprocessed using rolling ball background subtraction and segmented to create a binary for GFP nuclei using general analysis 3 (GA3). These images were then further processed in Imaris. The Surface feature was used to identify all the binaries. If there were errors with binary creation, the surfaces were created manually. The surfaces were then filtered such that only those with a volume greater than 200 µm were selected to eliminate small fragments of nuclei. Next, the spot feature with optimal threshold was used with x–y diameter of 0.25 to determine the spots for each probe (Kwon et al., 2017). The threshold was kept constant for all the sections of an experiment. Only the spots that were inside or touching the GFP surfaces were preselected. Using rtmatlab 9.6 in Imaris, the number of spots for each surface was calculated. The nucleus was said to express a probe only if the number of spots were greater than 3 for Opal620 and greater than 5 for Opal690.
Statistical analyses
All the statistical tests were performed in GraphPad Prism (version 9.3.1). Chi-square test was used for the comparison of the number of V2a neurons between different cervical regions. This was then followed up by Fisher’s exact test to compare V2a neurons between two different cervical regions at a time. A z-test of proportions was done to compare the proportion of neurons that are changed before and after injury. All data is represented as mean ± SD. The p value was considered significant if < 0.05. Figure 1 was made using BioRender.com.
Single-nucleus RNA sequencing of nuclei enriched for V2a neurons. A, Methodology for enrichment of V2a nuclei for RNA sequencing. Cervical spinal cord from V2a-Sun1GFP mice were harvested to enrich for V2a neurons using FANS. The sorted nuclei were then subjected to single-nuclei RNA sequencing and analyzed using computational tools. B, Confocal image of a representative cervical spinal cord section from a V2a-Sun1GFP mouse following immunohistochemistry with antibodies to Chx10 (red) and GFP (green). Scale bar, 100 µm. C, D, Boxed area in B showing coexpression of GFP and Chx10 in most cells. Yellow arrows denote GFP+/Chx10+ nuclei and white arrows denote GFP+/Chx10− nuclei. Scale bars, 50 µm. E, UMAP from Seurat showing clusters in uninjured dataset (both cohorts) identified cell type based on expression of genes in B prior to removal of non-V2a neuron nuclei. F, Dot plot showing the expression of pan neuronal, excitatory, inhibitory, V2a, non-V2a, astrocyte, microglia, and oligodendrocyte markers in the uninjured dataset. Blue represents high expression and gray represents low expression. The size of the dot represents the percentage of nuclei that express a particular marker in a cluster.
Results
Adult V2a neuron subpopulations are defined by unique molecular markers
To assess the heterogeneity of adult V2a murine neurons, we applied single-nucleus RNA sequencing (snRNA-seq) to genetically defined V2a spinal cord neurons (Wu et al., 2019). As V2a neurons represent a rare subset of cells in the spinal cord, we used the INTACT (Mo et al., 2015) protocol to enrich our nuclei samples for V2a. To facilitate enrichment and isolation of V2a nuclei, we generated V2a-Sun1GFP mice (Chx10Cre/+; ROSALSL-SUN1sfGFP/+) in which the nuclear membranes of cells that express the V2a marker gene Chx10 during development or adulthood are fluorescently tagged by a nuclear membrane bound green fluorescent protein (GFP). Nuclear preparations from cervical cords from these mice were sorted to enrich for GFP+ (V2a) nuclei using flow-activated nuclei sorting (FANS) and subsequent RNA sequencing (snRNA-seq) using the 10x Genomics Chromium platform (Fig. 1A). Only ∼200–300 GFP positive cells were identified per cord, necessitating isolation of 2–4 cords per cohort. To confirm GFP expression in V2a nuclei, we performed immunohistochemistry using antibodies to GFP and Chx10 on tissue sections from the cervical spinal cord (Fig. 1B–D). GFP+ nuclei are observed in the intermediate laminae of the spinal cord, consistent with the location of V2a neurons. As expected from previous studies (Zhong et al., 2010; Azim et al., 2014; Romer et al., 2017; Hayashi et al., 2018; Jensen et al., 2019), we observed some GFP+ nuclei that no longer express detectable Chx10 due to the downregulation of this protein postnatally in a subset of V2a neurons. Our results confirm the expected GFP expression pattern in V2a neurons.
We compared the effects of using more versus less stringent sorting parameters for enriching GFP+ nuclei in two separate cohorts of uninjured cords (see Materials and Methods for details). The less stringent approach was used to ensure that we did not lose specific subsets of V2a neurons due to sorting parameters. Each cohort was sequenced separately and integrated using the Seurat integration pipeline. This analysis identified 23 initial distinct nuclei populations using conservative clustering resolution options (Materials and Methods). As expected, we observed a small number of contaminant nuclei clusters, comprising nine glial cells clusters, one inhibitory neuron population, and three other neuronal cell populations that do not express well-defined V2a neuronal markers (Fig. 1E,F; Data S1). We noted that the majority of contaminant non-V2a nuclei had oligodendrocyte markers. After excluding contaminant populations, we identified 10 clusters of V2a neurons from a total of 1,251 nuclei in the cervical cord. V2a subpopulations were distinguished based on both unique gene expression profiles and clear separation by low-dimensional embedding of these profiles (Fig. 2A,B; Data S2). All of the clusters comprised at least 30 nuclei and each cluster contained nuclei from both cohorts, demonstrating reproducibility between samples (Fig. 2C). We identified 4,908 genes differentially expressed between V2a clusters and expressed in at least 25% of the cells per cluster (Wilcoxon rank sum test p < 0.05, FDR). The 10 V2a clusters were defined by distinct markers of axon guidance and migration (Sema5a, Doc1, Nav3, Atp2b4, Hmcn2, Erbb4), neuronal development (Tcf4, Rnf220, Pou6f2, Klhl1, Nfib), and embryonic neuronal patterning (Dach1, Dach2, Tshz2; Fig. 2A). Additionally, we found that genes of functional classes such as ion channels and transmembrane signaling proteins were differentially expressed between different clusters (Fig. 3). The heterogeneity of adult V2a neurons is further highlighted by the top gene enrichment study of each cluster showing that different pathways are enriched in each subset of V2a neurons (Fig. 3). These data will be a useful resource for further study of propriospinal neuron diversity and function.
Adult V2a neurons can be divided into 10 clusters based on differential gene expression. A, Dot plot showing top 5 differentially expressed genes for each cluster. Blue is high expression and gray is low expression. The size of the dot represents the percentage of nuclei expressing a particular marker in the cluster. B, UMAP of integrated (uninjured and injured) dataset showing diversity of V2a neurons (clusters 0–9) based on differences in gene expression. Each dot is one nucleus, and different colors represent different clusters. C, UMAP showing the 10 identified V2a clusters in each cohort from uninjured and injured cords.
V2a neuron subtypes show differential expression of functionally important genes. A–E, Heatmaps showing differential expression in V2a neurons of functional classes of genes involved in (A) axon guidance, (B) cell migration, (C) transmembrane receptor protein tyrosine, (D) chloride channel complex, and (E) potassium channel activity. Yellow is high expression and blue is low expression. F, Heatmap showing enrichment of gene ontology terminologies for each cluster. Red represents high z-score, and blue is low z-score.
N- and Z-group V2a neuron markers are maintained in adulthood
Previous studies have indicated that V2a neurons can be differentiated into two broad groups during development, the N- and Z-groups (Hayashi et al., 2018; Osseward et al., 2021; Sagner et al., 2021). These two groups of V2a neurons are generated at different times during development and distributed differently within the spinal cord. Z-group V2a neurons are early born and laterally enriched in the spinal cord whereas N-group V2a neurons are late-born and medially enriched in the spinal cord (Hayashi et al., 2018; Delile et al., 2019; Sagner et al., 2021). We assessed whether these N- and Z-group subdivisions of V2a neurons could be discerned in adult V2a neurons based on gene expression patterns of Z-group markers (Zfhx3, Zfhx4, Nefl, Shox2, Foxp2, and Otp) and N-group markers (Nfib, Tcf4, Prox1, Sp8, and Neurod2) identified in neonatal datasets (Hayashi et al., 2018; Osseward et al., 2021). We found that 7/10 V2a clusters expressed Z-group markers whereas 2/10 clusters expressed N-group markers (Fig. 4A). Cluster 4 was not categorized into either group because it did not express a definitive pattern of marker expression (i.e., only some nuclei in this cluster had low expression of the N-group marker Tcf4). Some neonatal N- and Z-group markers were poorly conserved in the adult dataset (Neurod2, Sp8, Prox1, FoxP2, Shox2) likely due to downregulation of these markers during postnatal development. In fact, Shox2, Sp8, and Foxp2 have been previously shown to be downregulated globally in the spinal cord (Morikawa et al., 2009; Dougherty et al., 2013). Based on the distance between clusters in the low-dimensional UMAP embeddings, our data shows that N-group clusters (2, 5) are molecularly distinct from Z-group clusters even into adulthood.
Limited conservation of neonatal and adult markers for V2a neurons. A, Dot plot showing average expression of markers in adult dataset previously used to distinguish N- and Z-group neurons in a neonatal dataset. Blue is high expression and gray is low expression. The size of the dot represents the percentage of nuclei expressing a particular marker in the cluster. B, Sankey plot showing correspondence of clusters identified in the uninjured adult V2a dataset to clusters identified from a published neonatal V2a dataset using the label transfer method. The thickness of flow lines indicates the number of nuclei/cells. C, Heatmap showing concordance matrix between clusters identified from adult dataset to clusters identified from neonatal dataset. Concordance ranges from 0 (black) to 1 (red). D, Sankey plot showing correspondence of clusters identified in the neonatal V2a dataset to clusters identified in the adult uninjured V2a dataset using the label transfer method. The thickness of flow lines indicates the number of nuclei/cells. E, Heatmap showing concordance matrix between clusters identified from the neonatal V2a dataset to clusters identified from the adult V2a dataset. Concordance ranges from 0 (black) to 1 (red).
We evaluated the concordance of V2a clusters between neonatal and adult time points using a single-cell label transfer approach in Seurat to define reciprocal populations from a published neonatal V2a dataset from cervical and lumbar cord (Hayashi et al., 2018) to our adult V2a dataset from cervical cord (Fig. 4B). This analysis found that 8/11 neonatal dataset clusters could be mapped onto our adult dataset with high confidence. Broadly, the N-group clusters from the neonatal dataset mapped to N-group markers in the adult dataset, whereas Z-group clusters in the neonatal mapped to Z-group markers in the adult. For example, the neonatal N-groups (N2 and N8) showed a high level of concordance with the adult N-group clusters (N2 and N5, respectively; Fig. 4B,C). However, individual neonatal Z-group clusters mapped onto multiple adult Z-group clusters and with a lower level of concordance for most clusters (Fig. 4C). Note that clusters Z0, Z1, and Z10 from the neonatal V2a dataset correspond to clusters that were previously identified as predominantly lumbar V2a neurons (clusters 1, 2, and 10 from Hayashi et al., 2018). We also performed label transfer using the adult markers onto the neonatal clusters and corroborated the finding that there is not a 1:1 correspondence of markers between neonatal and adult Z-group clusters (Fig. 4D,E; Data S3). These results are consistent with significant changes in gene expression occurring between neonatal and adult ages, including the loss of expression of many developmental markers.
We next validated the conservation of N- and Z-group marker expression in adult cervical spinal cord sections using in situ hybridization for Nfib and Zfhx4 mRNA (Fig. 5A–D). We found that 20 ± 4% and 72 ± 11% of the neurons expressed Nfib and Zfhx4 markers, respectively (n = 4 cords). In addition, we found that 9 ± 4% of V2a neurons expressed both Nfib and Zfhx4 mRNA (Fig. 5E). Limited coexpression of N- and Z-group markers was also found in neonatal spinal cord (Hayashi et al., 2018; Osseward et al., 2021). Next, we used immunohistochemistry to investigate protein expression and spatial location of these neurons in cervical spinal cord sections (Fig. 5F–O). We detected Nfib protein in 10 ± 2% of V2a (n = 7 cords) and Zfhx4 protein in 61 ± 5% of V2a (n = 7 cords). The relatively small discrepancies between the percentage of V2a expressing Nfib or Zfhx4 using snRNA-seq, RNAScope, and immunohistochemistry may be attributed to the differences in the sensitivities of each technique and/or differences in protein versus mRNA expression. To evaluate the distribution of Nfib+ and Zfhx4+ V2a, we mapped the location of these neurons in the cervical spinal cord and represented the data using contour plots (using the central canal as the origin). We found that Nfib-expressing V2a were distributed predominantly in the medial spinal cord and Zfhx4 expressing V2a were concentrated in the lateral cord, although with a small subset present in the medial cord (Fig. 5P). However, not all medial V2a expressed Nfib and some medially located V2a also expressed Zfhx4 (Fig. 5F–P). Thus, our results demonstrate that markers for N- and Z-group V2a are expressed in the adult spinal cord and follow a similar distribution pattern as in neonatal cords (Hayashi et al., 2018; Osseward et al., 2021).
N- and Z-group markers label medial and lateral V2a, respectively, in adult spinal cord. A–D, Confocal images of representative nuclei in cervical spinal cord tissue from V2a-Sun1GFP mice using in situ hybridization probes to detect RNA for Zfhx4 (white spots) and Nfib (red spots). Scale bars, 2 µm. E, Pie chart showing percentage of V2a nuclei expressing Zfhx4 and Nfib RNA in cervical spinal cord tissue assessed using in situ hybridization. F–O, Confocal images of cervical spinal cord sections from V2a-Sun1GFP mice immunostained with antibodies to Zfhx4 (red, F) or Nfib (red, K) and GFP (green, F–K) showing distribution of V2a in the medial (left) and lateral (right) spinal cord. G, H and L, M show individual neurons in the medial part of the spinal cord. I,J and N,O show individual neurons from the lateral part of the spinal cord. Yellow arrows denote V2a nuclei with Nfib or Zfhx4 expression white arrows denote V2a nuclei without marker expression. Scale bars: F and K = 100 µm, G–J and L–O = 50 µm. P, Contour plot (bottom) showing the density distribution of Nfib-expressing (red) and Zfhx4-expressing (blue) V2a nuclei in the cervical spinal cord based on immunohistochemistry. Contour indicates the density of neurons from 10th to 70th percentiles. Surface plot (top) represents the density in the mediolateral (M–L) plane. Q, Dot plot showing the clusters in which medial markers Nfib and Necab1 are expressed. Blue is high expression and gray is low expression. The size of the dot represents the percentage of nuclei expressing a particular marker in the cluster. R, Confocal image of cervical spinal cord section showing the expression of Necab1 (red) and Chx10 (white) detected by immunohistochemistry. Scale bar, 100 µm. Bottom left, Region of white dotted box showing Necab1+/Chx10+ (yellow arrows) and Necab1−/Chx10+ nuclei (white arrow). Scale bar, 5 µm. S, Contour plot (bottom) showing the density distribution of Nfib-expressing (red) and Necab1-expressing (black) Chx10+ neurons detected by immunohistochemistry. Contour indicates the density of neurons from 10th to 70th percentiles. Surface plot (top) represents the density in the mediolateral (M-L) plane. T, The percentage of Chx10+ neurons expressing Necab1 at C1-C2 (361 V2a neurons), C3-C5 (909 V2a neurons), and C6-C8 (580 V2a neurons) spinal segments. Fisher’s exact test: ***p = 0.003, ****p < 0.0001.
Our analysis reveals that the N-group can be further divided into 2 clusters (N2 and N5) distinguishable based on cluster-specific markers, such as the N5 marker Necab1 (Fig. 5Q). To validate that Necab1 is expressed in a subset of medial V2a neurons, we performed immunohistochemistry on cervical spinal cord tissue using antibodies against Necab1 and Chx10 (to label most V2a neurons; Fig. 5R). We found that 8 ± 2% (n = 4 cords) of Chx10+ neurons expressed Necab1, comparable with the 12% of neurons predicted by our RNA sequencing data. Necab1+/Chx10+ neurons were concentrated in the medial part of the cord, as expected (Fig. 5S). Interestingly, a higher fraction of the Chx10+ cells in the lower cervical (C6-C8) cord were Necab1+ compared with the upper cervical (C1-C2 and C3-C5; Fig. 5T), indicating that the distribution of this subset of V2a neurons varies by segmental level. This segmental difference was not observed in Zfhx4+ V2a neurons (p = 0.84, C1-C2: 401/653, C3-C5: 1229/2035, and C6-C8: 1051/1718 Zfhx4+ V2a/total V2a). Thus, our results show that V2a subtypes may be enriched in specific regions along both the medial-lateral and rostrocaudal axes of the spinal cord.
Specific V2a subsets are lost after injury
We examined gene expression changes in V2a neurons in the subacute phase (5 d) after a C2 hemisection injury. Similar to our analysis of nuclei from uninjured cords, spinal cords from lesioned mice were collected below the injury site, nuclei isolated and sorted by FANS, snRNA-seq performed using the 10x Genomics chromium platform (using 2 separate cohorts of injured cords), and quality control steps performed to eliminate low-quality neurons and non-V2a neurons (see Materials and Methods for details). As new neuronal subsets following spinal cord injury were not identified in prior studies (Russ et al., 2021; Kathe et al., 2022; Matson et al., 2022; Yadav et al., 2023), we projected snRNA-Seq profiles of cells from uninjured cords into the integrated (uninjured + injured) dataset using a label transfer approach (Fig. 6A, Data S4). While the majority of uninjured V2a specific clusters were retained in injury (at least 30% similarity in their proportion), we surprisingly find that three of our clusters were entirely lost (Tshz2+/Sema5a+ cluster Z7) or almost entirely lost (Gpc5+/Tshz2+ cluster Z8 and Atp2b4+/Kcnmb4+ cluster Z9; Fig. 6A,B). All three lost clusters were found to highly express the multifunction RNA-binding protein and translational regulator Rbms3 and the axon growth regulator Nav3 (Fig. 2A). Rbms3 is a notable regulator of cell cycle and apoptosis (Górnicki et al., 2022). Clusters Z7 and Z8 (but not Z9) express the axon guidance molecule Sema5a+. Thus, specific subsets of V2a neurons below the site of injury are either lost or change their cluster identity within 5 d of injury.
Specific subsets of Z-group V2a neurons are lost both ipsilateral and contralateral to injury. A, UMAP plots showing V2a neurons from uninjured (left) and injured (right) cords integrated into the same dataset. Each dot is one nucleus, and each color represents different clusters. Red boxes highlight clusters whose neurons are predominantly found in uninjured but not injured cords. B, Bar graph showing the percentage of V2a neurons from uninjured (blue) or injured (red) cords found in each cluster. C–H, Representative confocal images of GFP+ nuclei (green) in cervical spinal cord tissue from V2a-Sun1GFP mice using in situ hybridization probes to detect RNA for Pde3a (red spots) and Sema5a (white spots) (C–F) or Syt10 (red spots; G,H). Scale bars, 2 µm. I, Bar graph showing the percentage of Pde3a−/Sema5a+, Pde3a+/Sema5a−, Pde3a+/ Sema5a+, and Syt10+ V2a nuclei in the cervical spinal cord based on in situ hybridization. Data is represented as mean ± standard deviation (n = 4 animals, 538 V2a for Pde3a/Sema5a injured, n = 3 animals, 698 V2a for Pde3a/Sema5a uninjured, n = 3 animals, 515 V2a for Syt10 uninjured, and n = 3 animals, 555 V2a for Syt10 injured). J, Fraction of Pde3a-/Sema5a+ (Z7 and Z8 markers) V2a neurons found on the ipsilateral and contralateral cord in uninjured (blue) and C2Hx injured (red) mice. K, Top 10 gene enrichment terms (Gene Ontology and Pathway) for marker genes of clusters Z7, 8, and Z9. L, Violin plot showing distribution of regeneration module scores for V2a in each cluster before and after injury. Points represent model-based mean scores and error bars show 95% confidence intervals. Z7-9 are not included due to the low number of nuclei after injury.
We tested whether the lost V2a subsets were indeed missing in vivo using in situ hybridization to detect Sema5a+ V2a (Fig. 6C,D). As Sema5a is not only expressed in the injury ablated populations (clusters Z7 and Z8) but also in cluster Z6 which is retained post injury, we simultaneously monitored a selective cluster Z6 marker, Pde3a (Fig. 6E,F). Further, we examined Syt10+ V2a (cluster U4) that are also retained after injury (Fig. 6G,H). Prior to injury, Sema5a (11 ± 1% of V2a), Pde3a (5 ± 1%), and Syt10 (3 ± 1%) mark subsets of V2a as predicted from our uninjured adult V2a snRNA-Seq analyses. As expected for Z-group neurons, Pde3a+ and Sema5a+ V2a neurons were observed in both lateral and medial spinal cord. However, upon injury the number of Pde3a-/Sema5a+ V2a (clusters Z7 and Z8) is significantly decreased relative to the uninjured spinal cord (p < 0.001), with no significant change in number of cluster U4 (Syt10+) or Z6 (Pde3a+/Sema5a+) V2a (Fig. 6I). We tested the hypothesis that the Pde3a−/Sema5a+ V2a population might be selectively reduced on the side ipsilateral to injury but instead found that there is a reduction of Pde3a−/Sema5a+ V2a both ipsilateral and contralateral to injury (Fig. 6J). A two-way ANOVA showed that there were main effects of side (ipsilateral vs contralateral to injury, F(1,10) = 10.77, p = 0.0083) and group (uninjured vs injured, F(1,10) = 39.93, p = 8.7 × 10−5). There was no interaction between groups and side (F(1,10) = 1.07, p = 0.32). Holm–Sidak post hoc analysis revealed significant differences between injured and uninjured on the left (ipsi) side (p = 0.015) and right (contra) side (p = 0.0019). There was no significant difference between left and right sides within uninjured (blue; p = 0.0502) or injured (red; p = 0.12) groups. These results indicate that the reduction in Pde3a−/Sema5a+ V2a occurs both ipsilateral and contralateral to injury and is consistent with our snRNA-Seq data showing a near complete loss of cluster Z7 and Z8 V2a rather than the ∼50% loss expected if only ipsilateral V2a were lost.
To distinguish between loss of V2a neurons via cell death versus a change in cluster identity, we performed antibody staining for cleaved caspase-3 (a marker of apoptotic cells) 2 d following a C2Hx in V2a-Sun1GFP mice. Out of 3,097 V2a neurons counted below the site of injury (C3-C8, ipsilateral and contralateral to injury), none were labeled with cleaved caspase-3 (n = 5 mice). However, we did observe cleaved caspase-3 staining in sections at/near the site of injury at this time point, including 4 V2a out of 431 counted ipsilateral to injury (n = 5 mice), demonstrating that our antibody was effective in detecting dying neurons at this time point. These data, together with the observation that Z1, Z3, and U4 increase in number after injury (Fig. 6B), suggest that V2a neurons in clusters Z7-9 alter their identity rather than undergo programmed cell death.
We next examined the genes/pathways enriched in the cells most vulnerable to changes after injury (lost clusters Z7, Z8, and Z9). We identified a total of 152 genes enriched in clusters Z7, Z8, and Z9 compared with the other clusters with p adjusted value <0.05, where the top 5 enriched genes were Grm8, Galntl6, Plxdc2, Nkain2, and Samd5 (ranked by log2 fold change). Performing gene enrichment analysis on the identified 152 genes, we found neuron projection, synaptic transmission, ionotropic glutamate receptor complex, glutamatergic synapse, and cell–cell adhesion terminologies to be associated with these clusters (Fig. 6K). This analysis suggests that differences in expression of synaptic transmission related genes in these clusters could potentially make them more vulnerable to gene expression changes (e.g., shift in cluster identity) or cell loss after injury.
We hypothesized that damage caused by injury to V2a neurons with ascending projections to the brainstem could alter expression of RAGs. Using a list of 525 RAGs identified in a previous study (Saraswathy et al., 2024), we identified 373 RAGs expressed in V2a neurons. We calculated a regeneration module score for V2a neurons in each cluster in both injured and uninjured datasets (Fig. 6L, Data S5). We used a linear mixed model that accounts for variation between replicates, to test the effects of injury, cluster identity, and their interaction. The analysis showed that regeneration scores differ significantly between clusters (p = 2.2 × 10−16), indicating that some clusters naturally have higher or lower baseline expression of RAGs. Clusters Z7, Z8, and Z9, however, did not differ significantly from baseline (cluster Z0, p = 0.13, 0.96, 0.82, respectively). Further, there was no significant overall effect of injury (p = 0.94), and post hoc pairwise comparisons within each cluster showed that injury did not significantly change the regeneration score in any cluster. This suggests that baseline differences between clusters are more pronounced than any changes induced by injury. Thus, the C2Hx injury did not strongly alter the expression of RAGs in V2a at 5 d postinjury.
Different subsets of V2a neurons show distinct gene expression changes following spinal cord injury
While three of the snRNA-Seq defined V2a subsets were effectively lost after injury, in the remaining V2a clusters we aimed to quantify transcriptional perturbations in order to assess which V2a subsets were most impacted by C2 hemisection. We applied bioinformatics approaches to assess both global (Augur; Skinnider et al., 2021) and cell type-specific (cellHarmony; DePasquale et al., 2019) perturbations. Augur applies machine learning to quantify the differences between unperturbed and perturbed cells within a high dimensional space and represents the degree of perturbation with a value denoted as the AUC. Using this method, we found that some V2a clusters were more perturbed after spinal cord injury than others (Fig. 7A). For example, cluster N5 (Necab1+) showed the highest perturbation with AUC value of 0.802 whereas cluster N2 was least perturbed with value of 0.622. We found high and low AUC clusters within both N- and Z-groups. Thus, even the V2a clusters that are not reduced in cell number after injury show differences in the degree of change that they experience after injury.
Diverse transcriptional changes in V2a subtypes following spinal cord injury. A, UMAP showing cell type prioritization by Augur where red is highest AUC score, and blue is lowest AUC score. B, Bar plot showing the number of genes altered following injury in different V2a clusters from cellHarmony. Fold change cutoff of 1.2 and adjusted p value of 0.05 was used. Blue is downregulated and red is upregulated. C, Heatmap showing changes in the gene expression after injury for each V2a cluster. Each line is a gene whose expression is altered after injury. Blue is downregulated after injury and yellow is upregulated after injury. Gene Ontology enrichment per cluster is represented on the left side. D, WikiPathways (WP2267) visualized in AltAnalyze for differentially expressed genes in Cluster Z0 (injury vs uninjured). Blue is downregulated and red is upregulated. E, Network of differentially expressed genes in cluster Z0 (injury vs uninjured). Red dots represent upregulation, blue dots represent downregulation, blue arrow is inhibitory regulation, and red arrow is transcriptional regulation. F, Dot plot showing the downregulation of genes with largest fold changes in V2a neurons following injury. The size of the dot represents the percentage of nuclei expressing a particular marker in the cluster. G, Dot plot showing upregulation of splicing factors across all conditions. The size of the dot represents the percentage of nuclei expressing a particular marker in the cluster.
We next used cellHarmony to assess cell type-specific gene expression perturbations and impacted regulatory networks in different V2a clusters. For this analysis we also included cluster Z8, despite having only seven cells, in order to identify perturbations in these cells. We found that in all clusters, more genes are downregulated than upregulated after injury (Fig. 7B), consistent with results using different spinal cord injury models (Siebert et al., 2010; Garcia et al., 2018; Matson et al., 2022). Cluster Z0 (March1+/Col25a1−) was found to have the greatest number of genes altered following injury whereas cluster U4 (Col25a1+/Syt10+) to have the least (Fig. 7B). Altered synaptic vesicle signaling dynamics is particularly evident within cluster Z0, where many synaptic vesicle genes are downregulated or, in the case of Dnm2 (involved in endocytosis), upregulated (Fig. 7C,D). Furthermore, additional gene alterations following injury within this cluster are intricately associated with ion channel activity and cytoskeletal organization (Fig. 7C,E; Data S6). Notably, cluster Z8 (Gpc5+/Tshz2+) was found to have mostly upregulation of genes post injury.
We performed a comprehensive analysis of genetic perturbations following injury ordered by p value rankings and arranged by individual V2a clusters to reveal differential gene alterations between V2a clusters after injury (Fig. 7C). Only a minority of genes are up- or downregulated in all V2a clusters after injury. For example, we found an upregulation of splicing regulators in most V2a subsets, indicating alterations in mRNA diversity following injury (Fig. 7F). Likewise, we found that the most highly downregulated genes in V2a neurons were neuroprotective genes Lrp1b (Spuch et al., 2012; Benoit et al., 2013), Fgf1 (Tsai et al., 2015; Li et al., 2018), and Ube3a (Olabarria et al., 2019; Tonazzini et al., 2019; Fig. 7G). However, the majority of genes are altered only in specific clusters of V2a neurons, with some genes/pathways altered in opposite directions in different clusters. For example, cluster Z0 has a decrease in expression of genes involved in synapses and cell–cell junctions, whereas Z1 has an increase in synapse assembly and cell–cell junction assembly genes. Consistent with our analyses using Augur, these results show that different subsets of V2a neurons show different changes in gene expression patterns following spinal cord injury.
Discussion
The changes that occur in neurons below a spinal cord injury due the loss of descending and ascending connections between the spinal cord and brain are complex and not well understood. We focused on one class of propriospinal neuron (the V2a class) which prior studies have shown are important for recovery of function after spinal cord injury (Garcia-Ramirez et al., 2021; Skinnider et al., 2021; Huang et al., 2022; Squair et al., 2023; Van Steenbergen et al., 2023) to assess the gene expression changes at the single cell level 5 d after a C2 hemisection injury. We used single-nucleus sequencing, bioinformatics analysis techniques, and in situ hybridization to demonstrate that specific subsets of V2a are lost below the site of injury. We also show that the remaining clusters of V2a neurons show distinct changes in gene expression after injury, including upregulation of RNA splicing factors. These findings establish the importance of examining changes in gene expression at the single neuron level and provide a resource for assessing which neurons and pathways are most critical for recovery.
Our analyses of V2a nuclei from uninjured adult cervical spinal cord builds on work investigating gene expression in developing V2a neurons (Hayashi et al., 2018; Osseward et al., 2021; Sagner et al., 2021) and provides markers to facilitate characterizing the electrophysiological properties, location, and connectivity of adult V2a neuron subtypes. Our findings revealed a conservation of N- and Z-groups in adult spinal cord and identified specific markers for N-group (e.g., Nfib, Necab1) and Z-group (e.g., Zfhx4, Sema5a) clusters that are expressed in adulthood. Consistent with neonatal studies, we found an enrichment of N-group V2a in the medial spinal cord and Z-group V2a in the lateral cord. Intriguingly, the concordance between individual adult and neonatal clusters was not as strong as between N- versus Z- groups as a whole, since most adult Z-group clusters contained markers found in multiple neonatal clusters. This result could reflect the greater diversity found within the Z-group neurons and/or Z-group V2a may undergo more diverse changes in gene expression during postnatal development than N-group V2a. It is also important to note that our analysis compared adult V2a from cervical cord to neonatal V2a from both cervical and lumbar cord. Not surprisingly, one predominantly lumbar V2a cluster from the neonatal dataset (Z10) was not assigned a corresponding cluster in the adult cervical cord. However, the other two predominantly lumbar V2a clusters (Z0, Z1) were assigned similarity to adult cervical clusters, and two clusters containing cervical neonatal V2a (Z5, Z6) were not assigned corresponding adult V2a clusters. We hypothesize this result is due to the many changes in gene expression that occur between neonatal and adult stages and perhaps less reliance on positional cues in the adult. Our analyses provide a framework for identifying adult markers for specific developmental neuron subsets (and vice versa), which will be useful in future studies examining their specific roles in motor behaviors and recovery from injury.
Surprisingly, we observed differences in the distribution of Necab1+ (cluster N5) V2a along the rostrocaudal axis of the cord. Enrichment of N5 neurons in upper cervical segments could provide clues to their function. For example, upper cervical neurons have been described that project to both phrenic and intercostal muscles to control breathing (Lipski and Duffin, 1986; Lu et al., 2004). Conversely, long descending propriospinal neurons that project to lumbar cord (and include V2a neurons) are more prevalent in lower cervical and thoracic segments (Ni et al., 2014; Ronzano et al., 2021). Future experiments may test the hypothesis that N5 V2a neurons are more likely to contribute to the former rather than the later population of propriospinal neuron.
Following spinal cord injury, we observed downregulation of many genes in V2a neurons, more so than upregulation of genes. Prevalent downregulation of genes has also been noted in previous studies using different spinal cord injury models (Siebert et al., 2010; Garcia et al., 2018; Matson et al., 2022). Notably, we found an upregulation of splicing regulators 5 d after injury in V2a neurons. Splicing of mRNA is involved in generating transcriptomic diversity and plays important roles in synaptic plasticity, axon guidance, neuronal differentiation, maintenance of synaptic properties, apoptosis, and many disease pathologies (Raj and Blencowe, 2015; Iijima and Yoshimura, 2019; Traunmüller et al., 2023). A previous study of rat sacral cord tissue showed increased retention of introns after injury, consistent with changes in splicing factor expression (Hart et al., 2022). Specific changes in splicing have also been reported in rodent SCI models, including altered splicing of the serotonin 2c receptor and the NR1 subunit of the NMDA receptor (Prybylowski et al., 2001; Nakae et al., 2013; Stamm et al., 2017). The identification of specific splicing factors upregulated in V2a neurons after injury will provide a basis for testing which RNAs are altered by each splicing factor and whether those changes are protective or maladaptive.
Consistent with our hypothesis, we found that different subsets of V2a neurons showed disparate changes in response to injury. Remarkably, three clusters of V2a neurons (clusters Z7-9) are largely absent in the cervical cord following injury. Our results showing no detectable cleaved caspase-3 staining in V2a neurons 2 d after injury indicates that the V2a clusters are not lost due to programmed cell death. However, we cannot rule out that V2a might die by other mechanisms, such as necroptosis or ferroptosis. Alternatively, the loss of clusters Z7-9 may reflect changes in gene expression that result in altered cluster identity. This would be consistent with our previous study showing no significant loss of V2a below a C2Hx 1 d or 2 weeks after injury (Jensen et al., 2024) as well as a study showing no significant loss of glutamatergic neurons below a C2Hx 8 weeks after injury (Brezinski et al., 2025). The increase in the proportion of nuclei in clusters Z1, Z3, and U4 supports the hypothesis that neurons from the lost clusters are assigned to new clusters after injury.
Our findings show that changes in V2a neuron gene expression are not restricted to the side to injury. For example, rather than a 50% loss of V2a in clusters Z7-9, we see very few cells mapping to this cluster 5 d after injury. Further, we validated that Pde3a−/Sema5a+ V2a (Z7-8) are reduced on both the ipsilateral and contralateral side of the cord. Our results show that injury to one side of the cord affects neurons below the injury on both sides of the cord and thus it should not be assumed that circuits contralateral to injury are “normal.” Further, this result argues against axon damage to ascending V2a being a main cause of changes in gene expression, since this would affect only V2a ipsilateral to injury. We postulated that axon damage could activate RAGs in V2a after injury. However, we did not see a difference in RAGs in V2a 5 d after injury, consistent with a previous study showing that V2a are not among the neuron types that express a RAG signature (albeit a different RAG gene set) 1–3 weeks after a thoracic contusion injury (Matson et al., 2022). Intriguingly, we found differences in RAG expression score between clusters prior to injury, suggesting that some V2a subtypes may have inherently more capacity for regeneration than others, but this hypothesis will need to be tested in future studies.
Gene enrichment analysis of marker genes for the lost clusters identified terminologies related to glutamatergic synapses as well as neuron projection and cell–cell adhesion. Glutamatergic signaling is particularly interesting in this context as it has been implicated in promoting secondary damage following SCI (Rowland et al., 2008; Alizadeh et al., 2019). Downregulation of genes involved in glutamatergic signaling after injury could potentially serve as a protective mechanism against excitotoxicity that results in Z7-9 V2a converting to a different V2a subtype. However, future experiments are warranted to assess how this change in identity and/or altered glutamatergic signaling might affect the function of V2a within spinal circuits.
We used two different methods (Augur and cellHarmony) to assess the degree and types of gene expression changes that occur following injury in the V2a subsets that are preserved in both injured and uninjured tissue. Although the exact rankings of which clusters show the most perturbation (Augur) or number of genes altered (cellHarmony) are slightly different due to the different methodologies, they both found that cluster Z1 is one of the most perturbed and N2 and U4 are least perturbed by injury. Gene enrichment analyses show that the processes that are enriched in each cluster prior to injury differ, as do the processes that are up- or downregulated in each cluster after injury. Thus, our data shows that there are clear differences between V2a clusters on the impact of injury on transcriptional alterations.
Our results demonstrate the utility of using an in-depth approach to investigate the diversity of gene expression changes in neurons following injury and provide a foundation for investigating genes and pathways that underlie the differential vulnerability to injury, regenerative capacity, and/or plasticity of propriospinal neurons.
Data Availability
Single-nuclei RNA-seq data have been deposited at GEO (GSE269023; reviewer token: kpkxquqofvorxun) and are publicly available as of the date of publication. Any additional information required to reanalyze the data reported in this paper is available from synapse.org (syn53216432).
Footnotes
This research was made possible through the support of Craig H. Neilsen Foundation grants #598928 (S.A.C.), the L.B. Research and Education Foundation (I.W.), and the National Institutes of Health R01NS112255 (S.A.C.). We thank Dr. Mark Baccei and Elizabeth Serafin for generously sharing the nuclei isolation protocol and providing helpful advice; Kassidy Grover and Aditi Tarkar for their invaluable assistance in antibody testing; the dedicated team at the Research Flow Cytometry Facility (RRID:SCR_022635), Division of Rheumatology; Dr. Matthew Kofron and Sarah McLeod at Bio-Imaging and Analysis Facility (RRID:SCR_022628); Kelly Rangel and Shawn Smith at Single Cell Genomics Facility (RRID:SCR_022653); and David Fletcher at Genomics Sequencing Facility (RRID:SCR_022630) at Cincinnati Children’s Hospital Medical Center. The authors acknowledge that the content is solely their responsibility and does not necessarily reflect the official views of the NIH or other funding agencies.
The authors declare no competing financial interests.
This paper contains supplemental material available at: https://doi.org/10.1523/JNEUROSCI.1163-24.2025
- Correspondence should be addressed to Steven A. Crone at Steven.Crone{at}cchmc.org.
