Conserved genetic signatures parcellate cardinal spinal neuron classes into local and projection subsets (original) (raw)

• Journal List
• HHS Author Manuscripts
• PMC8612134

Science. Author manuscript; available in PMC 2021 Nov 24.

Published in final edited form as:

PMCID: PMC8612134

NIHMSID: NIHMS1748840

Peter J. Osseward, II,1,2,† Neal D. Amin,1,‡ Jeffrey D. Moore,3 Benjamin A. Temple,1,2 Bianca K. Barriga,1,4 Lukas C. Bachmann,1 Fernando Beltran, Jr.,1 Miriam Gullo,1 Robert C. Clark,1 Shawn P. Driscoll,1 Samuel L. Pfaff,1,* and Marito Hayashi1,†§

Peter J. Osseward, II

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

2Neurosciences Graduate Program, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92037, USA.

Neal D. Amin

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

Jeffrey D. Moore

3Howard Hughes Medical Institute, Department of Molecular and Cellular Biology, Center for Brain Science, Harvard University, Cambridge, MA, USA.

Benjamin A. Temple

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

2Neurosciences Graduate Program, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92037, USA.

Bianca K. Barriga

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

4Biological Sciences Graduate Program, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92037, USA.

Lukas C. Bachmann

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

Fernando Beltran, Jr.

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

Miriam Gullo

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

Robert C. Clark

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

Shawn P. Driscoll

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

Samuel L. Pfaff

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

Marito Hayashi

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

1Gene Expression Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92037, USA.

2Neurosciences Graduate Program, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92037, USA.

3Howard Hughes Medical Institute, Department of Molecular and Cellular Biology, Center for Brain Science, Harvard University, Cambridge, MA, USA.

4Biological Sciences Graduate Program, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92037, USA.

‡Present address: Department of Psychiatry and Behavioral Sciences, Stanford University, Stanford, CA 94305, USA.

§Present address: Howard Hughes Medical Institute, Department of Cell Biology, Harvard Medical School, Boston, MA 02115, USA.

†These authors contributed equally to this work.

Author contributions:

Conceptualization: M.H., P.J.O., and S.L.P. Methodology: M.H., P.J.O., J.D.M., and S.P.D. Software: S.P.D., J.D.M., and P.J.O. Validation: M.H., P.J.O., L.C.B., and R.C.C. Formal analysis: S.P.D., J.D.M., P.J.O., and M.H. Investigation: M.H., P.J.O., N.D.A., B.A.T., B.K.B., L.C.B., F.B., and R.C.C. Resources: M.G., S.P.D., J.D.M., P.J.O., and M.H. Data curation: S.P.D., J.D.M., and P.J.O. Writing, original draft: M.H., P.J.O., and S.L.P. Writing, review and editing: M.H., P.J.O., N.D.A., J.D.M., B.A.T., B.K.B., S.P.D., and S.L.P. Visualization: M.H., P.J.O., and S.P.D. Supervision: M.H., P.J.O., and S.L.P. Project administration: M.H., P.J.O., and S.L.P.

Supplementary Materials

Figure S1.

GUID: E035A152-BF6A-4163-887F-FA59040F4C7D

Figure S2.

GUID: 8531A6FA-C1C7-48D2-8464-0E6702E39F09

Figure S4.

GUID: DF2EE545-2953-4DEE-9E75-E7651D49CA78

Figure S3.

GUID: 8E470300-AE36-4428-B661-A5FB4050AD32

Figure S5.

GUID: B803F7D1-CC8C-4B26-B224-7DBC62BAB8C6

Figure S7.

GUID: 07E66302-4547-48CB-879E-5F46427C86F6

Figure S6.

GUID: 9B6753CE-02D3-4BF7-9ED9-ECB55849F735

Figure S8.

GUID: 332BC036-1226-4B37-B040-855C0F725A2E

Figure S9.

GUID: D3354DF4-01F8-4C32-A802-98DFDC9C1F51

Figure S10.

GUID: B9C18CB2-5B33-40A2-9092-DA0DFC45F5E5

Figure S11.

GUID: 35507D06-836D-46CE-A82D-299CE66CEC88

Figure S13.

GUID: 9AEDA017-A86C-4000-8ED9-FFAAB3107E85

Figure S14.

GUID: 5471B712-A42F-4BAE-8E1F-AC1734A0CF34

Figure S12.

GUID: 70DA6897-FDE6-44BD-AF1D-92E09AD7A82C

Figure S15.

GUID: B7F3ADE8-1EE1-4DDD-8BC7-FB6624B0F916

Figure S16.

GUID: 86DF7961-0326-4583-81B6-CA0D5564BCBF

Figure S17.

GUID: 1A78D076-9636-4380-85FF-732C871044AF

Figure S18.

GUID: 60E99676-0B95-43BF-B3CC-2F27D9A825E0

Figure S19.

GUID: 26435D8A-A07E-4F66-A873-FA5C5C6EC69B

Figure S21.

GUID: 700C66ED-3D02-407F-879B-C81B233D8C77

Figure S20.

GUID: EAEEB3B6-90A4-4C6B-8ECF-6D403DD3F781

Table S2.

GUID: B5D07433-F689-4E80-8410-1E858A52D5FB

Table S3.

GUID: F68CE93B-31D8-4580-B4DA-F37B45547BB2

Table S1.

GUID: 53F1CCC2-42F1-4A08-A009-49125F81C58E

Table S4.

GUID: 8C02FDE7-1D79-4D5B-835B-127AD2568BE0

Table S5.

GUID: F36BBFDC-80F1-4C22-A50A-D56112CDA01B

Supplementary Material.

GUID: E3B585D6-4747-468D-9CF5-DBAFE3AE5E4D

Data Availability Statement

All data are available in the main text or the supplementary materials. The accession number for the scRNA-seq data reported in this paper is NIH GEO: GSE171607, GSE171605, and GSE171606.

Abstract

Motor and sensory functions of the spinal cord are mediated by populations of cardinal neurons arising from separate progenitor lineages. However, each cardinal class is composed of multiple neuronal types with distinct molecular, anatomical, and physiological features, and there is not a unifying logic that systematically accounts for this diversity. We reasoned that the expansion of new neuronal types occurred in a stepwise manner analogous to animal speciation, and we explored this by defining transcriptomic relationships using a top-down approach. We uncovered orderly genetic tiers that sequentially divide groups of neurons by their motor-sensory, local-long range, and excitatory-inhibitory features. The genetic signatures defining neuronal projections were tied to neuronal birth date and conserved across cardinal classes. Thus, the intersection of cardinal class with projection markers provides a unifying taxonomic solution for systematically identifying distinct functional subsets.

From a functional perspective, spinal neurons can be divided along several axes, including motor-sensory, excitatory-inhibitory, and locally connected neurons for intracord processing versus projection neurons for communication with the brain (1–6). Neuron heterogeneity is often characterized on the basis of landmarks such as their neurotransmitter type, connectivity, cytoarchitecture, morphology, physiology, developmental origin, and genetic profile (7). In the spinal cord, separate cardinal neuron classes arise from molecularly distinct progenitor cells arrayed along the dorsoventral axis of the neural tube (2, 8). Neurons within each cardinal class share properties such as the same neurotransmitter identity and have been targeted in functional studies (1, 2). However, it is increasingly apparent that each cardinal class is itself composed of heterogeneous populations of neuron types, impeding a detailed understanding how spinal circuits function (9–13). We sought to determine whether there was a coherent logic for spinal neuron diversification in mice that could be used to systematically describe the heterogeneity within all the cardinal classes.

One perspective of neuronal diversity is that it arose through evolution in order to expand neural functions (14–16). In this view, primitive neuron types served as the precursors for more specialized subtypes, leading to the prediction that individual neuronal attributes may have emerged in a stepwise fashion with ordered hierarchical relationships (14, 16, 17). This analogy to animal speciation prompted us to investigate whether molecular and cellular correlates for this type of stepwise diversification could be detected among large heterogeneous populations of spinal neurons. Thus, rather than focusing on fine-grained molecular differences that may or may not correspond to functional cellular features, we used a top-down approach to identify the transcriptomic signatures linked to the emergence of different neuronal types. At each branch point in our established hierarchy, we defined the molecular, developmental, cytoarchitectural, neurotransmitter, and connectivity properties of the neurons. Our analysis uncovered molecular markers for local and projection neurons regardless of cardinal class or neurotransmitter. By combining markers for cardinal classes with those for conserved genetic signatures for projection status, we established a simple combinatorial matrix that systematically identified discrete subsets of spinal neurons.

Group-E and -H neurons divide the spinal motor-sensory circuitry

To enrich for neurons while minimizing technical and biological variability, we genetically labeled glutamatergic and γ-aminobutyric acid (GABA)–ergic neurons in postnatal day 0 (P0) mice, sorted lumbar spinal neurons separately, pooled them, and performed single-cell RNA-sequencing (scRNA-seq) on the combined sample (Fig. 1A and fig. S1). At P0, spinal neurons have been generated and the basic functional features of adult spinal circuitry have been established, but expression persists for many developmental markers of neuronal subtype diversity (2, 9, 10, 18). Using marker analyses, 6743 cells passed quality control and were identified as neurons (Fig. 1A). Glutamatergic and GABAergic neurons largely segregated in the uniform manifold approximation and projection (UMAP) space, and a graph-based comparison using the highly variable genes conservatively identified 45 distinct clusters (Fig. 1, ​B and ​C). As expected, many of these clusters were enriched for genes that have been previously identified as marker genes for cardinal classes (Fig. 1D and table S1) (1, 2).

scRNA-seq identifies relationships among groups of spinal neurons.

(A) Lumbar segments from P0 Vglut2:Cre; Ai14 and Vgat:Cre; Ai14 neonates (two females and one male from each litter) were microdissected. tdTomato+ cells were sorted and pooled as one sample before scRNA-seq. Using marker analyses and quality control measures, 6743 cells were identified as neurons (94.4% of total cells detected). These neurons had 1012 median genes per cell and 1606 median unique molecular identifiers (UMIs) per cell. (B) UMAP visualization of neurons identifies that glutamatergic neurons and GABAergic neurons largely segregate in UMAP space. (C) Forty-five identified clusters are color coded in UMAP space. (D) Colored dots with numbers correspond to clusters from Fig. 1C. As expected, clusters corresponded to known cardinal neuron classes, defined by their marker expression. V1 and V2b markers were found in the same clusters (12). Several clusters are not shown because they were not enriched for obvious cardinal class marker genes (table S1). (E) Schematic outlining methods to identify neuronal subtypes. (Top) Dimensionality reduction tools can be used to define the maximum number of genetically distinct clusters. Functions must be independently defined for each cluster. (Bottom) Hierarchical comparisons of neuronal populations can be used to define branch points. The functional differences at the branch points can be used to identify increasingly more specific neuronal attributes linked to function. (F) Schematic of analysis used to investigate the higher-order relationships across spinal neuron populations conducted with iterative _k_-means 2 divisive clustering. Data for each color-coded division are indicated for (G) to (K). (G to K) (Left) UMAP color coded according to identified division. (Right) Volcano plot depicting differentially expressed genes across each _k_-mean division. [(G) and (I)] Further analysis of highlighted genes are available in Fig. 2 and fig. S2. [(H), (J), and (K)] Glutamatergic- and GABAergic-related genes are highlighted. (G) Group-E: 3360 cells, Jaccard similarity [across bootstraps (supplementary materials, materials and methods)] = 0.960; group-H: 3383 cells, Jaccard similarity = 0.959. (H) Group-E (Glut): 1851 cells, Jaccard similarity = 0.960; group-E (GABA): 1509 cells, Jaccard similarity = 0.922. (I) Group-Z: 1707 cells, Jaccard similarity = 0.867; group-N: 1676 cells, Jaccard similarity = 0.809. (J) Group-N (Glut): 1040 cells, Jaccard similarity = 0.791; group-N (GABA): 636 cells, Jaccard similarity = 0.613. (K) Group-Z (Glut): 768 cells, Jaccard similarity = 0.769; group-Z (GABA): 939 cells, Jaccard similarity = 0.760.

To explore the conserved relationships that exist among lumbar spinal neurons, we developed a divisive clustering pipeline using the _k_-means algorithm to split spinal neurons into pairwise groups iteratively on the basis of their overall transcriptional composition (Fig. 1, ​E and ​F). We focused on transcription factors because they regulate neuronal identity and function (2, 14, 15). We anticipated that this approach would first reveal the glutamatergic and GABAergic neurotransmitter division in the cell population because it is such a prominent phenotypic difference, and the two cell types were isolated from separate litters. However, the neuronal populations split along an axis unrelated to neurotransmitter identity (Fig. 1, ​B and ​G) (11, 19). Differential gene expression of the first division indicated that transcription factor Ebf1 is highly expressed in a cell group we defined as group-E neurons, whereas transcription factor Hoxc10 is enriched in a different group we termed group-H neurons (Fig. 1G). Histological analysis of Ebf1+ and Hoxc10+ neurons identified a spatial division of these groups (Fig. 2, A and ​B, and fig. S2). Ebf1+ neurons were located in laminae I to III of the superficial dorsal horn, corresponding to the site where many exteroceptive sensory inputs terminate. By contrast, Hoxc10+ neurons were predominantly located in laminae IV to X, areas involved in proprioception and motor control (3, 4). Although group-E and group-H neurons broadly differ in motor-sensory functions associated with their positions, both types are composed of a mixture of glutamatergic and GABAergic neurons (Fig. 1, B and ​G, and fig. S3) (3, 4).

Branch points in the hierarchical comparisons define spinal neuron groups with discrete locations.

(A to D) (Top) Immunostaining against representative markers for group-E, -H, -N, and -Z neurons in lumbar segments at P0. All images are 20-μm cryosections. Scale bars, 100 μm. (Bottom) Spatial analysis of representative markers displayed on right half of lumbar spinal cord. Contours indicate density of representative markers at the 10th to 90th percentiles. y axis = 0 represents the top end of the central canal, and x axis = 0 represents the midline. (E) Matrix of pairwise mean colocalization rates between group-Z, -N, and -E markers. “N.M.” (blue boxes with X) indicates not measured because of antibody incompatibility (figs. S4 and S5). NeuN+ Nfib+ neurons were quantified to exclude Nfib+ glial cells. (F) Matrix of pairwise weighted Jaccard distance between group-Z, -N, and -E markers. 1.0, no overlap; 0, complete overlap. Hierarchical clustering dendrogram reveals closer spatial associations among marker genes within each group than between separate groups. (G) Immunostaining of Zfhx3 and NeuroD2 in lumbar spinal cord at P0 reveals low incidences of colocalization. Cryosection, 20 μm. Scale bar, 100 μm. This image is a composite of Fig. 2, C and ​D. (H) (Left) Transcriptomic hierarchy of lumbar spinal neurons. (Right) Schematic of spatial enrichments of group-Z, -N, and -E neurons.

Group-N and -Z neurons are arrayed along the mediolateral axis of the spinal cord

Using the hierarchical-divisive strategy outlined above (Fig. 1F), we further split the group-E neurons into two subgroups and found that this division identifies glutamatergic-related and GABAergic-related genes and their corresponding neurons (Fig. 1, ​B and ​H) (2, 4). Division of the group-H population produced a different result. Rather than splitting group-H neurons into glutamatergic and GABAergic cell populations, we found a subgroup of neurons enriched for Nfib, NeuroD2, and Prox1 (designated group-N) and a separate subgroup enriched for Zfhx3, Zfhx4, and FoxP2 (group-Z) (Fig. 1I and table S2). Histological analysis found that these two sets of genes spatially divided the spinal cord into large wedge patterns: group-N cells located dorsomedially and group-Z cells positioned ventrolaterally (Fig. 2, C and ​D, and fig. S4). Group-N and group-Z neurons occasionally intermingle (Fig. 2, C to ​H, and fig. S5), but few cells coexpressed markers for both group-N and group-Z neurons (Fig. 2, E and ​G, and fig. S5). Although group-N and group-Z neurons are defined by gene sets, we found that NeuroD2 and Zfhx3 were representative markers of group-N and -Z neurons, respectively (Figs. 1I and 2, C and ​D; and figs. S2 to S5). Our sequencing was conducted on cells from the lumbar spinal cord, but we found that NeuroD2 and Zfhx3 labeled neurons in the same mediolateral pattern along the entire rostrocaudal axis of the cord (figs. S6 and S7). Further _k_-means splitting of group-N and -Z neurons resolved the glutamatergic and GABAergic subsets in each population (Fig. 1, ​J and ​K, and fig. S3) (2, 4). Our findings unmask a stepwise hierarchy of genetic divisions among spinal neurons with layers that supersede neurotransmitter physiology (Fig. 2H).

Group-N and -Z neurons have distinct birth dates

Although the group-E–group-H division appeared to correspond to neurons involved in exteroception and motor control, respectively, we hypothesized that the group-N and -Z subsets derived from the group-H population might segregate along alternative features. To investigate a possible developmental distinction between group-N and group-Z neurons, we performed scRNA-seq of embryonic day 12.5 (E12.5) spinal cord cells composed of the embryonically defined cardinal cell classes (Fig. 3, A and ​B; fig. S8; and table S3). The E-H division was not evident at E12.5; however, the group-N and group-Z neuronal populations were detected (Fig. 3, C and ​D, and figs. S8 and S9). We birth-dated the neurons marked by NeuroD2 and Zfhx3 at P0 using EdU pulses embryonically. The peak of Zfhx3+ neuron birth was at E10.5, whereas the peak of NeuroD2+ neuron generation was later at E13.5 (Fig. 3E and fig. S10). To examine whether group-N and group-Z identity is dependent on neuronal connectivity, we genetically ablated proprioceptive sensory neurons or motor neurons. Despite the absence of sensory and motor neurons, Zfhx3 and NeuroD2 expression patterns were unaltered (fig. S11). Furthermore, we found that Zfhx3 and NeuroD2 expression was maintained in adult spinal neurons (P70), approximating the labeling pattern observed at P0 (Fig. 3F). Taken together, our findings indicate that the group-N and -Z populations have distinct timings of neurogenesis and are present embryonically through adulthood (Fig. 3G).

The group-N and -Z division is linked to neuronal birth date.

(A) Cells along the entire spinal cord were isolated from two E12.5 embryos from two litters and were subjected to scRNA-seq as separate samples. These datasets were merged post hoc. We bioinformatically identified 2945 neurons, which were selected for further analysis (fig. S8). These neurons had 3566 median genes per cell and 12,189 median UMIs per cell. Glutamatergic and GABAergic neurons are labeled on the UMAP. (B) Cardinal classes segregate in UMAP space. Twenty-seven clusters were identified (fig. S8). On the basis of the top 20 differentially expressed genes for each cluster, cardinal classes were identified. Clusters from the same cardinal class were merged (table S3). (C) Higher-order groups of spinal neurons were identified with _k_-means 2 divisive clustering (Fig. 1F). Embryonic neurons are color coded according to the first _k_-means 2 division. Population 1: 1971 cells, Jaccard similarity = 0.978; population 2: 974 cells, Jaccard similarity = 0.955. Further divisions are provided in fig. S9. (D) The first _k_-means 2 division of embryonic neurons largely corresponded to group-N and group-Z identity established from P0: 79% of population 1 [(C), blue cells] neurons expressed group-Z markers, and 84% of population 2 [(C), red cells] neurons expressed group-N markers. Thus, the _k_-means 2 division identifies group-N and group-Z neurons at an early embryonic stage. (E) EdU was injected into multiple pregnant females (E10.5 to E14.5), and lumbar segments were collected from their pups at P0. Representative markers Zfhx3 and NeuroD2 were used to birth date group-Z and group-N neurons, respectively. The peak birth date for group-Z neurons was E10.5, whereas the peak for group-N neurons was E13.5. Each dot represents the colocalization rate for a single animal. Welch one-way analysis of variance (ANOVA) comparing every time point to every other time point, with Dunnett T3 correction for multiple comparisons. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001. Cryosections, 20 μm. Scale bars, 100 μm. (F) Immunostaining of Zfhx3 and NeuroD2 in lumbar spinal cord at P70. The mediolateral location and marker expression patterns in adult P70 animals were similar to P0. Cryosection, 25 μm. Scale bar, 100 μm. (G) Summary. Group-Z neurons are born before group-N neurons. Molecular markers for both groups are expressed embryonically, postnatally, and into adulthood. Markers for group-Z and group-N neurons are not coexpressed.

Group-N and -Z neurons divide each cardinal neuron class into subsets

We set out to understand how group-E, group-N, and group-Z neurons identified here spatially and molecularly intersect with the previously characterized cardinal classes. Cre lines labeling the cardinal classes dI1, dI3, dI4/dILA, dI5/dILB, V1, V2a, and V3 were immunostained for Ebf1 (group-E), NeuroD2 (group-N), and Zfhx3 (group-Z) (Fig. 4A). Group-E neurons were found to be restricted to the dI4/dILA and dI5/dILB classes (figs. S12 and S13), whereas all tested cardinal classes were composed of a mixture of both group-N and group-Z neurons derived from the group-H population (Fig. 4, B to ​F, and figs. S12 and S13). Analysis of E12.5 scRNA-seq data confirmed that the dI2, dI6, and V2b classes also were a mixture of both group-N and group-Z neurons (fig. S13) (13). Moreover, the group-N neurons from each cardinal class were always medial to the group-Z cells (Fig. 4, C to ​E and ​H, and fig. S14) (9, 10, 20). Although each cardinal class contained both group-N and group-Z neurons, the ratio of these types varied within each class (Fig. 4, F and ​H, and figs. S13 and S15) (10, 20). Group-N and group-Z neurons do not obey cardinal class divisions; rather, all cardinal classes are a mix of group-N and group-Z cells with discrete mediolateral positions. Thus, the group-N–group-Z division is an orthogonal molecular and spatial axis that intersects with cardinal class identity (Fig. 4H).

Each cardinal class is a mixture of both group-N and -Z neurons.

(A) Summary of mouse lines used to label cardinal classes. Each cardinal class has a characteristic transcription factor expression, neurotransmitter, and axon laterality profile.

(B) Atoh1:Cre was crossed to a tdTomato reporter to indelibly label dI1 neurons. Immunostaining for NeuroD2 (group-N) and Zfhx3 (group-Z) reveals discrete dI1 populations marked by the two genes. Asterisk in the dorsal horn indicates ectopic sensory fibers from the DRG (29). Cryosection, 20 μm. Scale bar, 100 μm.

(C) (Left) Contours indicate density of dI1 neurons at the 10th to 90th percentiles. Mediolateral density displayed above a two-dimensional contour plot. (Right) Contours indicate densities of group-N and group-Z dI1 neurons at 50th to 90th percentiles. Group-N dI1 neurons were generally located medially, whereas group-Z dI1 neurons were located laterally.

(D) Contours and mediolateral densities of indicated cardinal classes. Contours indicate density of cardinal class neurons at the 10th to 90th percentiles. Dorsoventral densities are provided in fig. S14.

(E) Distributions of group-N and group-Z neurons within each cardinal class. Contours indicate density of interneuron subpopulations at the 50th to 90th percentiles. Each cardinal class was split into medial group-N and lateral group-Z neurons.

(F) Quantification of group-Z and group-N neurons within each cardinal class. Each dot indicates one animal. The ratio of group-N to -Z neurons varies across cardinal classes.

(G) Summary of the cardinal class types within the hierarchical comparisons of neurons.

(H) Stylized overview of cell spatial relationships defined by coexpression of cardinal and N-Z markers. Cardinal classes have stereotyped dorsoventral settling positions. Neurons composing each cardinal class are split into medial group-N and lateral group-Z cells.

Local versus long-range projection neuron types underlie the N-Z division

To further understand the cellular features of group-N and group-Z neurons, we explored the connectivity patterns of spinal neurons. We first tested whether the group-N and group-Z neurons represented premotor interneurons, using monosynaptic rabies virus tracing. We found that both the group-N and group-Z populations contained ipsi- and contralateral premotor neurons (fig. S16A) (21–23). Next, using immunostaining and mouse genetics to visualize presynaptic terminals, we found that both group-N and group-Z neurons receive a variety of inputs (for example, dorsal root ganglia sensory, descending serotonergic, dI1, dI4–6, V1, and V3) (fig. S16, B and C) (24). Gene Ontology (GO) term analysis revealed that group-Z neurons were enriched for genes involved in axon and synapse development, including high neurofilament (Nef-l, -m, and -h) expression (Fig. 5, ​A and ​B). Consequently, we speculated that group-Z neurons may contain one or more populations of long-range projection neurons (25). Subsets of neurons from multiple cardinal classes have been shown to share projection targets, further raising the possibility that group-Z neurons represent long-range projections across neuronal classes (7, 9, 26, 27). We labeled spinothalamic, spinocerebellar, spinomedullary, and lumbocervical projection neurons with retrograde tracer cholera toxin subunit β (CTB) in combination with group-E, -N, and -Z marker immunostaining (Fig. 5C). Although group-N and group-Z neurons are approximately equal in number (fig. S6), the long-range projection neurons in both the lumbar and cervical cord were preferentially enriched for Zfhx3+ neurons (Fig. 5, D to ​G, and figs. S17 to S20). Because the single marker Zfhx3 does not label all group-Z neurons, we also performed the supraspinal labeling and immunostaining with additional group-Z markers, including Zfhx4, Otp, and Foxp2. We found that multiple group-Z markers identify additional long-range projection neurons, providing further evidence that many group-Z cells project long axons (figs. S18 and S19). By contrast, markers for group-N (NeuroD2 and Prox1) and group-E (Ebf1) rarely (0 to 7.9%) labeled long-range neurons (Fig. 5, D to ​G, and figs. S17 to S20). As a positive control, injection of CTB directly into the cord labeled group-N neurons near the injection site (fig. S20). Taken together, our labeling studies indicate that group-Z neurons are composed of long-range projecting cells, whereas group-N and group-E neurons appear to be preferentially composed of local interneurons.

Projection and long-range neurons express group-Z markers.

(A) Neuronal morphology-related GO terms are enriched in group-Z compared with group-N neurons (table S4). (B) Heatmap of neurofilament genes expressed in spinal neurons showing enrichment in group-Z neurons (log2 scale). (C) Schematic of retrograde tracing analysis. CTB was injected into the forebrain (thalamus), cerebellum (vermis), medulla, and cervical spinal cord, respectively, to retrogradely label long-range projection neurons within the spinal cord. Injections were done in neonates, and tissue was collected 4 to 5 days later. (D to G) Immunostaining against group-Z (Zfhx3), -E (Ebf1), and -N (NeuroD2 and Prox1) markers combined with indicated CTB labeling. White arrowheads indicate group-Z long-range neurons. Each dot on the graph indicates the colocalization rate of one animal. Cryosections, 20 μm. Scale bars, 100 μm. Welch one-way ANOVA comparing Zfhx3 to every other marker, with Dunnett T3 correction for multiple comparisons. *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001. (D) Spinothalamic neurons within cervical segments (region of greatest abundance) expressed group-Z markers. Spinothalamic neurons within lumbar segments are shown in fig. S17. (E) Spinocerebellar neurons within lumbar segments expressed group-Z markers. Additional markers are provided in fig. S18. (F) Spinomedullary neurons within lumbar segments expressed group-Z markers. Additional markers are provided in fig. S19. (G) Cervicolumbar neurons expressed group-Z markers. Additional markers are provided in fig. S20.

Although most of our studies focused on the lumbar spinal cord, we found one counter-example to the group-N–group-Z correlation with axon projection in the thoracic region. Namely, a fraction of spinocerebellar projection Clarke’s column neurons was found to express the group-N marker NeuroD2 (fig. S18) (28). The reason for this discrepancy is unclear; however, these cells form a distinct medial nucleus with unusual molecular features, such as the expression of Vglut1 rather than Vglut2 (29, 30). Previous gene profiling of Clarke’s column neurons also revealed that they cluster separately from other thoracic spinal neurons (28). The non–Clarke’s column spinocerebellar neurons in the thoracic cord express group-Z genes that are typical for projection neurons (fig. S18) (28).

In principle, the correlation between group-Z markers and projection neuron labeling might have arisen as a byproduct of long-range neurons predominantly being located laterally within the spinal cord where group-Z markers are expressed for independent reasons. Therefore, we sought to define the axon projection pattern of group-Z neurons in atypical locations within the spinal cord. A sparse population of Zfhx3+ neurons were situated in the very superficial lamina of the spinal cord (Figs. 2D and ​3E). The superficial Zfhx3+ neurons were labeled after injection of CTB into the thalamus (Fig. 5D). These unusual Zfhx3+ projection neurons were further confirmed to represent spinothalamic cells by using Nk1r/Tacr1 labeling and were also born early like other group-Z neurons (Fig. 3E and fig. S17) (5, 31). Near the midline, the boundary between group-N and group-Z is not absolute, and there is sparse intermingling of NeuroD2+ and Zfhx3+ neurons. To further test the relationship between group-N and group-Z markers and their axonal projection patterns, we examined the marker profile of projection neurons in this region of the spinal cord. Retrograde fills of projection neurons preferentially labeled the Zfhx3+ neurons but not the NeuroD2+ neurons (Fig. 5, ​E to ​G). Thus, the molecular features that distinguish group-N from group-Z neurons reveal a core phenotypic difference: Group-N (and group-E) neurons are restricted to forming local connections within the spinal cord, whereas group-Z neurons can also form long-range connections with distant targets (Fig. 6).

A unified taxonomic index of spinal neurons.

(A) Summary of relationships of neuronal attributes. Neuronal subtypes can be viewed as a composite of distinct properties that fall into a range of categories (circles). Neuronal subtypes arise by joining different combinations of these attributes. Cardinal class and N-Z markers simplify the task of defining neuronal subtypes because both correspond to multiple neuronal features (magenta and green lines). (B) Summary of neuronal subsets defined by coexpression of cardinal class markers and N-Z labels. This matrix provides a simple labeling index to define neuronal subsets. The postsynaptic targets of group-Z neurons within each cardinal population are predicted from projection neuron–tracing studies (26, 29, 31, 37, 38). Local connections and collaterals are not shown (26, 27).

Discussion

In this study, we sought to identify a unifying logic for spinal neuron diversification. We used a top-down approach to identify the conserved genetic relationships among heterogeneous cell types and used the branch points of this hierarchical system to find marker genes that were linked to neuronal birth date, cytoarchitecture, neurotransmitter, and projection status. The divide between group-E neurons marked by Ebf1 and group-H cells marked by Hoxc10 corresponded to a split between sensory and motor circuits, respectively (Fig. 2H). The H-neuron population was further divided into group-N and -Z types, in which NeuroD2 labeled medial, late born, locally projecting interneurons and Zfhx3 marked lateral, early born, projection neurons (Figs. 2H, ​3G, and ​6, ​A and ​B).

The N-Z division giving rise to local versus projection neurons may have occurred early in evolution because the genetic signatures for these two groups are conserved across broad populations of spinal neurons. Consistent with this, group-N and -Z genes are conserved in many vertebrates, and functional studies of these factors have revealed roles in neural diversification (32–34). The correlation between birth date and N-Z identity is consistent with neuronal diversity being tied to neurogenesis (9, 22, 35, 36).

We anticipated that the most basic characteristic associated with neurons is their neurotransmitter phenotype (11, 15, 19) but instead found that this was superseded by the distinctions between motor-sensory and local-projection neuron types (Fig. 2H). Motor neurons in different species, although homologous cells, use different neurotransmitters, suggesting that there is flexibility and adaptability in the signals for neurotransmission (15). At a molecular level, the same neurotransmitter phenotype across cardinal classes can be specified by different transcriptional pathways (2). These findings suggest that spinal neuron neurotransmitter identity can arise through convergent evolutionary processes, causing the genetic signatures for neurotransmitter identity to appear lower in the hierarchical relationships because there are fewer cells with similar genetic characteristics.

Because the lineage-based markers of cardinal interneurons are orthogonal to the conserved markers for the N-Z neuronal sets, we found that the intersection of these two systems systematically labeled neurons with specific dorsoventral and mediolateral positions, birth dates, neurotransmitter phenotypes, ipsilateral-contralateral projections, and local–long-range connectivity (Fig. 6, A and ​B, and fig. S21). In isolation, neither cardinal class nor N-Z markers precisely define the postsynaptic targets of spinal neurons (2, 9, 10). However, the spinal neurons innervating the lateral reticular nucleus, thalamus, cerebellum, and long-range spinal targets (propriospinal) can be predicted on the basis of the intersection of their cardinal and N-Z genetic status (Fig. 6B and fig. S21). By identifying markers associated with diversification for the purpose of expanding cellular functions—in this case, evolution and development—it may be possible to create more unified and practical systematics for identifying functionally different cell types in many systems.

Supplementary Material

Figure S1

Figure S2

Figure S4

Figure S3

Figure S5

Figure S7

Figure S6

Figure S8

Figure S9

Figure S10

Figure S11

Figure S13

Figure S14

Figure S12

Figure S15

Figure S16

Figure S17

Figure S18

Figure S19

Figure S21

Figure S20

Table S2

Table S3

Table S1

Table S4

Table S5

Supplementary Material

ACKNOWLEDGMENTS

We thank K. Lettieri, B. O’Leary, H. Forman, A. J. Levine, M. J. Sternfeld, K. L. Hilde, members of the Pfaff laboratory, the FCCF core (C. O’Connor and C. Fitzpatrick; NIH-NCI CCSG: P30 014195), the NGS core (N. Hah; NIH-NCI CCSG: P30 014195, the Chapman Foundation and the Helmsley Charitable Trust), the UCSD IGM Genomic Center (K. Jepsen and other members), the Biophotonics core (NIH-NCI CCSG: P30 014195 and the Waitt Foundation), the GT3 core (J. Naughton, J. Marlett, and R. Armendariz; NIH-NCI CCSG: P30 014195, NINDS R24 Core Grant, funding from NEI), and the Janelia Viral Vector Core for technical support and advice; M. Goulding, G. Gatto, members of the Goulding laboratory, E. J. Kim and E. M. Callaway for providing reagents and advice; and S. Ackerman, K. Asahina, K. Baldwin, T. Komiyama, and Y. Zou for discussions and advice.

Funding:

P.J.O. was supported by the Christopher and Dana Reeve Foundation. N.D.A. is supported by the NIMH T32 Training Grant T32MH019938 J.D.M. is supported by NIH-NICHD K99 HD096512. B.K.B. is supported by the UCSD Cell and Molecular Genetics Training Program (T32 GM007240) and the Timken-Sturgis Foundation. M.H. was supported by the Japanese Ministry of Education, Culture, Sports, Science, and Technology Long-Term Student Support Program and the Timken-Sturgis Foundation. S.L.P. is a Benjamin H. Lewis chair in neuroscience. This research was supported by funding from the Howard Hughes Medical Institute, the Christopher and Dana Reeve Foundation, the Sol Goldman Charitable Trust, and National Institutes of Health (1 U19 NS112959-01).

Footnotes

Data and materials availability:

All data are available in the main text or the supplementary materials. The accession number for the scRNA-seq data reported in this paper is NIH GEO: GSE171607, GSE171605, and GSE171606.

REFERENCES AND NOTES

14. Arendt D et al., Nat. Rev. Genet 17, 744–757 (2016). [PubMed] [Google Scholar]

15. Hobert O, Kratsios P, Curr. Opin. Neurobiol 56, 97–105 (2019). [PubMed] [Google Scholar]

17. Zeng H, Sanes JR, Nat. Rev. Neurosci 18, 530–546 (2017). [PubMed] [Google Scholar]

18. Ladle DR, Pecho-Vrieseling E, Arber S, Neuron 56, 270–283 (2007). [PubMed] [Google Scholar]

21. Stepien AE, Tripodi M, Arber S, Neuron 68, 456–472 (2010). [PubMed] [Google Scholar]

22. Tripodi M, Stepien AE, Arber S, Nature 479, 61–66 (2011). [PubMed] [Google Scholar]

24. Alvarez FJ, Villalba RM, Zerda R, Schneider SP, J. Comp. Neurol 472, 257–280 (2004). [PubMed] [Google Scholar]

26. Pivetta C, Esposito MS, Sigrist M, Arber S, Cell 156, 537–548 (2014). [PubMed] [Google Scholar]

31. Nishida K, Ito S, Eur. J. Neurosci 46, 2608–2619 (2017). [PubMed] [Google Scholar]

35. Greig LC, Woodworth MB, Galazo MJ, Padmanabhan H, Macklis JD, Nat. Rev. Neurosci 14, 755–769 (2013). [PMC free article] [PubMed] [Google Scholar]

36. Doe CQ, Annu. Rev. Cell Dev. Biol 33, 219–240 (2017). [PubMed] [Google Scholar]

38. Ruder L, Takeoka A, Arber S, Neuron 92, 1063–1078 (2016). [PubMed] [Google Scholar]