Transcriptomic profiling of innate spinal cord repair in adult zebrafish. A Experimental pipeline to generate a single-cell atlas of zebrafish cells after SCI. SC tissue collection, nuclear isolation and single nuclear RNA-seq (snRNA-seq) were performed on wild-type fish at 1, 3 and 6 wpi. Uninjured SC nuclei were used as 0 wpi controls. 10x genomics sequencing with v3.1 chemistry was performed. Two biological replicates were used at each time point, and SC tissues from 50 animals were pooled into a single biological replicate. B Merged UMAP representation of the complete dataset. Two biological replicates, 4 time points and 58,973 cells were clustered into major spinal cell populations including Neurons A and B, glia/ependymo-radial glial cells (glia/ERGs), oligodendrocyte precursor cells (OPCs), oligodendrocytes (Oligos) A and B, Leukocytes A and B, and Erythrocytes. C Marker gene expression in the major cell types identified after SCI. Dot plot shows canonical markers are enriched in their respective cell types. Dot colors and diameters represent average gene expression and percent cells with at least one UMI detected per gene, respectively. Feature plots depicting the distributions of canonical markers corresponding to the major cell populations are shown. D Distribution of major cell types during SC regeneration. For each time point, cell proportions were normalized to the total number of cells present at that time point. E Heatmap showing top 10 differentially expressed (DE) genes in the major cell populations identified by snRNA-seq. Yellow and black represent highly enriched and low-expression genes, respectively. Source data are provided as a Source Data file.

Cell-cell interaction networks during spinal cord regeneration. A Relative strengths of outgoing signaling pathways in major spinal cell populations at 0, 1, 3 and 6 wpi. Bar graphs at the top of each heatmap show cumulative signaling strengths per cell population. Bar graphs at the right of each heatmap show cumulative signaling strength per pathway. Pathways highlighted in yellow are specifically enriched at one time point. Pathways highlighted in magenta are enriched after injury. B Cumulative strengths of signaling pathways outgoing from or received by major cell types. C Circle plots represent cell-cell interaction strengths of signaling pathways outgoing from Neurons A, B and OPCs at 6 wpi. Arrow thickness is proportional to the cumulative strengths of all predicted interactions between clusters. D Chord diagrams showing intercellular communication of signaling pathways that are enriched in neurons after spinal cord injury (SCI). Arc diameters are proportional to cumulative signaling strengths. Bar graphs at the bottom of each chord diagram represent the relative contribution of different ligand-receptor pairs predicted for each signaling pathway. Source data are provided as a Source Data file.

Recovery of excitatory/inhibitory balance during spinal cord regeneration. A UMAP plot showing 30 neuron subclusters identified with 0.6 resolution parameter. B Classification of neuron clusters based on neurotransmitter properties. Clusters enriched for at least one of the glutamatergic and GABA/glycinergic markers are identified as excitatory and inhibitory, respectively. Clusters that express serotonergic, cholinergic or a mix of glutamatergic and glycinergic markers are labeled as ?Other?. One neuron cluster that does not express any of these marker genes is labeled ?Unidentified?. C Pie charts show the proportions of excitatory, inhibitory and ?Other? neuron populations at 0, 1, 3 and 6 wpi. D in silico calculation of E/I ratios during SC regeneration. E, F Immunostaining for RFP, GFP, HuC/D and Hoechst in SC cross sections from Tg(vglut2a:RFP) and Tg(gad1b:GFP) zebrafish at 0, 1, 3 and 6 wpi. White arrows point to RFP+ HuC/D+ neurons in E and GFP+ HuC/D+ neurons in (F). Insets depict double channels for RFP and HuC/D in E and for GFP and HuC/D in (F). Single channel insets of vglut2a:RFP or gad1b:GFP are shown in dorsal and ventral SC domains, respectively. G, H Quantification of glutamatergic and GABAergic neurons after SCI. Percent vglut2a+ HuC/D+ neurons (ANOVA p: 0.0086) and gad1b+ HuC/D+ neurons (ANOVA p: 0.0043) were normalized to the total number of HuC/D+ neurons. I in vivo quantification of E/I ratios during SC regeneration. Ratios of glutamatergic to GABAergic neurons were calculated (ANOVA p: 0.0268). Data points in all bar charts indicate individual animals and sample sizes are indicated in parentheses. SC cross sections 450 µm rostral to the lesion were analyzed. Brown-Forsythe and Welch ANOVA were performed in (G, H) with Dunnett?s T3 multiple comparisons performed across time points with 95% CI. Ordinary one-way ANOVA with Tukey multiple comparisons was performed in (I). All statistical tests are two-sided. In all bar charts, data are presented as mean values +/- SEM. *p ? 0.05, **p ? 0.01; Scale bars, 50 µm. Source data are provided as a Source Data file.

Regeneration of excitatory and inhibitory neurons after spinal cord injury. A Experimental timeline to assess cell proliferation. SC cross-sections 450 µm rostral to the lesion were used for quantification. B Staining for RFP, GFP, EdU and Hoechst in Tg(vglut2a:RFP; gad1b:GFP) zebrafish at 0, 1, 3 and 6 wpi. Insets marked by dotted lines represent dual channels for RFP and EdU or for GFP and EdU. Yellow arrowheads indicate vglut2a+ EdU+ neurons. Green arrowheads indicate gad1b+ EdU+ neurons. C, D Quantification of vglut2a+ EdU+ neurons (ANOVA p < 0.0001) and gad1b+ EdU+ neurons (ANOVA p: 0.0449) at 0, 1, 3 and 6 wpi. Percent cells were normalized to the total number of nuclei in SC sections. E Experimental timeline to assess the cumulative profiles of regenerating neurons at 1 and 3 wpi. Tissue sections 150, 450 and 750 ?m rostral to the lesion were analyzed. F Staining for RFP, GFP, EdU and Hoechst in Tg(vglut2a:RFP; gad1b:GFP) zebrafish at 0, 1 and 3 wpi. Insets marked by dotted lines represent dual channels for RFP and EdU or for GFP and EdU. Yellow arrowheads indicate vglut2a+ EdU+ neurons. Green arrowheads indicate gad1b+ EdU+ neurons. G, H Quantification of regenerating glutamatergic neurons (vglut2a+ EdU+) (ANOVA p < 0.0001) and regenerating GABAergic neurons (gad1b+ EdU+) (ANOVA p < 0.0001). For each section, percent cells were normalized to the number of vglut2a+ neurons in (G) and to gad1b+ neurons in (H). Data points in all bar charts indicate individual animals and sample sizes are indicated in parentheses. Brown-Forsythe and Welch ANOVA were performed in (C, D) with Dunnett?s T3 multiple comparisons performed across time points with 95% CI. Two-way ANOVA was performed in (G, H) with Tukey?s multiple comparisons with 95% CI. In all bar charts, data are presented as mean values +/- SEM. All statistical tests are two-sided. *p ? 0.05, **p ? 0.01, ***p ? 0.001. Scale bars, 50 µm. Source data are provided as a Source Data file.

Emergence of an injury-responsive iNeuron signature during spinal cord regeneration. A Dynamics of neuronal subclusters at 0, 1, 3 and 6 wpi. Cell proportions for each cluster and time point were normalized to the total number of neurons analyzed at that time point. Clusters that either increased or decreased by >2-fold in proportion are highlighted in red and blue, respectively. B UMAP plot showing the predicted cluster identity for each neuron subpopulations. V0, V2a, V2b, interneurons, motor neurons are shown along with clusters identified as neuroblast-like or progenitor-like neurons. Cluster N20 was the only cluster identified as neuroblast-like neurons. C UMAP plot showing a calculated ?Regeneration Score? for each neuron subclusters. -log10(p value) of each cluster is represented in a gradient scale of gray (min = 0) to red (max = 6). Cluster N20 showed the highest regeneration score compared to all other neurons. D Feature plots show the expression of the regeneration associated genes gap43 and atf3 in cluster N20. E Gene ontology of N20 iNeuron DE markers. Twenty of the most enriched terms are shown. X-axis represents ?log10(p value) for each term. F Circle plots represent cell-cell interaction strengths at 1 wpi. Arrow thickness is directly proportional to the cumulative strengths of all predicted interactions between clusters. G Signaling pathways outgoing from or received by iNeurons at 1 wpi. Bar graphs at the top of each heatmap show cumulative signaling strengths per cell population. Bar graphs at the right of each heatmap show cumulative signaling strength per pathway. One-tailed hypergeometric test is performed in (C, E). Source data are provided as a Source Data file.

in vivo mapping of iNeurons in zebrafish spinal tissues. A syt11b HCR in situ hybridization and staining for HuC/D, EdU and Hoechst. SC cross sections from wild-type zebrafish subjected to SC transections and daily EdU injections are shown. Control SCs received daily EdU injections for 7 days before collection. Dotted lines delineate central canal edges. B, C syt11b integrated density was calculated in whole SC sections (ANOVA p: 0.0007) and in HuC/D+ neurons (ANOVA p: 0.0009). D vamp4 HCR in situ hybridization and staining for HuC/D, EdU and Hoechst were performed as described in (A). E, F vamp4 integrated density was calculated in whole SC sections (ANOVA p-value: 0.0012) and in HuC/D+ neurons (ANOVA p: 0.0043). G Numbers of iNeurons and proportions of EdU+ iNeurons at 0, 1 and 3 wpi. Statistical analysis was performed for iNeuron numbers (gray) (ANOVA p: 0.0040) and the proportions of EdU+ iNeurons (maroon) (ANOVA p: 0.1403). H Numbers of newborn neurons and proportions of syt11b+ iNeurons at 0, 1 and 3 wpi. Statistical analysis was performed for newborn neurons numbers (gray) (ANOVA p < 0.0001) and the proportions of syt11b+ iNeurons (maroon) (ANOVA p: 0.1403). I syt11b HCR in situ hybridization in Tg(vglut2a:RFP; gad1b:GFP) zebrafish at 0, 1 and 3 wpi. Violet and yellow arrowheads indicate vglut2a+ syt11b+ neurons and gad1b+ syt11b+ neurons, respectively. J Quantification of syt11b transcripts in vglut2a+ and gad1b+ neurons (ANOVA p < 0.0001). For all quantifications, SC cross sections 450 µm rostral to the lesion were analyzed. Data points in all bar charts indicate individual animals and sample sizes are indicated in parentheses. Brown-Forsythe and Welch ANOVA were performed in (B, C, E, F, G and H) with Dunnett?s T3 multiple comparisons performed across different time points with 95% CI. Two-way ANOVA was performed in (J) with Tukey?s multiple comparisons with 95 % CI. All statistical tests are two-sided. In all bar charts, data are presented as mean values +/- SEM. ns; p > 0.05, *p ? 0.05, **p ? 0.01, ***p ? 0.001. Scale bars, 50 µm. Source data are provided as a Source Data file.

Tracing and CRISPR/Cas9 mutagenesis of iNeuron marker genes in zebrafish.A?C Rostral-to-caudal neuronal tracing at 1 wpi. Biocytin was applied 4 mm rostral to the lesion and co-labeled with syt11b and HuC/D. SC cross sections 750 ?m caudal to the lesion are shown. Dotted lines delineate central canal edges. Biocytin-labeled neurons were normalized to total neurons in (B). Biocytin-labeled iNeurons were normalized to traced neurons in (C). D Pipeline for CRISPR/Cas9 mutagenesis of iNeuron marker genes and regeneration assessment. E Swim endurance in CRISPR/Cas9 targeted animals at 4 wpi. Uninjected control siblings were used (ANOVA p < 0.0001). F Swim behavior in targeted animals at 4 wpi (ANOVA p < 0.0001). Swim distance was tracked under water current velocities of 0, 10 and 20 cm/s. G Glial bridging in targeted animals at 4 wpi. Gfap+ bridges are shown at the lesion site of atf3, vamp4 and syt11a/b crispants relative to controls. Percent bridging represents the cross-sectional area of the glial bridge at the lesion (0 µm) relative to the intact SC 750 µm rostral to the lesion (ANOVA p: 0.0049). H Anterograde axon tracing in targeted animals at 4 wpi. Biocytin was applied rostrally and analyzed 600 ?m (proximal) (ANOVA p < 0.0001) and 1500 ?m (distal) (ANOVA p < 0.0001) caudal to the lesion. Representative biocytin traces are shown at the proximal level. Axon growth was first normalized to biocytin 450 and 750 µm rostral to the lesion, and then to labeling in wild-type controls. Data points in all bar charts indicate individual animal and sample sizes are indicated in parentheses. Brown-Forsythe and Welch ANOVA were performed in (G) with Dunnett?s T3 multiple comparisons performed across different time points with 95% CI. Ordinary one-way ANOVA was performed in (E) with Dunnett?s multiple comparison with 95% CI. Two-way ANOVA with Tukey?s multiple comparisons (95% CI) was performed in (F, H). In all bar charts, data are presented as mean values +/- SEM. All statistical tests are two-sided. *p ? 0.05, **p ? 0.01, ***p ? 0.001. Scale bars, 50 µm. Source data are provided as a Source Data file.