Methods
Use of Xenopus tropicalis was carried out in accordance with ethical regulations and under the approval and oversight of the IACUC committee and Office of Animal Welfare at UW, an AALAC-accredited institution. Generation of tadpoles was carried out according to published methods 33 . In this study, we used both wild-type frogs and frogs from the triple transgenic line Xtr.Tg( pax6 :GFP; cryga :RFP; actc1 :RFP) (RRID:NXR_1.0021 37 ) reared and purchased from the National Xenopus Resource ( https://www.mbl.edu/xenopus/ ; RRID:SCR_013731). pax6 transgenic matings were performed by crossing a heterozygous transgenic female frog to a wild-type male frog. These matings yielded clutches with 50/50 wt/ pax6 :GFP+ populations. Since Xenopus tropicalis are not identifiable by sex either morphologically or genotypically at this stage, tadpoles were not sorted by sex, and thus sex data was not collected. For amputation assay, NF stage 41 tadpoles were anesthetized with 0.05% MS-222 in 1/9× MR and tested for response to touch prior to amputation surgery. Once fully anesthetized, a sterilized scalpel was used to amputate the posterior third of the tail. Amputated tadpoles were removed from anesthetic media within 10 min of amputation into new 1/9× MR. Tadpoles were kept at a density of no more than 1 tadpoles per 2.5 mL. Tadpoles were fed daily sera Micron 109 at ¼ mL per 2.5 L, starting at 3 dpa until final timepoint collection.
To measure regenerated length of tail, samples were mounted and imaged with a Leica M205FA fluorescence Stereomicroscope. Using the Leica imaging software (LasX), measurements were taken from the vent to the most posterior point of the tail. Regenerating measurements were taken based off of morphological identification of the amputation plane.
Tadpoles were anesthetized with 0.05% MS-222 in 1/9× MR and tested for response to touch prior to amputation surgery. Once fully anesthetized, a sterilized scalpel was used to amputate and collect intended tissue (Fig. 1A ). Tail tissues were collected within 30 min of amputation and added to a vial with 72 µl trypLE (ThermoFisher Scientific:12604013) and 3 µl 100 mg/mL collagenase P (ROCHE:11213857001) for every 100 tail tissues. Samples were moved to a 28 °C water bath for dissociation and triturated with a p200 every 5 min until dissociated fully and no chunks was visible (30–60 min). To quench the reaction, 150 µl of ice cold 10% FBS was added for every 100 tail tissues. Samples were then spun down for 5 min at 1000 × g . Supernatant was disposed and cells were resuspended in 300 µl 10% FBS for FACS. For each sample, 22–294 tadpoles were used (tail tip: 294 tadpoles; un1 dpa: 222 tadpoles; un3 dpa: 133 tadpoles; un7 dpa: 22 tadpoles; mid-trunk: 294 tadpoles; 1 dpa: 208 tadpoles; 3 dpa: 144 tadpoles; 7 dpa: 25 tadpoles). Tail tip and mid-trunk tissues were collected from the same specimens.
Cell sorting and collecting were performed with an Aria III Cell Sorter (BD Biosciences) and a 70 µm nozzle. Live singlets were sorted for based size and granularity with FSC and SSC (Fig. S 1C, D ). For GFP gate setting, initial setting was done with positive control transgenic Xtr.Et(eef1a1:GFP) 16FMead tadpoles to test for GFP signal and autofluorescence. Selection areas were drawn to exclude most autofluorescent GFP- cells. Cells were collected into 500 µl solutions of 1× PBS with 1 mg/ml BSA (NEB:B9200S) and counted with a Countess 3 FL Automated Cell Counter (ThermoFisher Scientific). Single-cell mRNA libraries were prepared using the GEM-X single cell 3’ v4 kit (10× Genomics), with a target capture of 10,000 cells. Quality control and quantification assays were performed using a Qubit fluorometer (Thermofisher) and a Screentape Assay (Agilent). Libraries were sequenced on an Illumina NextSeq2000 using a 100-cycle, P4 XLEAP-SBS™ kit. Each sample was sequenced to an average read depth of 226 million total reads, which resulted in an average read depth of ~28,000 reads/cell after normalization. The raw FASTQ files generated during the current study are available under GSA #CRA033265 and on SRA under #PRJNA1448102. The processed Seurat datasets are available on GEO under # GSE327206 . Source Data of non scRNA-seq figure data are provided with this paper.
RNA sequencing reads were processed using CellRanger v9.01 Count and Aggregate by 10x Genomics. To generate the reference genome, the Xenopus tropicalis v10.0 reference genome and gene models were downloaded from Xenbase 110 . The output from CellRanger was processed using DoubletFinder (v2.0.6) 111 to label doublets using code adapted from https://github.com/ankita16lawarde/scRNA-Seq_endometriosis . The data was then run through the Seurat v3.1 pipeline 38 in R Studio v4.4.3. Cells expressing <200 genes were removed from downstream analysis, with a mitochondrial RNA read cutoff of 25%. The data was normalized, then scaled by Sample and nFeature_RNA before clustering. Data was integrated by Sample using RPCA Integration. Three poor-quality clusters indicated by a high percentage of cells with abnormally low UMI count (nFeature <500) were manually removed. Cell cycle was predicted following the Seurat Vignette “Cell-Cycle Scoring and Regression”. Cell types were annotated using common cell markers as well as via analysis of differentially expressed cluster markers (Fig. S 4 ). Post analysis, the final dataset contained 35,852 cells (Fig. S 1D ). ggplot2 (v4.0.1) was used to generate plots 112 .
The neural dataset was subset based on the expression of sox2 , pax6 , and tubb2b in neural progenitors and/or neurons. The neural subset was re-normalized, scaled by Sample and nFeature_RNA, then integrated by Sample with RPCA. Five clusters were manually removed based on their presence in only one sample (1 cluster), their expression of sclerotome and/or muscle markers (3 clusters), or abnormally low UMI count (nFeature <500) (1 cluster). Post analysis, the neural dataset contained 7315 cells. Go.db (v3.21) was used to evaluate Gene Ontology 54 , 55 terms during regeneration, and circlize (v0.6.16) 113 was used to generate chord graphs of gene-term connectivity. Regeneration-specific subclusters were either identified based on a cluster population of 100% regenerating cells ( prdm14 − MNs) or manually based on cluster sub-groups of 100% regenerating cells ( prdm14 + MNs, V0 INs, V2a/b INs, RBNs). Groups of less than 10 cells were not included. CellChat (v1.4.0) 114 was used to evaluate atf3 + regenerating neuron cell-cell interactions. Only Xenopus genes with 1-to-1 ortholog mapping, as defined on Xenbase 110 , were mapped to mouse orthologs.
To analyze the zebrafish iNeuron dataset from ref. 32 , we followed the same Seurat pipeline as described above, with the removal of cells expressing <200 genes. The iNeuron cluster was identified by expression of atf3 . Top differentially expressed markers (Fig. S 9A–C ) were calculated using FindAllMarkers from Seurat (v5.4.0). All relevant scRNAseq code will be deposited on GitHub upon publication under https://github.com/angellswearera/scRNAseq_xenopus_regen_wills .
For HCR, tadpoles were fixed overnight in 1× MEM with 3.7% formaldehyde at 4 °C. HCR was done in whole mount from established protocols 115 , at a probe concentration of 16 nM. Post HCR protocol, samples were washed with 3 × 15 min 1× PBS, 1 × 15 min 15% sucrose in 1× PBS, then overnight at 4 °C. To identify regenerating tissue post-sectioning, fins were trimmed dorsal and ventrally only posterior to the amputation plane. Samples were then transferred into OCT (Thermofisher 23-750-571) and kept at −80 °C until sectioning with a Leica CM3050 S Cryostat. Post sectioning, slides were baked at 37 °C for 20 min before imaging using a Leica SP8 Confocal with a 63× objective. They were processed using FIJI image analysis software. To measure Dorsal/Ventral location, the center point of cells was measured from the bottom of the spinal cord floor plate, then normalized to spinal cord height, as measured from the bottom of the floor plate to the top of the roof plate. ggplot2 (v4.0.1) was used to generate hex plots 112 . HCR probe sequences can be located in the Source Data file.
In order to assess the level of proliferative neurogenesis during regeneration, we used the DNA-synthesis marker BrdU. Tadpoles at the indicated stage were injected with 16 nL of BrdU (10 mM, Fisher, B23151 ) in the gills, following the protocol established Xenopus for tracking neurogenesis during development 57 . After injection, tadpoles were allowed to regenerate for 2 days to allow for detectable differentiation into neurons via HuC/D staining. Staining proceeded as previously described in ref. 26 , (1° Abs: 1:200 Rb anti-BrdU polyclonal, Fisher, PA5-32256; 1:200 Ms HuC/D monoclonal, Invitrogen, 16A11) with the addition of a 10 min 2 N HCl acid wash following the initial permeabilization in PBS-Triton. Because the BrdU 10 min 2 N HCl wash destroys HCR signal, we iteratively stained for HCR and BrdU. Samples were injected with BrdU as described, then run through the HCR protocol, sectioned, and imaged as described above (“Imaging 1”, Fig. 5I, J ). Slides then underwent staining for BrdU and re-imaging using the described protocol, with solution volumes adjusted for slide immunohistochemistry (“Imaging 2”, Fig. 5I, J ).
Tadpoles were stained and imaged as previously described in ref. 26 , (1° Abs: 1:200 Ms Tyrosine Hydroxylase (TH) monoclonal, Immunostar, 22941 or 1:200 Ms HuC/HuD monoclonal, Invitrogen, 16A11). Co-staining with BrdU was conducted as described above.
Tadpoles were treated by immersion in 5 µM Dclk1-in1 (Tocris cat #7285; stock solution was 50 mM in DMSO, diluent for immersion was 1/9MR) either immediately following amputation or at the corresponding uninjured timepoint. 0.01% DMSO was used as a vehicle control. Inhibitor solution was refreshed on days 1, 3, and 5. For inhibition of proliferation by a combined 150 µM Aphidicolin (Sigma A0781; stock solution was 150 mM in DMSO, diluent for immersion was 1/9MR) and 20 mM HUA (Sigma H8627; stock solution was 10 M in water, diluent for immersion was 1/9MR). 0.1% DMSO was used as a vehicle control.
In order to track migration and differentiate between newly generated neuronal cells, we used a commercial pMT-HuC: Dendra2 tol2 plasmid system to label post-mitotic neurons 63 . The tadpoles were injected at the one-cell and two-cell embryo stages with a total of 2 nanoliters of 14.87 ng/µl Dendra2 plasmid ((Addgene Plasmid #80904) and 5.19 ng/µl tol2 RNA (pCS-TP (Tol2)), with final amounts of 29.74 pg Dendra2 and 10.38 pg tol2 RNA per animal. After injection, Xenopus tropicalis tadpoles were grown to stage 41. Photoconversions were done using a Leica DM5500 microscope by directing violet (405) light on a posterior section of the tail, including ~100 µm posterior to the amputation plane, for 3–5 min. Pictures were taken before and after photoconversion. Then the distal 1/3 of the tail was amputated using a scalpel. Animals were then kept in the dark. At 3 dpa, the tails were imaged in the 488 and 594 nm channels.
All statistics were performed in using R Stats package (v3.6.2).
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Results
Our principal goal was to catalog changes in neural cell type composition in the spinal cord during tail regeneration and compare these with changes over a corresponding developmental timeline. To this end, we performed single-cell RNA sequencing on injured and stage-matched uninjured tail tissue from 0 to 7 days post amputation (dpa), beginning at stage 41 (3 days post fertilization) (Fig. 1A ) 36 . By 7 dpa, tail length in injured tadpoles matched uninjured peers, and we considered regeneration complete (Fig. S 1A, B ). To enrich for neural cells, we used the transgenic line Xtr.Tg( pax6 :GFP;cryga:RFP; actc1 :RFP) 37 , which broadly marks NPCs with GFP in X. tropicalis 33 . At stage 41, we amputated tadpoles and collected the posterior third of the tail as our “tail tip” sample and ~100 microns of tissue proximal to the amputation plane as our “mid-trunk” sample to control for anterior/posterior differences in cell composition (Fig. 1A ). We then split tadpoles into uninjured and regenerating groups and collected tail tissue from both groups at 1, 3, and 7 dpa to represent early, intermediate, and late spinal cord regeneration, respectively. Permissive FACS gating conditions were used to enrich for high-expressing pax6 :GFP+ cells (such as NPCs) as well as low-expressing pax6 :GFP+ (such as differentiated neurons) and some pax6 :GFP- cells, allowing us to populate non-neural and non- pax6 clusters (Fig. S 1C, D , “Methods”). After quality control filtering, a total of 20,450 uninjured and 15,402 regenerating cells were included for downstream analysis with Seurat (Figs. S 1E and S 2A–I ) 38 . After dataset integration and UMAP dimensional reduction, we identified clusters, to which we assigned unique cell type labels based on previous scRNA-seq analyses in Xenopus and other relevant developmental markers (Figs. 1B and S 3 ) 39 , 40 . Fig. 1 Single-cell RNA sequencing reveals dynamic changes in non-neural cell populations during development and regeneration. A Experimental design of single-cell RNA sequencing experiment, starting with the tail tip (stage 41, posterior third of the tail) and mid-trunk (100 µm proximal tissue to amputation plane) before following up with developmental uninjured (un) or amputated samples at 1, 3, or 7 days post amputation (dpa). Samples were collected and sorted via FACS before processing and sequencing. B UMAP of 35,852 cells grouped into 21 neural and non-neural cell clusters. C UMAP of each timepoint, including developmental and regenerating samples. D Stacked bar charts showing the cluster makeup of each sample. + denotes expanded goblet cells and basal cells in development, ++ denotes expanded immune cells in regeneration, and +++ denotes expanded syndetome and sclerotome in regeneration. Colors consistently identify unique cell types across ( B , C ).
A Experimental design of single-cell RNA sequencing experiment, starting with the tail tip (stage 41, posterior third of the tail) and mid-trunk (100 µm proximal tissue to amputation plane) before following up with developmental uninjured (un) or amputated samples at 1, 3, or 7 days post amputation (dpa). Samples were collected and sorted via FACS before processing and sequencing. B UMAP of 35,852 cells grouped into 21 neural and non-neural cell clusters. C UMAP of each timepoint, including developmental and regenerating samples. D Stacked bar charts showing the cluster makeup of each sample. + denotes expanded goblet cells and basal cells in development, ++ denotes expanded immune cells in regeneration, and +++ denotes expanded syndetome and sclerotome in regeneration. Colors consistently identify unique cell types across ( B , C ).
We first examined global trends in cell type composition as tail development and regeneration progressed. Many of the patterns we observe for days 1–3 are similar to those observed in X. laevis 39 . While the overall cell type composition of tail tips and mid-trunk tissue was similar, we noted several striking changes in cell type abundance between the stage-matched developmental and regenerating timepoints by 7 dpa (Fig. 1C ). Examples include the expansion of goblet cells and basal cells in development but not regeneration (Fig. 1D + ), and the relative enrichment during regeneration of immune cells (Fig. 1D ++ ), which are known to be critical for regenerative success 41 , 42 , and immature connective tissue types like sclerotome and syndetome (Fig. 1D +++ ) during regeneration.
For the remainder of our study, we focused on changes in neural cell types. To interrogate spinal cord regeneration, we subset 7337 neural cells (3708 uninjured cells and 3629 regenerating cells). We then set out to use molecular markers to better define a molecular atlas of neuron types. Previous designations of larval neuron types in Xenopus have relied on their morphology and reactivity. These approaches have identified dorsal Rohon-Beard neurons (RBNs) 15 , intermediate excitatory and inhibitory INs, ventrolateral MNs, and ventral cerebrospinal fluid contacting Kolmer–Agduhr interneurons (KAs) 16 , 17 (Fig. 2A , reviewed in refs. 18 – 20 ). Three subtypes of INs have subsequently been identified based on their expression of vsx2, en1 , or evx1 16 , 43 , 44 . In our analysis, we first mapped clusters back to Xenopus neuron types based on known neurotransmitter types and markers, where available. We then used unique patterns of gene expression from zebrafish and mouse to assign identities to neuron classes that had not previously been molecularly classified (Fig. 2B, C and Table S1 ) 45 . For all neurons, the neurotransmitter phenotype was used in parallel with established markers in order to avoid the possibility of neurotransmitter plasticity post injury 46 . Fig. 2 Distinct major and cardinal neurons can be identified using conserved markers in the Xenopus tropicalis spinal cord. A Cross-section diagram showing the relative position of major spinal cord cell types in a stage 41 tadpole. B UMAP of the merged neural subset showing major NPC and neuron cell types. C Top 3 markers of each cell cluster. D UMAP and legend of subset cardinal interneuron types in dorsal/ventral order. In addition to differentiating interneurons (diff. INs), interneuron populations include KA’ and KA” clusters, a cluster corresponding to V2b inhibitory interneurons (Ventral Longitudinal Descending or VeLD in zebrafish), a cluster corresponding to V2a Excitatory interneurons (Circumferential Descending (CiD) interneurons in zebrafish), a cluster corresponding to V1 inhibitory interneurons (Circumferential Ascending (CiA) interneurons in zebrafish), a cluster corresponding to V0v excitatory interneurons (Multipolar Commissural Descending and Unipolar Commissural Descending or MCoD/UCoD in zebrafish), and a cluster corresponding to inhibitory dI6 interneurons (Commissural Local or CoLo in zebrafish). These types were determined by ( E ) the top 3 gene markers of each interneuron type and marker genes (Table S1 ). X denotes unidentifiable glycinergic interneuron groups. F Dot plot showing gene expression between KA”, intermediate KA (Int KA), and KA’ interneuron groups.
A Cross-section diagram showing the relative position of major spinal cord cell types in a stage 41 tadpole. B UMAP of the merged neural subset showing major NPC and neuron cell types. C Top 3 markers of each cell cluster. D UMAP and legend of subset cardinal interneuron types in dorsal/ventral order. In addition to differentiating interneurons (diff. INs), interneuron populations include KA’ and KA” clusters, a cluster corresponding to V2b inhibitory interneurons (Ventral Longitudinal Descending or VeLD in zebrafish), a cluster corresponding to V2a Excitatory interneurons (Circumferential Descending (CiD) interneurons in zebrafish), a cluster corresponding to V1 inhibitory interneurons (Circumferential Ascending (CiA) interneurons in zebrafish), a cluster corresponding to V0v excitatory interneurons (Multipolar Commissural Descending and Unipolar Commissural Descending or MCoD/UCoD in zebrafish), and a cluster corresponding to inhibitory dI6 interneurons (Commissural Local or CoLo in zebrafish). These types were determined by ( E ) the top 3 gene markers of each interneuron type and marker genes (Table S1 ). X denotes unidentifiable glycinergic interneuron groups. F Dot plot showing gene expression between KA”, intermediate KA (Int KA), and KA’ interneuron groups.
Through this combined approach of validation and conservation, we identified RBNs as well as several types of both known and new MNs, INs, and KAs in our merged dataset. We found two main subtypes of MNs in our dataset, a division that has not been previously reported: a prdm14+ cluster, and prdm14 − cluster (Fig. 2B, C ). Within the IN cluster, we confirmed the presence of en1 + inhibitory GABAergic/glycinergic INs, vsx2 + excitatory glutamatergic/cholinergic INs, and evx1 + excitatory glutamatergic INs as previously identified in Xenopus laevis 17 , 43 , 44 . These populations likely correspond with V1, V2a, and V0v INs, which have been morphologically identified as ascending or descending INs and commissural INs in Xenopus laevis , respectively (Fig. 2D, E and Table S1 ). In addition, we found populations of descending INs corresponding to V2b inhibitory INs and dI6 INs that likely represent commissural INs (Fig. 2D, E and Table S1 ). While other cardinal classes, such as V3 and dI1-5, could not be identified, they were identified in a recent preprint for X. laevis 47 , suggesting they may be present at low abundance at these stages and conditions, or in different stages or more anterior positions. Finally, one cluster was identified as differentiating INs and one cluster of glycinergic INs remained unidentifiable (Fig. 2D, E , “X”).
Interestingly, the inhibitory KA INs split into 2 mature clusters and one intermediate cluster. Two subtypes of KAs, KA’ and KA”, have been reported in zebrafish, mice, and other species but not previously in Xenopus 48 , 49 . In our dataset, both clusters express gad2 and pkd1l2 , but only one expresses KA” markers foxa2 , nkx6-2 , and th ( tyrosine hydroxylase ). Given that the th expressing cluster did not express other markers of catecholaminergic or serotonergic activity ( slc6a2, slc6a4, tph2 ) and did express slc6a3 , we provisionally identified these as GABAergic/dopaminergic KA” INs (Fig. 2F ).
Next, we examined how neuron populations changed over development and regeneration. First, we found that all major neuron types (RBNs, INs, and both MNs) identified in uninjured populations were present during regeneration in our single-cell dataset, indicating that neuronal diversity is largely regenerated after amputation (Figs. 3A and S 4A ). While this confirms the presence of these neurons in the regenerating spinal cord, we next wanted to investigate if the spatial neuron organization found in uninjured tissue also regenerates. To identify the position of each cell, we performed Hybridization Chain Reaction (HCR) in situ RNA detection using marker genes for major cell types ( chat −/ prdm14 +: RBNs, lhx1 +: INs, chat +/ prdm14 −: MNs, chat +/ prdm14 +: MNs) (Figs. 3B and S 4B ). We then mapped the XY location of each identified neuron, based on HCR punctae and DAPI nuclear staining, normalizing coordinates to control for spinal cord height and width (Figs. 3B, C and S 4C ). Data from multiple samples were integrated using hex plots to track the relative position of all cells across the spinal cord (Figs. 3C, D and S 4F ). We then compared the positional distribution of each cell type over developmental and regenerative timepoints (Figs. 3E and S 4C–E ). Fig. 3 Major neuron types can be found in spatial domains within regenerating tissue. A UMAP showing cell clusters during developmental and regenerating timepoints. B Max projection of representative HCR-labelled spinal cord cross-sections with prdm14 in yellow, lhx1 in green, and chat in magenta. Scale bar is 10 µm. C Normalized hex charts of the representative cross-sections shown in ( B ), with prdm14 + RBNs in yellow, lhx1 + INs in green, chat + /prdm14 − MNs in magenta, and chat +/ prdm14 + MNs in cyan. D Summed and normalized hex charts of spinal cord cross-sections at each timepoint. Color brightness indicates cell count. E Boxplot and jittered dotplot showing normalized distance from the floor plate for each cell. Across ( B – E ): St. 41/uninj n = 98 cells, 4 tadpoles; un3 dpa n = 131 cells, 4 tadpoles; 3 dpa n = 123 cells, 6 tadpoles; un7 dpa n = 330 cells, 3 tadpoles; 7 dpa n = 330 cells, 3 tadpoles. * p -value < 0.05 and *** p -value < 0.001 by a two-sided Wilcox test of means, + p -value < 0.05 by a two-sided Levene’s test of variance equality. Exact p -values provided in Fig. S3D, F . Data are presented as median values +/− Q1/Q3, with minimum and maximum. For all panels, white solid line denotes spinal cord boundary (outer) and central canal (inner). Source Data are provided as a Source Data file.
A UMAP showing cell clusters during developmental and regenerating timepoints. B Max projection of representative HCR-labelled spinal cord cross-sections with prdm14 in yellow, lhx1 in green, and chat in magenta. Scale bar is 10 µm. C Normalized hex charts of the representative cross-sections shown in ( B ), with prdm14 + RBNs in yellow, lhx1 + INs in green, chat + /prdm14 − MNs in magenta, and chat +/ prdm14 + MNs in cyan. D Summed and normalized hex charts of spinal cord cross-sections at each timepoint. Color brightness indicates cell count. E Boxplot and jittered dotplot showing normalized distance from the floor plate for each cell. Across ( B – E ): St. 41/uninj n = 98 cells, 4 tadpoles; un3 dpa n = 131 cells, 4 tadpoles; 3 dpa n = 123 cells, 6 tadpoles; un7 dpa n = 330 cells, 3 tadpoles; 7 dpa n = 330 cells, 3 tadpoles. * p -value < 0.05 and *** p -value < 0.001 by a two-sided Wilcox test of means, + p -value < 0.05 by a two-sided Levene’s test of variance equality. Exact p -values provided in Fig. S3D, F . Data are presented as median values +/− Q1/Q3, with minimum and maximum. For all panels, white solid line denotes spinal cord boundary (outer) and central canal (inner). Source Data are provided as a Source Data file.
In the uninjured samples, all cell types were identified and conformed to their expected dorsal/ventral positions, with dorsal RBNs, intermediate INs, and ventral MNs (Fig. 3D, E ). During regeneration, these positional identities were mostly re-established. There were no significant changes in mean position for any of these cell types at 3 dpa. However, the position of both prdm14 − and prdm14 + MNs shifted significantly ventral by 7 dpa relative to their uninjured counterparts, while INs increased significantly in their range of distribution (Figs. 3E and S 4C–E ). Overall, we conclude that positional identity across the dorsal-ventral axis is largely re-established during regeneration, with ventral positioning becoming more exaggerated for MNs.
After confirming the regeneration of major neuron types and their spatial organization, we next wanted to see if neurons made at different times in development regenerate to the same extent in our scRNA-seq dataset. First, we categorized neurons by their developmental birth window: “early” to indicate neurons established by late neurula (primary neurogenesis), “mid” for neurons arising after neurulation but no longer being generated by stage 41 (secondary neurogenesis), and “late” for neurons still undergoing neurogenesis by stage 41 (also secondary neurogenesis) (Fig. 4A ). We assigned a preliminary primary or secondary origin to neuron categories based on what is known from X. laevis and zebrafish (Fig. 4B and Table S1 ) 20 , 50 – 53 , then confirmed this identity based on the developmental population dynamic data (Fig. 4C ). We found that during our developmental timepoints, primary early neurons decreased in proportion (RBNs, V2b INs), while secondary or late neurons increased in proportion beginning either at day 0 (MNs, V0v INs, KA’s) or day 1 (V2a INs, V1 INs) (Fig. 4C ). We identified “secondary mid neurons” (KA”s, dI6 INs, cerulean and yellow in Fig. 4B, C ) based on previous evidence indicating a secondary origin 15 , 16 but a decrease in proportion in our dataset. Interestingly, in comparison to KA”, KA’ neurons increase during our sampled developmental stages, illustrating a previously unidentified divergence between the two cell types (Fig. 4C ). Fig. 4 Neuron subtype repopulation after injury depends on developmental identity. A Schematic showing waves of primary and secondary neurogenesis ( y axis) over developmental time (developmental stages over x axis), additionally annotated as early neurogenesis (occurring during neurulation), mid neurogenesis (occurring between neurulation and stage 41), and late neurogenesis (occurring during stage 41) 95 . B Literature analysis from Xenopus laevis or Danio rerio (Table S1 ). C Proportional bar charts showing changes in the % neural population over developmental time, as labeled as decreased (−) or increased (+). D Proportional bar charts per neuron type showing changes in the % neural population over regenerative time, as labeled as decreased (−) or increased (+). Arrow denotes expanded prdm14− MN cluster at 7 dpa. E Stacked bar charts showing the proportion of neurons at each time point. Arrow denotes expanded prdm14− MN cluster at 7 dpa.
A Schematic showing waves of primary and secondary neurogenesis ( y axis) over developmental time (developmental stages over x axis), additionally annotated as early neurogenesis (occurring during neurulation), mid neurogenesis (occurring between neurulation and stage 41), and late neurogenesis (occurring during stage 41) 95 . B Literature analysis from Xenopus laevis or Danio rerio (Table S1 ). C Proportional bar charts showing changes in the % neural population over developmental time, as labeled as decreased (−) or increased (+). D Proportional bar charts per neuron type showing changes in the % neural population over regenerative time, as labeled as decreased (−) or increased (+). Arrow denotes expanded prdm14− MN cluster at 7 dpa. E Stacked bar charts showing the proportion of neurons at each time point. Arrow denotes expanded prdm14− MN cluster at 7 dpa.
If spinal cord regeneration fully recapitulates development, we would expect all neuron types to be newly generated in parallel from NPCs. Instead, neuron populations changed dynamically during regeneration based on their window of neurogenesis and identity. Most “early” and “mid” neurons, which have exited their neurogenic window by the time of amputation, proportionally decreased after injury. Most late neurons, which are still being generated at stage 41, proportionally increased during at least one regenerating timepoint (Fig. 4D ). However, several of these late neuron types diverged from their developmental trajectory over regenerative time. The most striking were the late prdm14 + MNs, V0v INs, and V2A INs, which expanded immediately after injury but then decreased markedly at later timepoints, and the late KA” neurons, which declined steadily over regenerative time (Fig. 4D ). Only two neuron types ( prdm14 − MNs and V1 INs) increased in proportion from 3 to 7 dpa. These prdm14− MNs undergo a unique and dynamic pattern: first expanding by 1 dpa, then declining at 3 dpa, and expanding again by 7 dpa (Fig. 4D , arrow), such that by 7 dpa, the large majority of neurons are prdm14 − MNs (Fig. 4E ).
The dynamic and cell-type specific changes in neuronal subtypes we observed suggested that waves of neurogenesis might be occurring between 0–1 dpa and 3–7 dpa, and so we next set out to investigate if these could be attributed to NPC proliferation. To assess the balance of progenitors and post-mitotic cells, we calculated the percent population of cells in S, G2M, or G1 phase using Seurat-based cell cycle prediction analysis. During development, the proportion of neural cells in G2M declines steadily from 0 to 7 days (Figs. 5A , left panel and S5A ). This is accompanied by a decrease in NPCs and in cell division score (GO:0051301) 54 , 55 (Fig. 5B, C , left panels and Fig. S 5A, B, D, E ). By contrast, in regeneration there is a marked increase in G2M at 3 dpa accompanied by an expansion of NPCs and an increase in the cell division score at the same timepoint (Fig. 5B, C , right panels and Fig. S 5B, D, E ). At 1 dpa, the proportion of NPCs actually decreases, as does the proportion of cells in G2M, in agreement with our previous findings (Figs. 5A, B and S 5A ) 33 . Neuron differentiation score (GO:0030182) 54 , 55 was only enriched at 7 dpa, not 1 dpa (Figs. 5D and S 5C, D, F ). This was matched by the upregulation of classic neuron differentiation genes ( ascl1, neurog2, neurog1 ) 56 at 7 dpa (Fig. S 5G ). Fig. 5 Regeneration proceeds through a wave of proliferation and then proliferative neurogenesis. A Bar chart of G1, G2M, and S phase cells as a proportion of entire neural population. B Bar chart of NPCs, Differentiating (Diff.) Neurons, and Neurons as a proportion of entire neural population. C GO:0051301 cell division score over developmental and regenerative time. D GO:0030182 neuron differentiation score over development and regenerative time. For ( C and D ) **** p -value < 0.0001, ** p -value < 0.01 compared to tail tip for developmental timepoints and mid-trunk for regenerating timepoints via a two-sided pairwise Wilcox Test. Exact p -values in Fig. S5D . Data are presented as median values +/− Q1/Q3, with minimum and maximum. E Dotplot and boxplot of the average number of HuC/HuD+ neurons in regenerated tissue over regeneration. * p -value 3–7 dpa = 0.025, 3–8 dpa = 0.020, via a two-sided Wilcox Test; 1 dpa n = 4, 2 dpa n = 9, 3 dpa n = 10, 4 dpa n = 7, 5 dpa n = 6, 6 dpa n = 10, 7 dpa n = 9. F Experimental diagram of BrdU/HuC/HuD Immunohistochemistry. Data are presented as median values +/− Q1/Q3, with minimum and maximum. G Whole mount image showing double HuC/HuD and BrdU+ cell at 6 dpa after injection at 4 dpa. White dashed lines indicate spinal cord border; magenta dotted lines indicate HuC/HuD+ neurons; green dotted lines indicate BrdU+ cells; yellow dotted line indicates double positive cell. Scale bar = 5 µm. H Dotplot and boxplot showing the percent of BrdU/HuC/HuD+ cells within HuC/HuD+ cells 48 h post injection. ** p -value = 0.003, *** p -value = 0.0003, via a two-sided Wilcox Test; for ( G , H ): 0–3 dpa n = 9, 1–3 dpa n = 7, 2–4 dpa n = 6, 3–5 dpa n = 5, 4–6 dpa n = 7, 5–7 dpa n = 6, 6–8 dpa n = 7. Data are presented as median values +/− Q1/Q3, with minimum and maximum. I Experimental diagram showing BrdU injection, sample collection, and iterative co-staining with HCR and BrdU. J Cross-section image showing chat +/ prdm14− BrdU+ MN at 6 dpa after injection at 4 dpa. Sample was imaged 1st with HCR (left), 2nd with BrdU (middle), then aligned to identify double-positive cells (right). White dashed lines indicate spinal cord border; green dotted lines indicate BrdU+ Neurons; magenta dotted lines indicate chat +/ prdm14 − cells; yellow dotted line indicates double positive cell. Scale bar = 5 µm. n = 2 chat +/BrdU+ cells (6% of BrdU+ cells; 5% of neurons) from 9 sections generated from 2 tadpoles. Total BrdU+ cells = 36; total neurons = 39; overall chat +/ prdm14 − cells = 30, prdm14 + cells = 1, chat +/ prdm14 + cells = 4, lhx1 + INs. All statistical tests done via Wilcoxon test. Source Data are provided as a Source Data file.
A Bar chart of G1, G2M, and S phase cells as a proportion of entire neural population. B Bar chart of NPCs, Differentiating (Diff.) Neurons, and Neurons as a proportion of entire neural population. C GO:0051301 cell division score over developmental and regenerative time. D GO:0030182 neuron differentiation score over development and regenerative time. For ( C and D ) **** p -value < 0.0001, ** p -value < 0.01 compared to tail tip for developmental timepoints and mid-trunk for regenerating timepoints via a two-sided pairwise Wilcox Test. Exact p -values in Fig. S5D . Data are presented as median values +/− Q1/Q3, with minimum and maximum. E Dotplot and boxplot of the average number of HuC/HuD+ neurons in regenerated tissue over regeneration. * p -value 3–7 dpa = 0.025, 3–8 dpa = 0.020, via a two-sided Wilcox Test; 1 dpa n = 4, 2 dpa n = 9, 3 dpa n = 10, 4 dpa n = 7, 5 dpa n = 6, 6 dpa n = 10, 7 dpa n = 9. F Experimental diagram of BrdU/HuC/HuD Immunohistochemistry. Data are presented as median values +/− Q1/Q3, with minimum and maximum. G Whole mount image showing double HuC/HuD and BrdU+ cell at 6 dpa after injection at 4 dpa. White dashed lines indicate spinal cord border; magenta dotted lines indicate HuC/HuD+ neurons; green dotted lines indicate BrdU+ cells; yellow dotted line indicates double positive cell. Scale bar = 5 µm. H Dotplot and boxplot showing the percent of BrdU/HuC/HuD+ cells within HuC/HuD+ cells 48 h post injection. ** p -value = 0.003, *** p -value = 0.0003, via a two-sided Wilcox Test; for ( G , H ): 0–3 dpa n = 9, 1–3 dpa n = 7, 2–4 dpa n = 6, 3–5 dpa n = 5, 4–6 dpa n = 7, 5–7 dpa n = 6, 6–8 dpa n = 7. Data are presented as median values +/− Q1/Q3, with minimum and maximum. I Experimental diagram showing BrdU injection, sample collection, and iterative co-staining with HCR and BrdU. J Cross-section image showing chat +/ prdm14− BrdU+ MN at 6 dpa after injection at 4 dpa. Sample was imaged 1st with HCR (left), 2nd with BrdU (middle), then aligned to identify double-positive cells (right). White dashed lines indicate spinal cord border; green dotted lines indicate BrdU+ Neurons; magenta dotted lines indicate chat +/ prdm14 − cells; yellow dotted line indicates double positive cell. Scale bar = 5 µm. n = 2 chat +/BrdU+ cells (6% of BrdU+ cells; 5% of neurons) from 9 sections generated from 2 tadpoles. Total BrdU+ cells = 36; total neurons = 39; overall chat +/ prdm14 − cells = 30, prdm14 + cells = 1, chat +/ prdm14 + cells = 4, lhx1 + INs. All statistical tests done via Wilcoxon test. Source Data are provided as a Source Data file.
This analysis suggested that only the expansion of neurons we see from 3 to 7 dpa (which included V1 INs and prdm14 − MNs) could be explained by proliferative neurogenesis from NPCs. To track neurogenesis in vivo, we stained for the pan-neuronal marker HuC/HuD at each day of regeneration. We found that the average number of neurons per regenerate is static from 1 to 5 dpa, then increases from 5 to 7 dpa (Fig. 5E ). To test if this increase from 5 to 7 dpa is due to proliferative neurogenesis, we injected tadpoles with BrdU according to established protocols, using a 2-day chase (Fig. 5F ) 57 . We found that little to no proliferative neurogenesis in the regenerated tissue prior to 3 dpa. Instead, there is a peak in neurogenesis at 4–6 dpa, with ~12% BrdU+/HuC/D+ cells, before declining to baseline levels (Figs. 5G, H and S 6 ). We conclude that proliferative neurogenesis begins between 3 and 4 dpa, in agreement with our scRNA-seq data (Fig. 5A–C ). This likely fuels the formation of new neurons we observe 5–7 dpa (Fig. 5E ).
To test which neuron types were generated by proliferative neurogenesis, we did iterative staining of HCR and BrdU (Fig. 5I ). The only double positive cells we found at 7 dpa were BrdU+ chat +/ prdm14 − cells, indicating that the most common BrdU+ HuC/HuD+ neurons generated 5–7 dpa are either V2a INs or prdm14 − MNs (Fig. 5J ), echoing the strong increase in prdm14 − MNs at 7 dpa compared to most other neuron types (Fig. 4 ). However, proliferative neurogenesis cannot explain the increase in neurons observed at 1 dpa (Fig. 5B ) or the increase in V0v, V2a, prdm14+ and prdm14− MNs at that timepoint (Fig. 4D ).
In zebrafish, a transient group of injury-induced neurons (iNeurons) arises soon after spinal cord injury and plays a central role in recovery. iNeuron markers were shown to colocalize with various neuron markers such as hb9, isl1, pax2, gad1b , and vgluta , indicating they arise from various neuron types after injury 32 . Given the lack of both proliferation and neurogenesis at 1 dpa as well as previous evidence of leptin ( lep )+ and leptin receptor ( lepr )+ regeneration-specific neurons in Xenopus 33 , 39 , we hypothesized that there might be a similar population of injury-responsive neurons at this timepoint.
To test this hypothesis, we identified subclusters of neurons that appeared transiently after injury (Fig. 6A–C ). These transient subclusters arose in V2a INs, prdm14 + MNs, prdm14 − MNs, V0v INs, and RBNs primarily at 1 and 3 dpa (Fig. 6B, C ). Despite this temporal commonality, each subcluster showed unique patterns of gene expression (Fig. S7 ). The V2a IN, prdm14 − MN, and prdm14 + MN subclusters all strongly expressed atf3 , a gene previously implicated in neuron regeneration and a marker of zebrafish iNeurons (Fig. 6D ) 32 , 58 . A significant difference, however, is that our clusters retain greater transcriptional similarity to their neuron type of origin and are distributed in several clusters, whereas iNeurons form one cluster after injury. The atf3 + MN subclusters also express lep and lepr , as previously identified (Fig. 6D ) 33 , 39 , 59 . Re-analysis of the zebrafish iNeuron dataset identified that, though most top marker genes are not shared between the two species (Fig. S 8A–C ), top signaling pathways share several commonalities (e.g., L1CAM signaling, NCAM, NRXN, NRG, SEMA3, Fig. S 8F–K ), and lepr is also expressed in the iNeuron cluster at 1 week post injury (Fig. S 8D, E ), indicating that zebrafish iNeurons and Xenopus lepr + transient MNs may share conserved regeneration-specific signaling features. Finally, we used HCR to track expression of lepr , atf3 , and the MN marker chat at 3 dpa, and saw localization of both lepr + /atf3 + /chat+ and lepr + /atf3 + /chat− cells anterior and posterior to the amputation plane in regenerating animals (Fig. S 9A–C ), confirming this is a broadly-activated injury response in both MN and non-MN neurons in vivo. Fig. 6 Regeneration-specific neuron subclusters arise in specific neuron types and show distinct patterns of gene expression and cell processes. A Combined UMAP showing sample distribution; black circles show regeneration-specific subclusters. B , C Regeneration specific subclusters (black cells and circles) appear near V2a INs, prdm14 + MNs, prdm14− MNs, V0v INs, and RBNs by 1 dpa, but not near Neural Progenitor Cells (NPCs)/ diff N. (Differentiated Neurons), KA’, KA”, V2b, dI6, V1, or unknown (X) groups, as shown by combined UMAP ( B ) and proportional bar charts ( C ). D Gene expression of NPCs, Differentiating (Diff) Neurons, Other Neuron types, and regeneration-specific subclusters. Some subclusters express iNeuron markers ( atf3, gap43 ), others express Xenopus lep + MN markers ( lepr, lep ), and others have unique expression patterns ( sfrp4, rtn4rl1, ppp1r17 ). Arrows highlight atf3 + Regenerating Neurons. Y -axis dots represent clusters gene is expressed. E Gene ontology of atf3 + Regenerating Neurons compared to other Neuron types at 1 dpa. Enriched terms are sorted by strength of adjusted p -value, which were calculated and corrected for multiple testing using “enrichGO” from clusterProfiler. F Boxplot of GO: 0007411 axon guidance score at 1 dpa in different cell clusters in uninjured (left bar) or regenerating (right bar) timepoints. Black denotes regeneration-specific subclusters. * p -value = 0.028, ** p -value = 0.004, *** p -value = 0.00026, **** p -value = 1.8 × 10 −15 via pairwise Wilcox test. Data are presented as median values +/− Q1/Q3, with minimum and maximum. G Plot of GO:0007411 axon guidance score and associated modules at 1 dpa (bottom), compared to upregulated genes representing scores. All statistical tests done via Wilcoxon test.
A Combined UMAP showing sample distribution; black circles show regeneration-specific subclusters. B , C Regeneration specific subclusters (black cells and circles) appear near V2a INs, prdm14 + MNs, prdm14− MNs, V0v INs, and RBNs by 1 dpa, but not near Neural Progenitor Cells (NPCs)/ diff N. (Differentiated Neurons), KA’, KA”, V2b, dI6, V1, or unknown (X) groups, as shown by combined UMAP ( B ) and proportional bar charts ( C ). D Gene expression of NPCs, Differentiating (Diff) Neurons, Other Neuron types, and regeneration-specific subclusters. Some subclusters express iNeuron markers ( atf3, gap43 ), others express Xenopus lep + MN markers ( lepr, lep ), and others have unique expression patterns ( sfrp4, rtn4rl1, ppp1r17 ). Arrows highlight atf3 + Regenerating Neurons. Y -axis dots represent clusters gene is expressed. E Gene ontology of atf3 + Regenerating Neurons compared to other Neuron types at 1 dpa. Enriched terms are sorted by strength of adjusted p -value, which were calculated and corrected for multiple testing using “enrichGO” from clusterProfiler. F Boxplot of GO: 0007411 axon guidance score at 1 dpa in different cell clusters in uninjured (left bar) or regenerating (right bar) timepoints. Black denotes regeneration-specific subclusters. * p -value = 0.028, ** p -value = 0.004, *** p -value = 0.00026, **** p -value = 1.8 × 10 −15 via pairwise Wilcox test. Data are presented as median values +/− Q1/Q3, with minimum and maximum. G Plot of GO:0007411 axon guidance score and associated modules at 1 dpa (bottom), compared to upregulated genes representing scores. All statistical tests done via Wilcoxon test.
Previously, zebrafish iNeurons were found to have a role in cell signaling and neural plasticity through cell mechanisms like axon regrowth and migration 32 . To assess if atf3 + subclusters play a similar central signaling role in Xenopus , we investigated what GO terms were enriched in atf3 + regenerating neurons compared to non-regeneration specific neurons at 1 dpa. Similar to zebrafish iNeurons, the top-enriched GO term was axon guidance, with other top terms including intracellular signaling transduction, neuron projection development (also increased in zebrafish iNeurons), and nervous system development (Figs. 6E, F and S 10A–E ). To see what genes were driving this enrichment, we plotted genes associated with nervous system development against axon guidance and other top overlapping associated GO terms. Many nervous system development genes were associated with axon guidance and specifically the semaphorin-plexin signaling pathway and/or cell-cell adhesion (Fig. 6G ), suggesting that like iNeurons, atf3 + cells prioritize neural plasticity via axon guidance/neurite outgrowth during early regeneration (Fig. S 10D, E ).
Though all neuron types exist in regenerated tissue, proliferative neurogenesis doesn’t occur until 4 dpa (Fig. 5 ) and is absent for some neuron types, notably Kolmer–Aghdur neurons (Fig. 4 ). This raises the question of how these neurons repopulate. A possible mechanism is movement of post-mitotic, differentiated neurons from uninjured tissue into new tissue, as has been reported in urodeles 29 . We also noted that cell migration was a significantly enriched GO term in all KA populations at 1 dpa (Fig. S 11A ). We queried the contributions of cell division and cell migration to KA” neurons, using Tyrosine hydroxylase (TH) as a unique marker of this cell type (Fig. 2F ). First, we asked if the number of KA” neurons present at amputation increases over regeneration. We found that KA” do not increase significantly in regenerating tails (Fig. S 11B ), suggesting that the pool of KA” neurons present at injury is now just spread over the longer distance of the regenerated tail. We also found zero instances of Brdu+ KA neurons in the regenerate ( n = 503 BrDU+ cells) (Fig. 7A ). We then used Hydroxyurea (HUA) and Aphidicolin (Aph) 60 to inhibit all cell proliferation, first confirming that this treatment during the first 72 h caused a full regeneration arrest, as would be expected for blocking all cell proliferation (Fig. S 11C ). We then treated regenerating tails with HUA/Aph only during the 3–6 dpa window when proliferative neurogenesis is active (Fig. 5 ). We found that the total number of neurons in the treated regenerating tail fell by 82% (Fig. 7B–E , magenta cells in E compared to D). However, the number of KA neurons was unaffected by HUA/Aph treatment (Fig. 7B , white cells in E compared to D), making them now the most abundant neuron type in the regenerated tail (Fig. 7C ). These results suggest KA neurons do not rely on proliferation to repopulate the tail. Fig. 7 Post-mitotic neurons are displaced from uninjured tissue to regenerating tissue. A Image of 6 dpa tadpole co-stained for TH (green) and BrdU (magenta). Green dotted lines represent TH + KA” INs, white solid lines show spinal cord border. Number of positively-labeled cells indicated in bottom right corner, collected from 2 tadpoles. B Boxplot and dotplot of the total number of HuC/D+ neurons or TH+ neurons in the regenerated spinal cord of 6 dpa tadpoles treated from 3 to 6 dpa with either 0.1% DMSO (control, n = 6 for HuC/D; and n = 3 for TH) or 20 mM Hydroxyurea and 150uM Aphidicolin (HUA/Aph) ( n = 10 HuC/D; n = 4 tadpoles for TH). *** p -value = 0.000002 by two sided t-test Data are presented as median values +/− Q1/Q3, with minimum and maximum. C Boxplot for the percentage of TH+ neurons as a fraction of all HuC/D+ neurons at 6 dpa ( n = 3 for control/DMSO, n = 4 tadpoles for HUA/Aph). ** p -value = 0.003 by two-sided t -test. Data are presented as median values +/− Q1/Q3, with minimum and maximum. D , E Images of 6 dpa tails stained for TH (green) and HuC/D (magenta) after treatment from 3 to 6 dpa with DMSO ( D ) or HUA/Aph ( E ). Note the abundant single stained HuC/D+ magenta cells in DMSO (magenta arrowheads) relative to double-stained HuC/D + /TH+ white cells (white arrowheads), with much fewer magenta cells in HUA/Aph treatment (magenta asterisk). F Boxplot and dotplot (left) and images (right) of TH + KA” cells in the regenerating tissue in 0.1% or 5 µM DCLK inhibitor. *** p -value = 0.0004 by two-sided t -test. DMSO n = 10, DCLK inh n = 14 tadpoles. Data are presented as median values +/− Q1/Q3, with minimum and maximum. G Top: Diagram of experimental design: Tol2 HuC:Dendra2 plasmid and tol2 RNA was injected into the single-celled embryo. At stage 41, tadpoles were photoconverted and amputated. At 3 dpa, tadpoles were imaged. Bottom: Representative image of double-positive neuron, dendrite, and axon in the regenerated tail. White dotted line represents amputation plane. White box represents insets of the 488 (top, green) and 594 (bottom, magenta) channels. Solid white lines represent spinal cord border. Sample size detailed in Fig. S12 . Source Data are provided as a Source Data file.
A Image of 6 dpa tadpole co-stained for TH (green) and BrdU (magenta). Green dotted lines represent TH + KA” INs, white solid lines show spinal cord border. Number of positively-labeled cells indicated in bottom right corner, collected from 2 tadpoles. B Boxplot and dotplot of the total number of HuC/D+ neurons or TH+ neurons in the regenerated spinal cord of 6 dpa tadpoles treated from 3 to 6 dpa with either 0.1% DMSO (control, n = 6 for HuC/D; and n = 3 for TH) or 20 mM Hydroxyurea and 150uM Aphidicolin (HUA/Aph) ( n = 10 HuC/D; n = 4 tadpoles for TH). *** p -value = 0.000002 by two sided t-test Data are presented as median values +/− Q1/Q3, with minimum and maximum. C Boxplot for the percentage of TH+ neurons as a fraction of all HuC/D+ neurons at 6 dpa ( n = 3 for control/DMSO, n = 4 tadpoles for HUA/Aph). ** p -value = 0.003 by two-sided t -test. Data are presented as median values +/− Q1/Q3, with minimum and maximum. D , E Images of 6 dpa tails stained for TH (green) and HuC/D (magenta) after treatment from 3 to 6 dpa with DMSO ( D ) or HUA/Aph ( E ). Note the abundant single stained HuC/D+ magenta cells in DMSO (magenta arrowheads) relative to double-stained HuC/D + /TH+ white cells (white arrowheads), with much fewer magenta cells in HUA/Aph treatment (magenta asterisk). F Boxplot and dotplot (left) and images (right) of TH + KA” cells in the regenerating tissue in 0.1% or 5 µM DCLK inhibitor. *** p -value = 0.0004 by two-sided t -test. DMSO n = 10, DCLK inh n = 14 tadpoles. Data are presented as median values +/− Q1/Q3, with minimum and maximum. G Top: Diagram of experimental design: Tol2 HuC:Dendra2 plasmid and tol2 RNA was injected into the single-celled embryo. At stage 41, tadpoles were photoconverted and amputated. At 3 dpa, tadpoles were imaged. Bottom: Representative image of double-positive neuron, dendrite, and axon in the regenerated tail. White dotted line represents amputation plane. White box represents insets of the 488 (top, green) and 594 (bottom, magenta) channels. Solid white lines represent spinal cord border. Sample size detailed in Fig. S12 . Source Data are provided as a Source Data file.
We then asked if KA neurons might be repopulating the tail through migration. Here, we treated the regenerating tail with an inhibitor of DCLK signaling, a critical regulator of neuronal migration 61 , 62 . DCLK1 inh1 treatment resulted in a significant decrease in the number of KA neurons in the regenerated tail at 7 dpa (Fig. 7F ), consistent with a possible role for cell migration. Finally, to directly ask whether neurons from the uninjured tissue migrate into the regenerated spinal cord, we injected 1-cell stage embryos with tol2 mRNA and a donor plasmid containing sequence for the photoconvertible protein Dendra2 under the control of the zebrafish HuC promoter 63 . This promoter marks all differentiating and mature neuron types, not only KA neurons, for which we lacked a cell type-specific driver. Dendra2 is subsequently expressed mosaically in a small fraction of HuC+ neurons. At stage 41, the entire tail was photoconverted, amputated, and tadpole allowed to regenerate in the dark (Fig. 7G , S 12A, B ), thus neurons present before the injury would appear double positive (S 12C ). At 3 dpa, we found double-positive Dendra2 cells in the regenerating spinal cord, confirming that neurons move from uninjured tissue into the regenerating spinal cord (Figs. 7G and S 12D, E ). Our results are therefore most consistent with a model in which KA neurons, and possibly other neurons or immature neuroblast types, repopulate the tail by migration during 3–7 dpa, but do not preclude contributions of proliferation or differentiation at other times or stages.
Discussion
Our study creates a molecular atlas of neurons in the Xenopus tadpole, defining marker genes for both known and previously unconfirmed neuronal types. We establish that the overall spatial organization of regeneration largely recapitulates developmental organization, but that the cell type composition of the regenerating spinal cord differs in many respects from its uninjured, developing counterpart. Proliferative neurogenesis occurs both less and later than expected, giving rise principally to prdm14 − MNs. Some neuron types give rise to a transient population of neurons that prioritize axon guidance and neurite outgrowth as early as 1 dpa. Other neuron types, including KA” INs, repopulate the spinal cord without proliferation, implicating post-mitotic movement (Fig. 8 ). These findings have important implications for understanding the mechanisms that may enable plasticity in regenerative conditions, or limit repair in non-regenerative conditions. Fig. 8 Neuron regeneration proceeds dependent on neuron identity. A During development, primary and secondary neurogenesis gives rise to early, mid, and late neuron types, respectively, resulting in a complex population of neurons by stage 41. B At 1 dpa, most late neurons increase in proportion in comparison to early and late neurons. Putative iNeuron clusters arise from prdm14 + MNs, prdm14− MNs, and V2a INs and drive neural plasticity. C At 3 dpa, NPCs prioritize proliferation, although putative iNeuron clusters remain present. Few new neurons are generated, but neurons from uninjured tissue move into regenerating tissue. D New neurons—primarily prdm14− MNs—are generated by proliferative neurogenesis.
A During development, primary and secondary neurogenesis gives rise to early, mid, and late neuron types, respectively, resulting in a complex population of neurons by stage 41. B At 1 dpa, most late neurons increase in proportion in comparison to early and late neurons. Putative iNeuron clusters arise from prdm14 + MNs, prdm14− MNs, and V2a INs and drive neural plasticity. C At 3 dpa, NPCs prioritize proliferation, although putative iNeuron clusters remain present. Few new neurons are generated, but neurons from uninjured tissue move into regenerating tissue. D New neurons—primarily prdm14− MNs—are generated by proliferative neurogenesis.
One of our principal goals was to better index the neuron types that exist in the early Xenopus free-swimming tadpole. By stage 41, tadpole motor circuits for locomotion and escape must be fully functional, but molecular identifiers for all types of neuron populations were sparse 19 . We were able to confirm the presence of and identify markers for each major neuron type and seven conserved cardinal classes of INs that have previously been identified in zebrafish and mouse 45 . In addition, we identified several previously unidentified subtypes; we were especially interested to find that there are two molecularly distinct populations of Kolmer–Agduhr neurons, similar to the KA’ and KA” populations found in zebrafish and mice 16 , 17 , 48 , 49 . We also found two previously unseparated subpopulations of MNs ( prdm14 − and prdm14 +) 64 , 65 . Our study did not identify dl1–5 or V3 INs, though a recent preprint establishes that these do exist in frogs, and are either too low in abundance for us to capture or are only present at later timepoints or more anteriorly 47 . As previously suggested 21 , Rohon-Beard sensory neurons specifically became very rare by later stages of development and regeneration. This raises the question of how sensation is transduced in this transition period before the dorsal root ganglia develop. Functionally, the regenerated tail is still able to respond to light touch by activating an escape reflex. It seems most likely that anterior sensory neurons expand their connectivity to include posterior sensory processes during development and regeneration, though this requires direct evidence through analysis of sensorimotor circuits and behavior, tools whose development we are watching with interest 66 .
A common hypothesis is that regeneration in developing organisms will recapitulate developmental mechanisms to replace lost tissue 67 . Our study supports some aspects of this model. For example, we find that spatial organization of post-mitotic neurons of the spinal cord is broadly restored after injury, a pattern that also agrees with our previous finding that spatially restricted Shh signaling re-patterns the dorsal-ventral axis of the regenerating tadpole spinal cord, in a recapitulation of development 26 . This result recalls the persistent spatial identity that is maintained in mammalian neural progenitors even as they undergo changes in competence 68 , and the similar restoration of positional identity that is governed by Wnt signaling in cnidarians and planarians 69 – 71 . While we didn’t query anterior-posterior patterning in this study, our work and others have made it clear that Wnt, FGF, Hif1a, and Hox genes contribute to neural regeneration of the Xenopus tail 39 , 72 – 74 . An interesting next direction will be to profile how neural cell type composition varies before and after injury along the A/P axis, and how developmental signals like Wnt, which are generally activated by injury but also are posteriorizing, are interpreted at more anterior positions.
In comparison to the perseverance of spatial identity, we find that the composition and dynamics of neuron types in the regenerating spinal cord differ in several crucial aspects from their uninjured, developing counterparts. While all uninjured neuron types are present in the regenerating spinal cord, neurons still undergoing secondary neurogenesis at stage 41 regenerate much more abundantly than those known to be established through primary neurogenesis. Specifically, two mid-secondary neuron types that expand during mid-development (dI6 and KA” INs) instead decline in regeneration, compared to increases in late secondary neurons (V0v, V2a, prdm14 + MNs, and prdm14 − MNs). Interestingly, several secondary neurons that increase throughout development (V0v, V2a, KA’ INs, and prdm14 + MNs) decrease during late regeneration, compensated by exaggerated increases in prdm14 − MNs. Overall, this indicates that regenerative repopulation depends closely on developmental identity: early primary and mid secondary neurons appear no longer able to dynamically repopulate in comparison to late secondary neurons, possibly due to the exhaustion of their unique progenitor pools. This may reflect the ability of the regenerating larval spinal cord to simply sustain secondary neurogenesis to replace some lost neurons, instead of reactivating neurula-stage mechanisms.
The disparate responses to injury by neuronal cell type suggest that neuron types are differentially primed to adapt post-injury, even amongst birthdate groupings. Several possibilities exist for why this diversity exists- it may be that certain neurons are crucial to the early regenerative response, or it may be a secondary effect of identity or location. For example, KA neurons are known to play a role in the development of the Reissner Fiber, a thread of glycoprotein in the central canal that contributes to body axis morphogenesis during development 75 , 76 , spinal curvature detection 77 , and tail straightening during regeneration 78 . If KA”s are unable to regenerate via proliferative neurogenesis due to progenitor limitations, migration of KA” neurons might be necessary to restore guide body axis restoration via the Reissner Fiber during regeneration. On the other hand, KA neurons are also medial-ventrally positioned amongst NPCs next to the central canal. Their migration might be a secondary effect of this unique location. More work is necessary to identify the unique roles different neuron types may play during regeneration, and if this identity is correlated with mechanism.
Another major motivation for our study was to understand how different classes of neurons are restored over the course of regeneration. Previously, we had described that NPCs undergo an early cell cycle exit after injury, followed by later proliferation 33 . Here, the increased duration of our study and the inclusion of stage-matched developmental timepoints allowed us to generate a much clearer picture of how neurogenesis proceeds after injury. In contrast to the axolotl limb, single-cell analysis of the regenerating Xenopus tail does not indicate a clear “blastema” cell identity, nor is there the large mesenchymal mass of cells characteristic of blastema formation found in either the axolotl or the anole 79 – 81 . Our findings bear conceptual similarity to observations of the zebrafish lateral line 82 , in which there are temporally discrete phases of regeneration. The early wave of regeneration-specific neuron production we observe, characterized by atf3 and lepr expression, dominates in V2a INs, prdm14 + MNs, and prdm14 − MNs, but is largely absent in other classes. Meanwhile, the later wave of proliferative neurogenesis contributes strongly to repopulation of prdm14− MNs, but very little to INs or KA” neurons, with the latter not relying on proliferation at all, but rather activating their own program of cell migratory gene expression and movement. We favor active migration as the mechanism enabling KA” neuron repopulation, though they may also be more passively “dragged along” with other cell types. The regeneration-specific differentiation we see is very much like that described for zebrafish iNeurons, with atf3 as a notable shared marker, but unlike zebrafish, Xenopus r egeneration-specific neurons are split into multiple subclusters of both MN groups and V2a INs 32 .
Although our study did not directly interrogate connectivity, similar studies in other regenerative vertebrates, including lamprey and zebrafish, suggest that functional activity and connective plasticity are features of early regeneration, particularly in iNeurons 32 , 83 – 85 . The emphasis we observed on neurite outgrowth, axonogenesis, and axon repair are consistent with similar mechanisms being emphasized in tadpoles. Axon sprouting and neurite outgrowth in mature neurons that are hallmarks of plasticity in fish 32 , 85 , as are functional responses like calcium propagation and muscle activity 85 . Atf3, the marker most strongly associated with axonogenesis and sprouting, is a highly conserved pro-regenerative transcription factor that promotes axon outgrowth and repair, even in mammals 58 . The shared gene expression programs between our regeneration-specific neurons and iNeurons may indicate that these tadpole neurons are also associated with functional plasticity or reinnervation of new tissue, and may point to potential mechanisms underlying why these vertebrates are able to regenerate while others cannot. Loss of these neurons, of migratory capacity, or the coupling of these mechanisms to proliferative neurogenesis, therefore, represent intriguing lines of inquiry for loss of regenerative capacity in adult frogs, as well as potential therapeutic avenues 85 .
Considering the broader implications of our work on the regeneration field, we think our study represents an opportunity to integrate two principal models of neural regeneration: stem cell (or progenitor)-driven proliferative neurogenesis and neuronal-intrinsic activation. In both transection and amputation models of spinal cord regeneration, sox2- expressing neural stem cells (including radial ependymal glia) have long been recognized as a critical source of neurogenesis. These cells are required for regeneration in zebrafish, axolotls, anoles, and Xenopus 86 – 93 . Stem-cell driven neurogenesis is a deeply conserved mechanism for neural regeneration: invertebrates, including planaria 94 , hydra 95 , and acoels 96 , also rely on pluripotent stem cell proliferation to replace neurons. However, recent work has begun to highlight injury-induced activation mechanisms that are intrinsic to injured neurons and may be equally important to regeneration. These include the injury-induced neurons of zebrafish 32 , as well as the long-distance activation of a specific population of dpErk+/ etv + neurons in axolotls 97 . More broadly, long-range activation or priming of regeneration is still an emerging area of regeneration biology even outside of the CNS, with work in axolotls and Xenopus suggesting that injury to one part of the body activates long-range signals that may affect or potentiate regenerative responses in distant tissues 98 , 99 . Even in mammals, a conditioning lesion to peripheral axons of dorsal root ganglion neurons can potentiate regeneration from their central axons 100 , 101 .
Our work suggests that both stem cell-driven neurogenesis and neuronal-intrinsic injury-induced activation are important to regeneration, and that they are partitioned to some extent by cell type. We specifically propose that stem-cell driven proliferation is most relevant for secondary neurogenic cell types, especially MNs, and that these are more easily generated from stem cells than primary neurogenic fates. A similar model governs regenerative outcomes in the retina. In zebrafish, retinal Müller glial cells are able to produce all neuron types after injury 102 , 103 , whereas in the chick and mouse, Müller glial cell potency is limited 104 – 106 . Reprogramming of the Müller glia by expression of Ascl1 and histone modifiers can increase their potency to include a wider range of neuron types, including photoreceptors 107 , 108 . Our model of differential stem cell potency makes several useful parallel predictions. For one, it may be that in animals that retain their regenerative capability throughout life, neural stem cells have wider potency. Within Xenopus , we predict that stem cell potency may vary by axial position and be restricted with age. Finally, reprogramming of ependymal glia, like Müller glia, may represent a path for increasing regenerative capability in non-regenerative animals. Overall, our study highlights that the restoration of neuronal diversity in the spinal cord integrates several cell-type-specific responses, with no single strategy ensuring regeneration of all neurons.