Ontogeny of the spinal cord dorsal horn

Abstract

The dorsal horn of the mammalian spinal cord is an exquisite example of form serving function. It is comprised of diverse neuronal populations stacked into laminae, each of which receives different circuit connections and plays specialized roles in behavior. An outstanding question is how this organization emerges during development from an apparently homogeneous pool of neural progenitors. Here, we found that dorsal neurons are diversified by time, with families of related cell types born as temporal cohorts, and by a spatial-molecular gradient that specifies the full array of individual cell types. Excitatory dorsal neurons then settle in a chronotopic arrangement that transforms their progressive birthdates into anatomical order. This establishes the dorsal horn laminae, as these neurons are also required for spatial organization of inhibitory neurons and sensory axons. This work reveals essential ontogenetic principles that shape dorsal progenitors into the diverse cell types and architecture that subserve sensorimotor behavior.

Keywords: Spinal cord; Zic; cell types; development; dorsal horn; interneurons; lamina; ontogeny; structure; temporal specification.

Figure 1:. A single-cell transcriptomic atlas of the mouse spinal cord post-neurogenesis.

(A) Neuronal birthdate labeling in spinal cord sections from E16.5 embryos that received an EdU pulse (yellow) to mark dividing cells at the indicated timepoints between E10.5 and E13.5, with DAPI staining (blue). (B) Neuron counts per section of inhibitory (Pax2+, red) or excitatory (Lmx1b+, green) neurons that were EdU positive following EdU pulses at the indicated timepoints. (C) UMAP plot of single-cell RNA-sequencing data taken from E14.5 and E16.5 wild-type mouse embryos showing all neurons, and annotated by Group (1 = Early-born and ventral, 2 = dILA, 3 = dILB) and by neuronal classes – including six dILA and six dILB families, each shown in a distinct color. (D) Dot plot showing gene expression in each neuron class. Dot size indicates the percent of each cluster expressing each gene and dot shade indicates the average expression level of each gene in cells of each cluster. (E) Diagram illustrating the relationships between embryonic neuron classes and postnatal/adult neuron classes from Russ et al., (2021). Line thickness is inversely proportional to centroid-centroid distance in principal component (PC) space, with the strongest relationships emphasized with black (versus gray) lines. (F) A dendrogram of centroid-centroid distance in PC space for refined neuron types among all dIL neurons. Each terminal branch number corresponds to named cell-types in Figure S3. Selected branch points are annotated with top DEGs. Where branches encompass large parts of dIL families they are shaded with colors representing each dIL family (see legend). Numbers: (A,B) EdU birthdating data was derived from 79 E16.5 mouse embryos given EdU on E10.5 (n=13 from four litters), E11.0 (n=10 from three litters), E11.5 (n=15 from four litters), E12.0 (n=7 from three litters), E12.5 (n=12 from four litters), E13.0 (n=10 from three litters) or E13.5 (n=12 from four litters). (C-F) Single-cell RNA-sequencing data was derived from fifteen litters of mouse embryos (n=3 E14.5, n=12 E16.5), and of which data from N=94048 neurons were extracted. Scale bars: (A) 250 μm Abbreviations: DEG – differentially-expressed gene, Exc – excitatory, Inh – inhibitory, UMAP – uniform manifold approximation and projection.

Figure 2:. dIL families represent temporal cohorts.

(A) UMAP feature plots of gene expression of each indicated gene in single-cell RNA-sequencing data, including temporal transcription factors (Onecut1, Zfhx3, Nfib), a marker of migrating and commissural neurons (Robo3) and a marker of neural progenitors (Sox2). (B) PHATE plot of dIL neurons, together with their early-born dI4/5 counterparts and immature neurons. (C) Dot plot of temporal transcription factors among dIL families, with dot size representing the percent of cells in each type expressing the gene and shade representing average expression level. (D) Composition over time of dILA families (top) and dILB families (bottom), together with immature dILA/B neurons (black) in single-nucleus RNA-sequencing data from precisely-timed mouse embryo dissections. (E-H) EdU birthdating of dILA and dILB families. (E,G) Examples of dILA (E) and dILB (G) family neurons (indicated by marker RNA or protein in green and red) co-labeled with EdU (blue) injected at the indicated timepoints. White arrows indicate examples of dIL family subtypes co-labeled with EdU. Dashed white lines indicate approximate anatomical boundaries of the spinal cord, including the pial outer boundary, the gray/white boundary, and the superficial dorsal horn. Boxes in dashed-white lines indicate the field of view magnified in the inset. (F,H) Cumulative percent of each dILA (F) or dILB (H) family co-labeling with EdU pulsed at the indicated embryonic time point (x-axis) showing mean +/− standard deviation for all litters at each timepoint, connected by sigmoidal fit lines for each family. The imputed 50% cumulative birthdating timepoint for each family is indicated by a colored vertical line beneath the x-axis. Numbers: (A-C) Single-cell RNA-sequencing data is the same as Figure 1C (see legend). (D) Enriched spinal cord neuron single-nucleus RNA-sequencing data (N=349,306 neurons) was retrieved from Qiu et al., (2024) and refined cell type labels were transferred from our data (Figure 1C–F); from this, 74,764 dILA and 71,834 dILB neurons were analyzed across 24 time points. (E-H) EdU birthdating data was derived from 34 E16.5 mouse litters from a separate cohort to Figure 1A,B, and data from individual embryos per litter were averaged. Pregnant mothers were given EdU on E10.5 (n=3–5 litters), E11.0 (n=4 litters), E11.5 (n=5 litters), E12.0 (n=5 litters), E12.5 (n=4–5 litters), E13.0 (n=5 litters) or E13.5 (n=2–5 litters). Scale bars: (E) 250 μm, 25 μm for insets. Abbreviations: TTF – temporal transcription factor

Figure 3:. dILB neurons settle into a chronotopic map.

(A) Images of embryonic spinal cord sections taken from separate E16.5 embryos that received EdU (blue) pulses at the indicated timepoints between E10.5 and E13.5, with immuno-labeling for the excitatory marker Lmx1b (green) and the inhibitory marker Pax2 (red). Hemisections (top row) and magnifications (bottom row) are shown, with the approximate region of Lmx1b+/EdU+ nuclei indicated with brackets (yellow). Solid and dashed white lines indicate approximate anatomical boundaries of the spinal cord. (B) Density plots of coordinates for birth-dated dILA neurons (Pax2+/EdU+) or dILB neurons (Lmx1b+/EdU+), for EdU pulse experiments at each of the indicated timepoints. Separate perimeters were generated for each replicate per timepoint and for five proportions of labeled cells (50%, 25%, 15%, 5%, 2.5%), and were colored in order of decreasing opacity. Replicates consist of all embryos of a single litter labeled with EdU at the indicated timepoint. (C) Scatter plots of neurons in data from Xenium spatial transcriptomics, annotated based on predicted dIL family identity. Adjacent histograms demonstrate the frequency of neurons belonging to dILA or dILB families by depth. Probe sets used were the mouse brain panel v1.1 (248 genes) plus 100 custom genes based on single-cell RNA-sequencing generated here (Table S9). Numbers: (A,B) n=79 embryos, 7–15 per timepoint. Note that data is reprised from Figure 1 (A,B), see legend for details. (C) For spatial transcriptomics, cell coordinates were derived from six whole sections and seven hemisections (N=9.5 sections) among n=3 replicate wild-type E16.5 spinal cords. From these sections 54,166 cells were identified as neurons and plotted, of which 27,714 were annotated as dIL neurons. Scale bars: (A) 250 μm (top row), 50 μm (bottom row).

Figure 4:. dILB neurons are required for laminar organization of the dorsal horn.

(A,B) Images of E16.5 spinal cord dorsal horn hemi-sections from control and mutant embryos. (A) Sections from Ptf1a^+/+ and Ptf1a^Cre/Cre embryos showing loss of Pax2-immunolabeled neurons (red), tissue morphology (DAPI, gray), an example marker of an excitatory interneuron (Plscr5 RNA, gray), and proprioceptive sensory afferents (PV antibody, gray). (B) Sections from Gsx1^+/−; Gsx2^+/+ and Gsx1^−/−; Gsx2^−/− embryos showing loss of Lmx1b-immunolabeled neurons (green), tissue morphology (DAPI, gray), an example marker of an inhibitory interneuron (Osr1 RNA, gray), and proprioceptive sensory fibers (PV antibody, gray). (C,D) Spatial distribution of individual gene or protein markers or marker sets (as indicated at the top of each box), showing the normalized proportion of total signal per dorso-ventral bin (color intensity) where the pia is set as the dorsal boundary and approximate laminar positions are shown, for either (C) absence of dILA or (D) absence of dILB neurons. Biological replicates are grouped together as individual columns for either Ctl (control – Ptf1a^+/+ or Gsx1^+/−; Gsx2^+/+) or Mut (Ptf1a^Cre/Cre or Gsx1^−/−; Gsx2^−/−) embryos. (E,F) Scatter plots of neurons in data from Xenium spatial transcriptomics, annotated based on predicted dIL family identity. (E) Brachial control Ptf1a^+/+ embryo sections (n=3, left side) were compared to a single Ptf1a^Cre/Cre embryo section (n=1, right side). The control sections are also displayed in Figure S10. (F) Thoracic control Gsx1^+/+; Gsx2^+/+ embryo sections (n=3, left side) were compared to Gsx1^−/−; Gsx2^−/− embryo sections (n=3, right side). Control neurons were downsampled to match the number of mutant cells displayed. Cells not belonging to the indicated dIL cell type are plotted in gray. Solid and dashed white and yellow lines indicate approximate anatomical boundaries of the spinal cord. Numbers: (A,D) Representative images were derived from n=4–6 Ptf1a^+/+, n=3–7 Ptf1a^Cre/Cre, n=3–4 Gsx1^+/−; Gsx2^+/+ and n=3–4 Gsx1^−/−; Gsx2^−/− and histograms were derived from these same embryos (different numbers within each range were analyzed for each marker set depending on staining quality). (E,F) Spatial transcriptomics from E16.5 wild-type embryo brachial segments was previously described (Fig. S9). For E16.5 wild-type thoracic segments, cell coordinates were derived from five whole sections and one hemisection (N=5.5 sections), among n=3 replicate spinal cords. From these sections 27,758 cells were identified as neurons and plotted, of which 13,724 were annotated as dIL neurons. For E16.5 Ptf1a-null brachial segments, cell coordinates were derived from two whole sections (N=2 sections) among n=1 spinal cord. From these sections, 8,267 cells were identified as neurons, and of which 2,767 were dIL neurons. E16.5 brachial data was downsampled to 8,106 neurons of which 4,151 were dIL neurons to match. For E16.5 Gsx1/2-null thoracic segments, cell coordinates were derived from four whole sections and three hemisections (N=5.5 sections) among n=3 replicate spinal cords. From these sections, 24,737 cells were identified as neurons, and of which 12,217 were dIL neurons. E16.5 thoracic data was downsampled to 24,583 neurons with, also, 12,217 dIL neurons, to match. Scale bars: 250 μm Abbreviations: ctl – control, dc – dorsal columns, mut – mutant, PV - parvalbumin.

Figure 5:. A progenitor gradient of Zic transcription factors patterns dIL neurons.

(A) UMAP feature plots of Zic1-5 expression in E14.5/E16.5 single cell data. (B) Dot plots of Zic1 (green), Zic4 (red), and Zic2 (blue) RNA expression in refined cell types within each dIL family, with dot size representing the percent of cells in each type expressing the gene and shade representing average expression level. (C) Gradients of Zic gene expression align with UMAP structure, within the dILA2 refined types M-X. (D) Dot plot of dILA2 neuron subtypes, showing unique markers for each type as well as groups of gene expression which co-vary along the Zic axis. (E) dIL neurons (left) and progenitors (right) on the y-axis, ranked by Zic1 expression (x-axis), showing expression of each Zic gene. (F) Images of E11.5 spinal cord sections with in situ hybridization for Zic genes (left, merged; right single genes). (G) Image of E11.5 spinal cord section with immuno-labeling for Zic1 (green), a marker of newborn immature and migrating dIL neurons (Lbx1, red), and with Lhx2 (dI1, blue) shown for orientation. The ventricular zone is shown as the most medial area bounded by dashed white lines. Right panel shows higher magnification inset of newborn dIL neurons as they exit the ventricular zone. (H) Genetic strategy for perturbing the Zic gradient through double knockout (DKO) alleles of Zic1 and Zic4. (I) E12.5 spinal cords from control and Zic1/4^DKO embryos, showing in-situ hybridization for Satb2 (magenta) and DAPI (cyan). (J) Scatter plots of normalized coordinates of neurons identified in Xenium spatial transcriptomics data, annotated to show dILA2 refined subtypes in control brachial sections (n=3, left) and in Zic1/4^DKO sections (n=3, right). Spatial transcriptomics from E16.5 wild-type embryo brachial segments was previously described (Fig. S9). For E16.5 Zic1/4^DKO brachial segments, cell coordinates were derived from three whole sections and one hemisection (N=3.5 sections) among n=3 replicate spinal cords. From these sections, 17,505 cells were identified as neurons, and of which 7,376 were dIL neurons. E16.5 brachial data was downsampled to 17,347 neurons of which 8,852 were dIL neurons to match. (K) Quantification of the data from (J) showing the relative composition of neuron subtypes as a percentage of all neurons, including replicates and lines connecting control and mutant means. Numbers: (A-D) N=94,048 neurons from single-cell RNA-sequencing data generated here (see Figure 1C), (E) Single-cell RNA-sequencing data from dIL progenitors (N=1,963 cells) and dIL neurons (N=61,444 cells, see Figure 1C), (F-G) n=3 E11.5 wild-type embryos, (I) n=3 E12.5 wild-type and n=3 E12.5 Zic1/4^DKO embryos, (J-K) n=3 wild-type and n=3 Zic1/4^DKO E16.5 mouse embryos. Statistics: (K) Two-way ANOVA (neuron subtype × genotype). For each neuron subtype in dILA2, Tukey’s multiple comparisons test was performed for each pair of genotypes, the result of which is indicated by asterisks above the plot. *<0.05, **p<0.01, ***p<0.001 Scale bars: 100 μm, except 50 μm in the magnified panel of (G) at the right. Abbreviations: DKO – double knock-out, scRNA – single-cell RNA

Figure 6:. Four principles governing dorsal horn ontogeny.

(A) Mechanisms of neuronal diversification: (1) Temporal cohorts of neurons are progressively born from the pdL progenitor zone, producing six cohorts of dILA and dILB neurons each; (2) dorsoventral position within the pdL progenitor zone produces variant neuron types particular to each dIL family; (3) previously known Notch-mediated mechanisms allow dILA and dILB neurons to be simultaneously born from asymmetric terminal progenitor divisions. (B) Mechanisms of laminar formation: (4) Sequentially born excitatory dILB families settle into spatially adjacent laminae.