Abstract
The mammalian brain contains numerous neurons distributed across forebrain, midbrain, and hindbrain that project axons to the lower spinal cord and work in concert to control movement and achieve homeostasis. Extensive work has mapped the anatomic location of supraspinal cell types and continues to establish specific physiological functions. The patterns of gene expression that typify and distinguish these disparate populations, however, are mostly unknown. Here, using adult mice of mixed sex, we combined retrograde labeling of supraspinal cell nuclei with fluorescence-activated nuclei sorting and single-nuclei RNA sequencing analyses to transcriptionally profile neurons that project axons from the brain to lumbar spinal cord. We identified 14 transcriptionally distinct cell types and used a combination of established and newly identified marker genes to assign an anatomic location to each. To validate the putative marker genes, we visualized selected transcripts and confirmed selective expression within lumbar-projecting neurons in discrete supraspinal regions. Finally, we illustrate the potential utility of these data by examining the expression of transcription factors that distinguish different supraspinal cell types and by surveying the expression of receptors for growth and guidance cues that may be present in the spinal cord. Collectively, these data establish transcriptional differences between anatomically defined supraspinal populations, identify a new set of marker genes of use in future experiments, and provide insight into potential differences in cellular and physiological activity across the supraspinal connectome.
SIGNIFICANCE STATEMENT The brain communicates with the body through a wide variety of neuronal populations with distinct functions and differential sensitivity to damage and disease. We have used single-nuclei RNA sequencing technology to distinguish patterns of gene expression within a diverse set of neurons that project axons from the mouse brain to the lumbar spinal cord. The results reveal transcriptional differences between populations previously defined on the basis of anatomy, provide new marker genes to facilitate rapid identification of cell type in future work, and suggest distinct responsiveness of different supraspinal populations to external growth and guidance cues.
Introduction
The supraspinal connectome comprises a diverse and widely distributed set of neurons that project axons to spinal targets and which convey a wide range of motor, autonomic, and sensory modulatory commands. Work spanning more than a century has elucidated the anatomy of supraspinal projections in various model organisms, and in the mouse our recent work has provided a unified description of the supraspinal connectome in three-dimensional space (Wang et al., 2022). In contrast to the state of anatomic and physiological knowledge, much less is understood about the patterns of gene expression that characterize supraspinal neurons or which may distinguish different subtypes. Establishing the transcriptional identities of supraspinal cell types is foundational to understand descending communication from the brain to spinal cord. In addition, the transcriptional signatures of the diverse supraspinal cell types provide essential baseline information to interpret the response of supraspinal populations to injuries and disease states that impact these populations (Blackmore et al., 2021).
Single-cell RNA sequencing (scRNA-seq) technologies offer unprecedented insight into cellular diversity by profiling transcript abundance in individual cells and by classifying them based on multidimensional indices of similarity (Armand et al., 2021). In the murine nervous system, single-cell datasets have identified ever finer distinctions between subtypes of neurons within numerous regions, including spinal cord, retina, sensory ganglia, and cortex (Tran et al., 2019; Renthal et al., 2020; Russ et al., 2021; Yao et al., 2021). Apart from recent datasets in corticospinal tract (CST) neurons, however, little is known about the transcriptional profiles that characterize neurons that project axons to different spinal targets, and the degree to which anatomic distinctions between the locations of supraspinal cell bodies are associated with distinct patterns of gene expression (Golan et al., 2021; Sahni et al., 2021b).
Here we combined retrograde tracing from lumbar spinal cord axons with single-nuclei RNA sequencing (snRNA-seq) to transcriptionally profile diverse types of supraspinal neurons. We identified 14 discrete classes of supraspinal neurons and identified marker genes that distinguish each class from other supraspinal types. To link these transcriptional clusters to anatomically defined populations, we cross-referenced candidate marker genes with anatomically registered transcript expression in the Allen Brain Atlas (https://mouse.brain-map.org/). These correspondences were further validated using ISH detection of marker transcripts in retrogradely labeled supraspinal neurons. Finally, differential gene expression analyses between supraspinal cell types revealed differences in the expression of transcription factors (TFs) and in receptors for growth and guidance cues, hinting at mechanisms that maintain cellular identity and potential differences in responses to extracellular cues between populations. Overall, these data provide new insight into patterns of gene expression that typify and distinguish diverse classes of spinally projecting neurons.
Materials and Methods
Plasmid construction and cloning
DNA encoding mScarlet (Bindels et al., 2017) was synthesized (Genscript) and fused in frame without a linker to human H2B (H2BC11) (accession #NM_021058) and cloned into the pAAV-CAG-tdTomato (Addgene #59462) using the sites KpnI and EcoRI at the 5′ and 3′ end, respectively. AAV2-retro-H2B-mScarlet was made by the University of North Carolina Viral Vector Core with a titer of 4.3 × 1012.
Spinal cord injections
All animal testing and research were conducted in compliance with ethical regulations laid out by the National Institutes of Health's Guide for the care and use of laboratory animals, and all experimental protocols involving animals were approved by the Institutional Animal Care and Use Safety committee at Marquette University (protocol #AR-314). Mice were bred and raised under a 24 h light–dark cycle with 12 h of light and 12 h of darkness. Ambient temperature was maintained at 22 ± 2°C and humidity between 40% and 60%. All 6 mice were female. AAV2-retro particles (1 μl) were injected into the spinal cord with a Hamilton syringe driven by a Stoelting QSI pump (catalog #53311) and guided by a micromanipulator (pumping rate: 0.04 μl/min). Adeno-associated virus (AAV) viral particles were injected to L1 spinal cord, 0.35 mm lateral to the midline, to depths of 0.6 and 0.8 mm, with 0.5 μl delivered to each site.
Dissection, FACS, and library preparation
Animals were anesthetized with isoflurane and promptly decapitated. The brain and brainstem were dissected out and placed in ice-cold slushy aCSF (Hearing et al., 2013) for 1 min. Brains were then sectioned in the sagittal plane at 500 μm intervals using Adult Mouse Brain Slicer Matrix on ice (Zivic Instruments BSMAS005-2). A total of nine sections were created from each animal (six replicate experimental animals plus four additional for FACS optimization and other pilot studies), including a midline section and four sections moving lateral into each hemisphere. Retrogradely labeled neurons were then microdissected using a stereomicroscope and fluorescence adapter (Nightsea SFA-GR). The brain regions collected from each section were recorded for future reference. The three sections from the left hemisphere and one midline section were dissected and flash frozen in a 1.5 ml Eppendorf DNA LoBind Microcentrifuge tube. The three sections from the right hemisphere were then collected and flash frozen with identical conditions, as mentioned above. Samples were then stored at −80°C until FACS.
To prepare cell nuclei on the day of library preparation, frozen tissue was promptly transferred from −80°C storage to a chilled 2 ml Dounce with chilled Nuclei EZ Lysis Buffer (Sigma-Aldrich N3408). Sample was Dounced 25 times with pestle A and 20 times with pestle B. The Dounced sample was then transferred to a chilled 15 ml conical tube and an additional 2 ml of Nuclei EZ Lysis Buffer was added and gently mixed by inversion. The sample incubated on ice for 5 min and then was centrifuged at 500 × g at 4°C. The supernatant was removed, and the pellet was resuspended in 4 ml of Nuclei EZ Lysis Buffer. The sample again was incubated on ice and centrifuged at 500 × g at 4°C for 5 min. Following the centrifugation, the supernatant was removed, and the pellet resuspended in 500 μl of Nuclei Suspension Buffer (2% BSA, 40 U/μl RNase Inhibitor; Invitrogen Ref. AM2684, 1× PBS). The resuspended solution was then filtered through a 20 μm filter and moved directly to the FACS machine. The collection tube for the sorted nuclei was coated with 5% BSA and contained the 10× Genomics master mix. We initially gate out debris using side scatter area (SSC-A) versus forward scatter area (FSC-A) and forward scatter width (FSC-W) versus forward scatter area (FSC-A). To eliminate doublets, the nuclei that pass through the gates above were then gated by side scatter width (SSC-W) versus side scatter height (SSC-H). Nuclei that passed through all the above gates were then filtered by level of fluorescence so that only the brightest were collected. With the above parameters, the FACS machine was set to collect 5000 nuclei in 5-7 min. The collected nuclei were then prepared into libraries, using Chromium Next GEM Single Cell 3′ Reagent Kits version 3.1 (PN-1000269), according to the manufacturer's protocol (10× Genomics, CG000204 Rev D).
Experimental design and statistical analyses
Libraries were sequenced with an Illumina NovaSeq 6000 at University of Wisconsin–Madison Biotechnology Center and then processed with CellRanger using default parameters to produce a Unique Molecular Identifier matrix for all nuclei-containing droplets, which were then analyzed in the Seurate version 4.1.0 R package. Median reads per cell were 125,600, 237,677, and 124,791; and mean unique genes (nFeature) per cell were 3190, 5017, and 4366 for Samples 1, 2, and 3, respectively. Nuclei with <1000 unique genes were excluded, as were nuclei with high levels of Atf3 and Creb5, indicative of cellular stress that likely resulted from the dissection procedure. The three samples were merged using the 2000 most variable genes as input for the “anchor.features” of the FindIntegrationAnchors() function. Clustering was performed according to Seurat recommendations based on 30 principal components as selected by the elbow plot heuristic and using FindNeighbors() and FindClusters() functions. Marker genes for each cluster were identified using FindAllMarkers() with default parameters, which utilizes a Wilcoxon rank-sum test comparing gene expression of cells within a given cluster versus all other cells. To identify variable TFs, we first extracted from the unified dataset a list of all TFs and associated average RNA counts based on TFs as identified by Lambert et al. (2018). For each TF, we then calculated a normalized expression value within each cluster by dividing that cluster's average RNA count value by the sum of all clusters and multiplying by 100 (i.e., in a hypothetical case in which a gene is expressed at the same level in all 14 clusters, each would receive a normalized score of 100/14 = 7.14.) Variable TFs were defined as scoring <1 (underrepresented) or >20 (overrepresented) in any cluster. Scripts used to analyze data are available (https://github.com/RegenerationLab/Supraspinal-Uninjured-Lumbar), and raw data are available through the NCBI Gene Expression Omnibus (Accession GSE212409).
ISH, imaging, and quantification
ISH was performed using RNAscope Multiplex Fluorescent Detection Kit version 2 from ACDbio (catalog #323110). Mouse brains were dissected, fresh frozen in OCT (VWR catalog #95057-838), and cryo-sectioned at 30 μm. Cryosection slides were dried at room temperature for 1 h before storage at −80°C. Slides were removed from −80°C and promptly began the ISH protocol according to the manufacturer's protocol. Slides were initially fixed in fresh 4% PFA for 1 h at room temperature and washed with 1× PBS (Alfa Aesar, J62036) twice for 2 min each. The slides were then dehydrated in 50%, 70%, and 100% ethanol for 5 min each and 100% repeated. Slides were allowed to dry at room temperature, followed by addition of hydrogen peroxide for 10 min. Slides were then washed with nuclease-free water, twice for 2 min each. We then applied Protease III (ACDbio catalog #322337) for 8 min at room temperature and washed with 1× PBS (Alfa Aesar, J62036) twice for 2 min each. Probes were then applied for 2 h at 40°C and washed with 1× wash buffer (ACDbio catalog #320058) twice for 2 min each. We then performed three consecutive amplification steps, washing with wash buffer in between, as instructed by the manufacturer. The corresponding probe HRP signal was then developed. All probes were detected using a 1:750 dilution of TSA plus fluorescein 488 (Akoya Biosciences NEL741001KT). Slides were then dried overnight. All probes were ordered from ACDbio (catalog numbers provided below). Probes used were Plagl1 (catalog #462941), Ttc6 (catalog #1125881), Emx2 (catalog #319001), Prdm6 (catalog #456891), Lhx4 (catalog #1130251), Pard3b (catalog #832241), Sox14 (catalog #516411), Kit (catalog #314151), and Col19a1 (catalog #539701).
Following ISH, IHC was used to amplify the retrograde label. Sections were washed in 1× PBS for 5 min and then blocked (10% serum, 1% PBST) for 1 h at room temperature. The primary antibody (2%-3% serum, 0.4% PBST, 1:500 Rockland 600-401-379) was applied and incubated for 90 min at room temperature. Sections were then washed with 1× PBS 3 times, 5 min each. Secondary antibody (2%-3% serum, 0.4% PBST, 1:500 Invitrogen A11035, 1:500 DAPI) was applied and incubated for 1 h at room temperature. Sections were then washed with 1× PBS (Alfa Aesar, J62036) 3 times, 5 min each, then mounted on slides with mounting media (Fluoro-Gel with Tris Buffer catalog #17 985-10).
Slides were imaged within 2 weeks from completion of ISH/IHC. Slides were imaged with Keyence BZ-X810 microscope and 10× Plan Apochromat Objective (BZ-PA10). Additional images were taken with Nikon AR1+ Confocal Microscope and 60× Apochromat Oil DIC N2 objective (MRD71600, NA 1.4); 60× images were gathered in a single z plane. Quantification was performed manually on digital images by first identifying all supraspinal cell nuclei in anatomic ROIs, as defined by expression of H2B-mSc, and then counting the number of visible puncta from RNAscope signal that overlapped with each nucleus. Nuclei were scored as positive if three or more puncta colocalized with nuclear signal. Images were obtained from 3 replicate animals, and a minimum of 60 individual cell nuclei were quantified from each region.
Results
Single-nuclei analysis of supraspinal neurons reveals transcriptionally distinct cell types
To label supraspinal neurons, adult mice received lumbar injection of AAV2-retro expressing mScarlet fluorophore that was localized to the nucleus by fused histone protein 2B (AAV2-retro-H2B-mSc) (Fig. 1A). We have shown previously that this procedure labels tens of thousands of supraspinal neurons distributed through forebrain, midbrain, and hindbrain (Wang et al., 2018, 2022). Two weeks after injection, animals were killed, a sagittal series of brain slices were rapidly prepared and placed in ice-cold aCSF and observed under fluorescent light. Based on anatomic location, we preliminarily identified fluorescent cells in cortex, hypothalamus, midbrain, dorsal pons, and hindbrain. Regions with fluorescent label were rapidly microdissected and flash frozen, taking care to exclude as much unlabeled tissue as possible. For each animal, the anatomic location of microdissected tissue from all sections was recorded. Tissue from 2 animals were pooled for each sample, from which cell nuclei were extracted and processed by fluorescent activated nuclei sorting (FANS) to purify the supraspinal subset. In initial studies, nuclei were sorted into PBS and inspected visually to confirm mScarlet signal in >98% (Fig. 1C). In subsequent experiments, nuclei were sorted directly into RT Master Mix followed immediately by GEM formation and library preparation according to 10× chromium instructions. Pilot studies revealed that library quality was highly dependent on rapid sorting, such that run times >10 min led to unacceptable levels of ambient RNA in the final sample. Thus, to minimize run times, it was essential to initiate FANS sorting with microdissected samples highly enriched for labeled nuclei. Five thousand nuclei were gathered from each sample; and of these, based on the volumes transferred to GEM formation, an estimated 3400 nuclei entered initial library creation. This procedure was repeated 3 times, yielding three independent replicates, each derived from 2 pooled animals. Libraries were sequenced to a minimum depth of 100 million reads, followed by clustering and analysis by Seurat (Butler et al., 2018).
Retrograde labeling and FANS purification of supraspinal nuclei. A, AAV2-retro-H2B-mScarlet was injected to lumbar spinal cord, followed 2 weeks later by microdissection of fluorescently labeled brain regions in sagittal section. Middle panels, Schematics represent the anatomic ROIs. Right, Tissue slices with retrograde fluorescent label. Red lines indicate regions of microdissection. B, Representative image of gates used for isolation of mScarlet-labeled nuclei. C, Visual inspection of nuclei after FANS shows mScarlet label in nearly all nuclei. LHA, Lateral hypothalamic area; PVH, paraventricular hypothalamus; MB, midbrain; DP, dorsal pons; Ret, medullary reticular formation. Scale bars: A, 1 mm; C, 50 µm.
In initial quality control steps, clusters of nuclei were identified with discrepantly low feature counts and with high levels of stress-related transcripts, such as Atf3 (Hai et al., 1999). These ∼1000 nuclei were presumed to be responding to damage or loss of nuclear membrane integrity during sample preparation and were removed from subsequent analysis. The remaining nuclei, which numbered 7748 across the three samples, were reclustered to yield 18 initial groups (Fig. 2A). Clusters were highly consistent between the three samples, an important indication that all samples contained a consistent set of putative cell types and that no cluster was the result of artifacts that were specific to any one sample (Fig. 2B). Based on the method of long-distance retrograde labeling, we predicted that sorted nuclei should derive from long distance projection neurons and not from non-neuronal cell types. Indeed, all clusters expressed high levels of neuronal markers Syt1 and Rbfox3/NeuN (Duan et al., 2016; Park and Ryu, 2018) and at most trace amounts of various markers from non-neural cells, including oligodendrocytes, astrocytes, microglia, and ependymal cells (Eng and Ghirnikar, 1994; Marques et al., 2016; Chen et al., 2017; Konishi et al., 2017; Galland et al., 2019; Sock and Wegner, 2021) (Fig. 2C). We examined markers for neurotransmitter subtypes, including glutamatergic (either Vglut1/Slc17a7 or Vglut2/Slc17a6) (Herzog et al., 2001; Kaneko and Fujiyama, 2002), inhibitory (Gad2) (Erlander et al., 1991), serotonergic (Tph2) (Hendricks et al., 1999; Ren et al., 2019), or noradrenergic neurons (Slc2a6) (Mulvey et al., 2018), all of which are known to project from brain to lumbar spinal cord. All 18 clusters were marked by expression of one of these transmitters, and none displayed expression of multiple transmitters (Fig. 2C,D). Interestingly, nearly all clusters expressed either Slc17a6 or Slc17a7, while only three very small clusters expressed the inhibitory, serotonergic, or noradrenergic markers: these nonglutamatergic clusters comprised <3% of the analyzed nuclei (Fig. 2D,E). This level of glutamatergic enrichment is likely disproportionate to brain-lumbar projection types, as inhibitory neurons may comprise more than one-third of brainstem-spinal projection neurons (Holstege, 1991). One likely explanation is that the tropism of AAV2-retro favors glutamatergic neurons over others (Tervo et al., 2016; Wang et al., 2018; Ganley et al., 2021). In summary, retrograde viral labeling and FANS produces a high purity neuron population, predominantly glutamatergic, which readily segregates into transcriptionally distinct clusters.
snRNA-seq analysis indicates that samples are mostly glutamatergic neurons. A, UMAP clustering of 7748 supraspinal nuclei identifies 18 distinct groups. B, Separation of the three merged samples confirms that each contributes to all 18 clusters, supporting replicate consistency. C, Dotplot represents expression of neuronal markers in all clusters and very low expression of non-neuronal markers, as expected from the retrograde labeling strategy. D, E, Violin plots (D) and feature plots (E) represent expression of glutamatergic markers in most clusters, and expression of GABAergic, serotonergic, or noradrenergic markers in only small clusters.
As a first step in classifying supraspinal populations, we sought to identify CST neurons, a prominent and relatively well-characterized projection from layer V of cortex to lumbar spinal cord. Indeed, a large and distinctly positioned set of cells displayed high expression of layer V markers genes Bcl11b, Crym, and Fezf2 (Arlotta et al., 2005; Greig et al., 2013; Fink et al., 2015) (Fig. 3A-C). Examination of transcripts that were enriched in this group revealed additional transcripts localized to layer V cortex according to in situ data from the Allen Brain Atlas, including Bcl6, Pdlim1, Cacna1h, Mylip, and Kcng1 (Fig. 3G–K). This group also expressed Slco2a1, a marker for layer Vb neurons that project to brainstem and spinal cord, but not Slc30a3, a marker for intracortical projection neurons or Npsr1 and Hpgd, markers for layer Vb neurons that project to the thalamus (Economo et al., 2018; Zhang et al., 2021). As an additional check, we also compared transcript abundance in this cluster to recent preprint data that used a single-cell approach to examine differential gene expression between lumbar- and cervically-projecting CST neurons (Golan et al., 2021). An important point of distinction is that our data included only lumbar-projecting neurons, preventing a direct comparison. It is notable, however, that the prior study reported Bcl6, Tox2, Rreb, Slc16a2, Bclr01, and Atp2b4 as enriched in lumbar-projecting CSTs; and consistent with this, we detected expression of each in the putative CST cluster. Interestingly, however, these transcripts were broadly expressed in other supraspinal clusters described below, indicating that, although these markers were reported to distinguish subtypes of CSTs from one another, they may not differentiate lumbar-projecting CSTs from other lumbar-projecting cell types (for average expression of all transcripts across all clusters, see Extended Data Fig. 9-1). Overall, we conclude that Seurat clusters 0 and 2, which together contain 3823 nuclei (49.9% of the total gathered), correspond to CST neurons. Another small cluster contained 88 nuclei and displayed high expression of general markers for cortex (Slc17a7, Satb2, Camk2a) but low levels of layer Vb markers and high levels of Slc30a3, a recently identified marker of intrahemispheric cortical neurons (Zhang et al., 2021). This cluster likely reflects trace contamination (∼2.2%) by non-CST neurons in the cortical sample and was removed.
Corticospinal neurons identified by established layer V markers. A-C, Feature plots of Bcl11b (A), Crym (B), and Fezf2 (C) show high enrichment in Clusters 0 and 2 (red circles). D-K, ISH detection from the Allen Brain Atlas of transcripts with enrichment in layer V of cortex (red arrows). Corresponding violin plots indicate high expression of transcripts in Clusters 0 and 2 (dotted red rectangles), consistent with corticospinal identity.
We next aimed to identity the remaining 15 Seurat clusters, which presumably corresponded to various supraspinal populations of subcortical origin, most of which lack established markers. Starting from transcripts that were highly enriched in each cluster, we adopted a strategy of manual curation that compared ISH data from the Allen Brain Atlas to the known locations of supraspinal neurons. This approach was aided by our recently created 3D atlas of supraspinal populations (Wang et al., 2022), which registers retrogradely labeled cell nuclei to a digital neuroanatomical atlas based on the Allen Brain Atlas. In this way, we could systemically examine locations in Allen Brain Atlas images that we knew to harbor supraspinal neurons, in search of putative marker genes from the scRNA-seq data. Below we outline the results of this initial classification in approximate rostral-to-caudal sequence.
Two clusters likely derived from distinct regions of the hypothalamus, a known source of supralumbar input (Hosoya, 1980). The first expressed high levels of Sim1, a known marker for the paraventricular hypothalamic region (Michaud et al., 1998) (Fig. 4A–D). The second expressed high levels of Plagl1, which in Allen Brain data are expressed strongly in the lateral hypothalamus (LH) (Fig. 4E–H). Next, a prominent cluster selectively expressed Ttc6, which we found to be highly expressed in the red nucleus in Atlas images (Fig. 4I–L). Cartpt is a known marker for the Edinger–Westphal nucleus (EW) (Xu et al., 2014; Topilko et al., 2022), and was selectively expressed in a small cluster (Fig. 4M–P). Another small cluster expressed Slc6a2, a well-established marker for noradrenergic neurons, and thus likely corresponds to the locus coeruleus (Mulvey et al., 2018) (Fig. 4Q–T). As noted previously, it is possible that the small size of this cluster may not reflect the true abundance of this projection, but rather the predominantly glutamatergic tropism of AAV2-retro. Finally, a larger group of cells selectively expressed the TF Prdm6, a transcript that localized in the Allen Brain Atlas to the dorsal pontine area (Fig. 4U–X). Our prior registration data showed a prominent cluster of supralumbar neurons in this region, which are distributed across several adjacent subnuclei, such as Barrington's nucleus. Consistent with this, the current data showed expression of Crh, a known marker for Barrington's nucleus (Verstegen et al., 2017), in a subregion of the Prdm6 cluster. Accordingly, we preliminarily assigned the Prdm6 cluster, which likely contains Barrington's nucleus among others, to the dorsal pontine region.
ISH data from the Allen Brain Atlas enables preliminary identification of transcriptional clusters. A, E, I, M, Q, U, Violin plots (left) and feature plots (right) of transcripts with high enrichment in single clusters. B, F, J, N, R, V, ISH data from the Allen Brain Atlas. Red boxes represent ROIs. C, G, K, O, S, W, Corresponding anatomic registration to regions that are known to supply supraspinal input. D, H, L, P, T, X, Higher-magnification views with detection of putative marker genes. Arrows point to areas of specific gene expression: LHA, Lateral hypothalamic area; PVH, paraventricular hypothalamus; RN, red nucleus; LC, locus coeruleus; PCG, pontine central gray. Scale bars: V, 1 mm; W, X, 200 µm.
The remaining nine subcortical clusters, which comprised ∼34% of all cells, were marked by expression of Hox3 and Hox4 gene clusters (Krumlauf and Wilkinson, 2021) (Fig. 5A,B). The Hox genes are important for patterning and segmentation of rhombomeres that give rise to the medulla and pons and serve as canonical markers of brainstem neurons (Chambers et al., 2009). Interestingly the specific Hox gene clusters (Fig. 5A) are mostly expressed in the caudal part of the developing hindbrain, suggesting that these clusters are anatomically located in the medulla. Consistent with this, most of the Hox+ clusters also expressed the homeobox protein Lhx4, previously noted for expression in reticulospinal neurons of brainstem origin (Cepeda-Nieto et al., 2005; Bretzner and Brownstone, 2013) (Fig. 5C). Also prominent in the Hox+ clusters was another homeobox gene, Vsx2/Chx10, previously detected in some populations of brainstem-spinal neurons (Fig. 5D) (Usseglio et al., 2020). We therefore concluded that these nine clusters likely derived from the caudal part of the brainstem. The Hox+ group also contained two small clusters, noted above, that expressed markers for inhibitory or serotonergic subtypes which likely correspond to trace amounts of brainstem GABAergic and raphe-spinal populations that took up AAV2-retro. In general, the remaining cell clusters that expressed Hox genes were not as sharply segregated as other supraspinal populations. Nevertheless, some transcripts were found to display clear concentration in subregions of the putative brainstem populations. These included Daam2 (Fig. 5E), Pard3b (Fig. 5F), Sox14 (Fig. 5G), Col19a1 (Fig. 5H), and Kit (Fig. 5I). Examination of data from the Allen Brain Atlas supported expression of these markers in brainstem regions, but unlike clusters assigned to the midbrain or forebrain, did not reveal strong segregation to annotated subregions of brainstem.
Identification of hindbrain populations. A, Dotplot of Hox transcript levels represents selective expression in a subset of HOX gene clusters, indicating hindbrain populations of caudal rhombomeric origin. B-I, Feature plots with putative hindbrain populations outlined and specific transcripts enriched in subregions indicated by red arrows.
Figure 6 summarizes assignments of identity to 14 cell populations, some of which are aggregates of original Seurat clusters, and lists “marker” transcripts that were concentrated in each. Supraspinal populations rostral to the hindbrain appear to be transcriptionally distinct, as evidenced by widely separated clustering via Uniform Manifold Approximation and Projection (UMAP) (Fig. 6A) and the presence of multiple transcripts in each that are highly specific (Fig. 6B). In contrast, distinctions between hindbrain populations are generally less definitive, with less separation by UMAP and fewer distinctive marker genes. As such, the hindbrain was divided into five populations, the first four of which were based on well-segregated markers (Daam2, Pard3b, Col19a1, and Sox14), and the fifth derived from a disparate collection of clusters, none of which displayed highly distinctive markers. Overall, single-nuclei profiling identified a total of 13 types of supraspinal neurons that are readily distinguishable, along with an additional set of neurons that derive from hindbrain regions, but which were not clearly distinguishable in the current dataset.
Anatomical assignment and marker genes for populations of supraspinal neurons. A, UMAP clustering of supraspinal nuclei with anatomic designations for each. B, Dotplot of gene expression showing examples of transcripts with enriched expression in each cluster relative to other supraspinal populations. HB, Hindbrain; PVH, paraventricular hypothalamus; RN, red nucleus; LC, locus coeruleus; DP, dorsal pons.
ISH verifies expression of putative marker genes in supraspinal populations
To validate the classification of populations outlined above, we combined retrograde labeling of lumbar-projecting supraspinal neurons with visualization of transcripts using ISH. As previously, adult mice received injection of AAV2-retro-H2B-mScarlet to lumbar spinal cord followed 2 weeks later by ISH of brain sections. ISH was performed using ACDbio probes for specific marker genes, followed by imaging to assess colocalization of mScarlet+ supraspinal nuclei with in situ signal. This analysis found strong concordance of the pattern of transcript detection with predictions from the single-nuclei data. For example, Plagl1 transcript was detected in supraspinal neurons located in the LH (Fig. 7A,B). Two candidate markers for the red nucleus, Ttc6 and Emx2, were both detected in rubrospinal neurons (Fig. 7D,E), Prdm6 was detected as predicted in supraspinal neurons clustered in the dorsal pons (Fig. 7G,H,J,K), and Lhx4 was detected as predicted in hindbrain supraspinal neurons (Fig. 7M,N). Importantly, inspection of other supraspinal populations within the same tissue sections showed an absence of transcript detection, indicating selective expression in the predicted region (Fig. 7C,F,I,L,O). Indeed, quantification revealed that, within expected regions, nearly all supraspinal neurons colocalized with signal for predicted transcripts (range 83.3%-98.9%), whereas in other regions almost no supraspinal neurons showed colocalization (range 0%-4.6%) (Fig. 7P). These data substantiate the ability of the single-nuclei data to identify marker transcripts that are ubiquitous within an anatomically defined supraspinal population and which are also highly selective for that population. It should be noted that, for all the probes tested, we also observed positive signal in nonsupraspinal brain regions. Thus, these markers can be considered to be specific within the context of the supraspinal connectome, but not in the broader context of the nervous system.
Visualization of candidate marker genes in subcortical supraspinal populations. Adult mice received lumbar injection AAV2-retro-H2B-mScarlet followed 2 weeks later by tissue sectioning and ISH detection of candidate marker genes. In all panels, red represents retrograde H2B-mSc; green represents transcript detection. A, D, G, J, M, Brain sections with regions of predicted expression indicated in blue and a predicted nonexpressing region in orange. B, E, H, K, N, Higher-magnification view of target populations with insets. C, F, I, L, O, Higher magnification of nonexpressing supraspinal populations. P, Percent of supraspinal cells in which transcript was detected, indicating selective expression of Plagl1 in the lateral hypothalamic area, Ttc6 Emx2 in the red nucleus, and Prdm6 in the dorsal region of the pons. LHA, Lateral hypothalamic area; HB, hindbrain; SN, solitary nucleus; DP, dorsal pons. Scale bars: M, 2 mm; N, 50 µm; N′, O, 10 µm. N = 3 animals; >60 cells per region. Error bars are SEM.
In hindbrain regions, the correspondence between anatomic division and marker expression was more complex. As expected, ISH confirmed expression of Pard3b (Fig. 8A–C), Sox14 (Fig. 8D–F), Kit (Fig. 8G–I), and Col19a1 (Fig. 8J–L) in supraspinal neurons in the hindbrain region. Col19a1, Sox14, and Pard3b were detected in supraspinal neurons that spanned the gigantocellular reticular and magnocellular reticular nucleus (Fig. 8A,D,J), whereas Kit appeared more laterally, likely including the paragigantocellular reticular nucleus, lateral part (Fig. 8G). Importantly, however, marker-positive supraspinal neurons were interspersed with nonexpressing supraspinal neurons (Fig. 8C,F,I,L: compare positive nuclei marked with arrows to adjacent negative nuclei marked with arrowheads). To quantify this phenomenon, for each probe, we selected medulla sections from 3 replicate animals with positive signal and then quantified the percent of total supraspinal neurons in that section, identified by retrograde mSc, with detectable expression. We found expression of Pard3b in 34.5% (±1.9% SEM), Col19a1 in 10.1% (±1.8% SEM), Sox14 in 25.7% (±3.8% SEM), and Kit in 5.9% (±3.2% SEM). Importantly, because the entirety of the hindbrain was not sampled, these values are not estimates for hindbrain-wide expression but rather illustrate the point that, even within discrete sample planes, markers are expressed in subsets of cells that are intermixed with nonexpressing cells. Overall, the ISH results support the identification of markers that distinguish broad classes of anatomically defined supraspinal neurons and support the assignments made previously in Figure 6, but indicate that, within the hindbrain region, the transcriptional distinctions in cell type do not map cleanly onto existing anatomic classification.
Visualization of candidate marker genes in hindbrain supraspinal populations. Adult mice received lumbar injection AAV2-retro-H2B-mScarlet followed 2 weeks later by tissue sectioning and ISH detection of candidate marker genes. A, D, G, J, Overview of brain sections with medullary ROIs (arrows). B, E, H, K, Higher-magnification views of target regions. C, F, I, L, High-resolution images of supraspinal cell nuclei (red) with transcript detection (green) in some neurons (arrows) but absent in nearby supraspinal cells (empty arrowheads). Scale bars: A, D, G, J, 1 mm; B, E, H, K, 50 µm; C, F, I, L, 20 µm.
TF expression, receptor/ligand profiles, and voltage-gated ion channels in supraspinal populations
Diverse families of TFs are essential for neuronal identity and connectivity during development and postnatal periods (Russ and Kaltschmidt, 2014; Perreault and Giorgi, 2019). To gain insight into the regulation and maintenance of different supraspinal subtypes, we first queried the data for differences in the expression of TFs. Starting from a curated list of all TFs (Lambert et al., 2018), we focused on TFs that displayed the greatest variability in expression between different supraspinal populations (see Materials and Methods). Figure 9A shows a set of 79 variable TFs, while Extended Data Figure 9-1 provides a spreadsheet with the average expression of all transcripts, including those identified as TFs, across all supraspinal cell types. As already noted, some TFs are highly specific to a single supraspinal type, such as Sim1 in paraventricular hypothalamic, Prdm6 in Dorsal Pons, Plagl1 in LH, and Sox14 in a subset of hindbrain neurons. Additional TFs with high specificity included Ebf2 in the midbrain midline nuclei (EW), Phox2a and -b in the Locus Coeruleus, Rfx4 in LH, and Rreb1 in the Red Nucleus. Numerous TFs differed broadly between CST and subcortical populations, with some factors highly enriched in cortex (Etv1, L3mbtl4, Satb2, and Zeb2) and others expressed in numerous subcortical populations but low or absent in CST (Cux2, Dach1, Glis1, Zfhx4, and Zbtb20). Other TFs were not specific to any one population but instead were common to a small subset and excluded from all others. One striking example of this is Tfap2b, which was found to be abundant in EW, LC, and a subtype of hindbrain but was not expressed elsewhere. Overall, these data provide foundational knowledge regarding the overlapping sets of TFs that typify single types and groups of supraspinal neurons.
Differential expression of TFs, growth factor, axon guidance, and growth inhibitory receptors in supraspinal populations. A, Dotplot represents the expression levels of selected TFs in each of 14 distinct types of supraspinal neurons. B, Dotplot represents the level of expression of selected receptors in supraspinal populations. Relevant ligands are listed to the left and include neurotrophins, axon guidance cues, and growth-inhibitory cues. C, Dotplot represents expression of voltage-gated ion channels. Extended Data Fig. 9-1 provides average expression of all transcripts in all populations.
We next considered the expression of receptors for growth factors and axonal guidance cues across supraspinal subtypes. This information is potentially informative in predicting differential responses to cues in the spinal cord environment. We first assembled a list of ligands of interest, based on neurotrophins, growth factors, and axon guidance cues that are either endogenous to spinal cord tissue or which have been exogenously applied in experimental paradigms of spinal injury (Hata et al., 2006; Giger et al., 2010; Goldshmit et al., 2011; Dun and Parkinson, 2017; Haenzi and Moon, 2017; Liu et al., 2017; Rosich et al., 2017; Badhiwala et al., 2018; Walker and Xu, 2018; Gao et al., 2019; Kitamura et al., 2019; Yamane et al., 2019). We then examined the expression of various receptors for each cue across different subtypes of supraspinal neurons (Fig. 9B). Many receptors displayed broad expression across all types of supraspinal populations (e.g., neurotrophin receptors Ntrk2/TrkB and Ntrk3/TrkC), and others were missing from all types (e.g., Ntrk1/TrkA). It is notable that receptors for GDNF, which was shown recently to stimulate growth from injured propriospinal neurons, are mostly absent from supraspinal populations in the brain, which is consistent with minimal levels of supraspinal regeneration that were reported (Anderson et al., 2018). Expression of some receptors was highly variable; for example, CST neurons expressed high levels of receptors for semaphorins, Npn1 and Plexa2, and high levels of Unc5d, which confers repulsive guidance signaling from netrin (Huber et al., 2005; Huettl et al., 2011; Helmbrecht et al., 2015; Alto and Terman, 2017; Dun and Parkinson, 2017). In contrast, the Dcc receptor that confers positive growth responses to netrin (Dun and Parkinson, 2017) as virtually absent from CST neurons but expressed broadly across other supraspinal types. In addition, some interesting differences between populations were detected. These data are consistent with the established responsiveness of CST neurons to semaphorin signaling and hint that CST neurons may potentially respond differently to netrin ligands than do other supraspinal populations. Important caveats to this approach include the fact that mRNA abundance does not necessarily scale with protein expression and that some receptors are known to be expressed by neurons but not transported into axons (Koseki et al., 2017). Thus, the presence of mRNA for a receptor is not definitive evidence of that the neuron's axon will respond to the corresponding ligand when presented in spinal tissue. Nevertheless, particularly in cases in which the mRNA for a receptor is not detected in a given supraspinal class, these comparative data provide the basis to rapidly generate hypotheses of differential response to ligands presented in the spinal cord in regeneration studies.
Finally, we examined the expression of voltage-gated ion channels, an important determinant of electrophysiological properties (Fig. 9C). Some channels (e.g., sodium channels Scna1 and Scna2, potassium channels Kcnb2 and Kcnd2, and L-type, N-type, and P-type calcium channels) were found to be broadly expressed across all supraspinal cell types. Others, however, showed more variable patterns of expression, with CST neurons emerging as possessing the most distinct channel profile. For example, compared with other supraspinal cell types, CST neurons expressed very low levels of the sodium channel Scn9a (Nav1.7) and potassium channel Kcnc2 (Kv3.2), but higher expression of the calcium-gated potassium channel Kcnn2 (SK2), sodium-gated potassium channel Kna1.2 (Slick), and, notably, T-type calcium channels (Cav3.2 and Cav3.3). To validate the channel profiles, we compared our channel detection in CST neurons to prior immunohistochemical or pharmacological inhibition/electrophysiology studies of layer V cortical neurons and found broad correspondence, corroborating expression of channels, including Nav1.1, Nav1.2, Nav1.6, Kv1.1, Kv1.2, Kv7.2, Kv7.3, Kna1.2, Cav3.2, and Cav3.2 (Kole and Stuart, 2012; Battefeld et al., 2014; Rizzi et al., 2016; Hu and Bean, 2018). Thus, these data provide a framework to generate hypotheses and to interpret potential differences in electrophysiological properties between different supraspinal cell types.
Discussion
We have identified 14 transcriptionally distinct populations of neurons that project axons from the murine brain to the lumbar spinal cord. In addition, using a combination of well-established and newly discovered marker genes, we have associated these transcriptionally distinct clusters with anatomically defined populations, thus linking new transcriptional insights to a rich neuroanatomical literature. These data establish new marker genes, offer insight into potential physiological differences between supraspinal cell types, and provide a population-by-population baseline for future study of the transcriptional impacts of injury and disease.
Neuronal diversity in the supraspinal connectome
scRNA-seq technologies are rapidly expanding awareness of cellular diversity in the nervous system. Single-cell and -nuclei data exist for various regions in the adult and developing mouse nervous system, including retina (Tran et al., 2019), DRG (Renthal et al., 2020), spinal cord (Russ et al., 2021), cortex (la Manno et al., 2021; Yao et al., 2021), and CST (Golan et al., 2021). Here we aimed to broaden understanding of supraspinal cell types spanning forebrain, midbrain, and hindbrain, and to test whether transcriptionally unique profiles associate with established anatomic classifications. Because supraspinal populations carry distinct functions and target distinct spinal circuits (e.g., motor vs autonomic), we expected to identify differences. Indeed, our findings broadly confirm the assumption that anatomically separated supraspinal neurons also differ from one another in their expressed transcripts. Major classes of supraspinal neurons existed as discrete UMAP clusters with some transcripts that were ubiquitous within that cluster but essentially absent elsewhere, a finding confirmed by ISH-based visualization and quantification. Thus, across major groups of the supraspinal connectome, these data enhance the collective understanding of neural diversity by linking the transcriptional state to the location of their cell bodies, and by extension to previously established understanding of physiological roles.
Hindbrain-spinal populations, however, stand as a partial exception. Although hindbrain neurons clustered separately from midbrain and forebrain and were recognizable by markers, such as Hox and Lhx factors, hindbrain subtypes remained closely adjacent by UMAP and mostly shared putative marker genes with at least one other cluster. These results are reminiscent of the ventral spinal cord, in which neuronal subtypes were not well delineated in early single-cell datasets (Sathyamurthy et al., 2018) but could later be distinguished in larger samples (Russ et al., 2021). Finer distinctions between hindbrain populations will likely emerge as additional datasets are created and aggregated. Interestingly, even for the few transcripts detected uniquely in one hindbrain cluster (e.g., Sox14, Pard3b, Col19a1, and Kit), ISH visualization showed marker-positive neurons to be distributed throughout the medulla and/or intermingled with marker-negative supraspinal neurons. These data support the presence of transcriptionally distinct subsets of medullary supraspinal neurons but indicate that they are not necessarily anatomically segregated or aligned with existing anatomic classifications. These data lay a foundation for future work to link transcription and neuroanatomy, for example, testing whether transcriptional subtypes share common projections, inputs, or functions (Usseglio et al., 2020; Ruder et al., 2021).
A foundation for cell-specific manipulation
In principle, starting from marker transcripts expressed exclusively in in one supraspinal cell type, it may be possible to recapitulate selective gene expression by identifying relevant promoter or enhancer sequences and incorporating them into vectors. This strategy has been used successfully to target layer V cortical neurons and specific inhibitory subclasses (Vormstein-Schneider et al., 2020; Graybuck et al., 2021). Importantly, as we have emphasized throughout, the markers described here are specific only within the context of the supraspinal connectome but can be found elsewhere in the brain. Thus, the preferred strategy would be to design vectors with cell type-specific enhancer elements and then deliver them to the spinal cord in retrograde fashion to avoid this off-target expression. Delivering optogenetic, chemogenetic, or Cre recombinase transgenes under the control of cell type-specific enhancers, either to WT animals or to transgenic animals with genes flanked by loxP sites, offers an attractive strategy to parse circuit functions or to selectively test gene function within individual supraspinal populations.
Implications for injury and diseases
Damage to the spinal cord disrupts many ascending and descending tracts that carry a wide range of motor, sensory, and autonomic functions. SCI research, however, has historically focused disproportionately on major motor pathways, such as the corticospinal, rubrospinal, and reticulospinal, resulting in more limited understanding of additional supraspinal pathways that carry functions also of value to individuals with spinal injury (Blackmore et al., 2021). scRNA-seq approaches now afford new opportunities to study a broader set of cell types impacted by spinal injury, in alignment with the needs of individuals with spinal injury. In this context, the current data serve as a needed foundation for at least three research directions.
First, the current data provide baseline information for future characterization of changes triggered by axotomy or disease. It has been known for decades that different supraspinal cell types display highly distinct regenerative responses after spinal injury and after treatment; for example, seminal work showed regeneration of brainstem-spinal axons, but not corticospinal, into peripheral nerve grafts (David and Aguayo, 1981; Richardson et al., 1984). One possibility is that underlying these different growth responses are variations in the transcriptional response to injury. The present data provide an initial classification of transcriptionally distinct supraspinal cell types and, critically, a population-by-population baseline to identify changes after damage.
Second, the present data can be mined to generate hypotheses regarding responsiveness to external cues. As one example, we have examined receptors for axonal growth and guidance cues, both endogenous and exogenous. We find that supraspinal populations generally express receptors for neurotrophins BDNF and NT-3 but not NGF, consistent with an extensive literature regarding their growth effects after spinal injury (for review, see Keefe et al., 2017). Supraspinal populations also generally lacked receptors for GDNF, which may explain its stimulation of axon growth from propriospinal but not supraspinal axons (Deng et al., 2013; Anderson et al., 2018). Interestingly, we also find that most CST neurons projecting to the lumbar segments express Crim1, a transmembrane protein that contributes to lumbar targeting in CST neurons (Sahni et al., 2021a), hinting at a wider role in directing other supraspinal axons to lumbar targets.
In contrast, some receptors were highly variable between supraspinal populations, hinting at differential responses to cues. For example, semaphorins are inhibitory cues expressed in injured spinal tissue (Shim et al., 2012; Ueno et al., 2020). CST neurons express high levels of semaphorin receptors Nrp1 and Plxna2, and interfering with semaphorin/Plxna2 signaling enhances postinjury sprouting and reduces axon retraction in CST neurons (Shim et al., 2012; Ueno et al., 2020). We find, however, that most other supraspinal populations expressed very low levels of semaphorin receptors, indicating that this inhibition may be specific to CST axons. Netrin-1, which can attract axons via the Dcc receptor or repulse axons via Unc5 receptors, is expressed by spinal oligodendrocytes after injury and has been studied both as an endogenous inhibitor of axon growth (Löw et al., 2008) and a potential reparative factor. Interestingly, CST neurons express very low levels of Dcc but high levels of Unc5a, hinting at a more negative response to Netrin than other supraspinal cell types.
Finally, an interesting future application for these data would be to harness existing knowledge of adhesive and synapse formation protein–protein interactions to build models of brain-spinal circuity. In principle, by examining differential cell-surface protein expression in the present data alongside similar data from diverse spinal cells, it may be possible to predict favored synaptic partnerships (Ximerakis et al., 2019; Sanes and Zipursky, 2020; Russ et al., 2021). Doing so could potentially clarify the complex descending connections and provide foundational knowledge for efforts to reconstruct them after injury.
Important caveats apply to these data. First, the proposed applications, particularly those involving ligand/receptor pairs, must consider the potential disconnect between mRNA abundance and protein expression and/or localization. Second, technical constraints excluded populations that were too small or too dispersed for rapid collection and sorting, notably the pontine reticular formation and a small population in the dorsal medulla near the solitary nucleus. A third caveat regards the tropism of AAV2-retro. The lumbar spinal cord receives substantial descending input from serotonergic, noradrenergic, and GABAergic populations, yet together they comprised <3% of transduced cells. The relative inefficiency of AAV2-retro in serotonergic and noradrenergic cell types has been noted previously (Tervo et al., 2016; Wang et al., 2018; Ganley et al., 2021). The tropism of AAV2-Retro likely favors glutamatergic neurons, highlighting the importance of developing additional retrograde vectors with wider tropism. Fourth, because neurons were traced from only the lumbar spinal cord, the data do not reveal subtype diversity associated with segmental targeting in the spinal cord, as in Golan et al. (2021). Expanded cell collection, larger sample sizes, and alternative methods of analysis might reveal additional neuronal subtypes within the main groups distinguished here. Nevertheless, the present data offer a significant advance regarding the transcriptional state of a wide range of neuronal populations that communicate with the lumbar spinal cord, offering new insights into potential physiological differences.
Figure 9-1
Average expression values of all detected transcripts in all clusters. The associated Excel file provides the average transcript abundance of all detected features for each of the 14 identified cell types. A column indicates whether each transcript is identified as a TF (Lambert et al., 2018) or voltage-gated ion channels. Download Figure 9-1, XLSX file.
Footnotes
This work was supported by the National Institute of Neurological Disorders and Stroke, the Bryon Riesch Paralysis Foundation, the Miami Project to Cure Paralysis, and the Buoniconti Fund. We thank Angela Schmoldt (University of Wisconsin–Milwaukee) for bio-analyzer assistance; the University of Wisconsin–Madison Biotechnology Center for sequencing; and James Choi, Dr. Dmitry Velmeshev, Dr. Kevin Park, and Dr. Jae Lee and Nirupa Chaudari for insight and critical evaluation of the work.
The authors declare no competing financial interests.
- Correspondence should be addressed to Murray G. Blackmore at murray.blackmore{at}marquette.edu
