Abstract
The corticospinal tract (CST) forms a central part of the voluntary motor apparatus in all mammals. Thus, injury, disease, and subsequent degeneration within this pathway result in chronic irreversible functional deficits. Current strategies to repair the damaged CST are suboptimal in part because of underexplored molecular heterogeneity within the adult tract. Here, we combine spinal retrograde CST tracing with single-cell RNA sequencing (scRNAseq) in adult male and female mice to index corticospinal neuron (CSN) subtypes that differentially innervate the forelimb and hindlimb. We exploit publicly available datasets to confer anatomic specialization among CSNs and show that CSNs segregate not only along the forelimb and hindlimb axis but also by supraspinal axon collateralization. These anatomically defined transcriptional data allow us to use machine learning tools to build classifiers that discriminate between CSNs and cortical layer 2/3 and nonspinally terminating layer 5 neurons in M1 and separately identify limb-specific CSNs. Using these tools, CSN subtypes can be differentially identified to study postnatal patterning of the CST in vivo, leveraged to screen for novel limb-specific axon growth survival and growth activators in vitro, and ultimately exploited to repair the damaged CST after injury and disease.
SIGNIFICANCE STATEMENT Therapeutic interventions designed to repair the damaged CST after spinal cord injury have remained functionally suboptimal in part because of an incomplete understanding of the molecular heterogeneity among subclasses of CSNs. Here, we combine spinal retrograde labeling with scRNAseq and annotate a CSN index by the termination pattern of their primary axon in the cervical or lumbar spinal cord and supraspinal collateral terminal fields. Using machine learning we have confirmed the veracity of our CSN gene lists to train classifiers to identify CSNs among all classes of neurons in primary motor cortex to study the development, patterning, homeostasis, and response to injury and disease, and ultimately target streamlined repair strategies to this critical motor pathway.
Introduction
The corticospinal tract (CST) is the major descending motor pathway responsible for driving fine coordinated movement in mammals (Lemon, 2008). Corticospinal neurons (CSNs) in layer (L)5b (L5) of sensorimotor cortex project their axons through the internal capsule, decussate in the medulla, and innervate every spinal segment along the neuroaxis. The CST is wired postnatally, with pioneer axons extending to spinal enlargements, followed by periods of exuberant gray matter terminal arborization and pruning via activity-dependent mechanisms during a protracted critical period that sculpts a mature motor pathway (Gianino et al., 1999; Martin, 2005). Complex wiring of the mature CST, together with its central role in voluntary and fine motor control, means that damage leads to significant and lasting functional impairments. Efforts to repair the damaged CST have focused on either nullifying the effects of the axon growth inhibitory environment of the mature CNS or recapitulating cell autonomous developmental mechanisms (Schwab and Strittmatter, 2014; Gutilla and Steward, 2016; He and Jin, 2016; Hilton and Bradke, 2017; Bradbury and Burnside, 2019; Fawcett, 2020). Although strides have been made in stimulating axotomized CSNs to regenerate or intact CSNs to undergo plasticity after injury, recovery of significant motor function remains to be realized (Park et al., 2008; Fink et al., 2017; Kauer et al., 2022). The inefficacy of current interventions could in part be because of an incomplete understanding of the molecular heterogeneity among CSNs. For instance, in vitro and in vivo screening approaches designed to identify axon growth modulators may overlook molecular heterogeneity within cortical neuron subclasses including CSNs, consequently diluting the most potent candidates among differentially sensitive neuronal subtypes (Blackmore et al., 2010; Buchser et al., 2010; Sekine et al., 2018).
Recent data from our laboratory supports differential sensitivity of CSN subdivisions. We showed that novel proplasticity factors identified in intact CSNs undergoing functional plasticity after unilateral pyramidotomy (PyX) stimulates growth of lesioned and intact forelimb (FL) but not hindlimb (HL) CST axons (Fink et al., 2017; Kauer et al., 2022). Proaxon growth candidates Lppr1 and Inpp5k were identified via retrograde labeling of sprouting CSNs in the cervical spinal cord, suggesting that retrograde labeling of plastic CSNs from the lumbar cord may identify a separate set of factors.
Single-cell RNA sequencing (scRNAseq) approaches offer the sensitivity to identify transcriptional heterogeneity among CSNs. Indeed, previous studies have revealed the rich phenotypic diversity among CNS cell types (Arlotta et al., 2005; Macosko et al., 2015; Saunders et al., 2018; Tasic et al., 2018; Milich et al., 2021; Yao et al., 2020; Muñoz-Castañeda et al., 2021). Central to the utility of exploiting transcriptional indices is understanding the functional or anatomic role of potential subclasses of cells. For instance, scRNAseq of corticocortical projection neurons showed that although these cells form a single genetic cluster, their anatomically traced subdivisions have diverse gene expression (Kim et al., 2020).
To identify CSN transcriptional heterogeneity, we combined retrograde adeno-associated virus (AAV) tracing from the cervical and lumbar spinal cord with scRNAseq of adult CSNs. We used machine learning to develop classifiers that go beyond cell marker characterization and provide robust tools for unbiased classification of CSNs among other neurons in primary motor cortex (M1) and among limb-specific CSNs. We exploited taxonomic definitions of projection neuron cell types in M1 to show that CSNs can be anatomically and molecularly defined based on supraspinal terminals in addition to their spinal projections. Critically, we present evidence showing that intact adult CSNs express a constellation of genes previously shown to be unique to nonspinally terminating projection neurons. Using intersectional viral tracing, we confirm that forelimb and hindlimb CSNs differentially innervate supraspinal structures, further supporting emerging data that CSNs influence motor output beyond directly modulating spinal motor circuitry (Nelson et al., 2021). We believe these sequencing data can be leveraged for in vitro and in vivo screening, targeting, and exploitation strategies to enhance our understanding of the development, patterning, housekeeping, and response to injury within the CST.
Materials and Methods
Mice
Three mouse lines were used for experiments—Rbp4 Cre mice (gift from David Berson, Brown University), Ai14 (The Jackson Laboratory), and C57BL/6J (The Jackson Laboratory). Rbp4 Cre mice were crossed with the cre-dependent tdTomato reporter line Ai14 to label deep-layer cortical neurons and to confirm the veracity of our CSN classifier (CSN-C).
Surgery
All procedures and postoperative care were performed in accordance with the guidelines of the Institutional Animal Use and Care Committee at Yale University.
Intraspinal retrograde AAV injections for scRNAseq
To complete retrograde labeling of cervical-projecting CSNs, adult mice (2–3 months old, n = 19, 10 males, 9 females) were anesthetized with ketamine (100 mg/kg; Covetrus) and xylazine (15 mg/kg; Covetrus) and placed in a stereotaxic frame (Stoelting). An incision was made over the cervical enlargement and the C5–C8 vertebrae revealed by blunt dissection of overlying muscle. A bilateral laminectomy was performed to expose the underlying C5–C8 spinal cord, and a small incision was made in the dura mater. The tip of a pulled glass capillary tube attached to a 5 µl Hamilton syringe loaded into a Micro4 infusion device (World Precision Instruments) was slowly inserted stereotaxically to a depth of 500 µm into the C6 level of the spinal cord and ∼600 µm lateral from the midline. Thirty seconds after the introduction of the capillary tube, 100 nl of retro-AAV-CAG-GFP (catalog #37825-AAVrg, Addgene) was infused into the spinal cord over 2 min. The tip was left in situ for an additional 30 s before removal.
This procedure was completed seven additional times bilaterally at the same coordinates at evenly spaced injection sites between C6 and C7 resulting in a total infusion of 800 nl (100 nl over eight sites). Muscle layers were sutured with VICRYL (Ethicon) and skin with monofilament suture. All animals received postsurgical antibiotics (Ampicillin, 100 mg/kg subcutaneously) and analgesia (Buprenorphine, 0.05 mg/kg subcutaneously) for 2 d after lesion. All animals recovered uneventfully.
To complete retrograde labeling of lumbar-projecting CSNs, adult mice (2–3 months old, n = 19, 10 males, 9 females) were anesthetized with ketamine (100 mg/kg) and xylazine (15 mg/kg) and placed in a custom-built spine stabilizer (Farrar et al., 2012). Using the last rib as a landmark, an incision was made, and a bilateral laminectomy was performed to expose L4–L5 spinal cord. The tip of a pulled glass capillary tube attached to a 5 µl Hamilton syringe loaded into a Micro4 infusion device was slowly inserted stereotaxically to a depth of 500 µm into the L4 level of the spinal cord and ∼500 µm lateral from the midline. Thirty seconds after introduction of the capillary tube, 100 nl of retro-AAV-CAG-GFP was infused into the spinal cord over 2 min. The tip was left in situ for an additional 30 s before removal. This procedure was completed seven additional times bilaterally at the same coordinates at evenly spaced injection sites between L4 and L5, resulting in a total infusion of 800 nl (100 nl over eight sites). All animals received postsurgical antibiotics (Ampicillin, 100 mg/kg subcutaneously) and analgesia (Buprenorphine, 0.05 mg/kg subcutaneously) for 2 d postlesion. All animals recovered uneventfully.
Injections for single molecule fluorescent in situ hybridization validation of scRNAseq
For confirmation studies, adult mice (2–3 months old, n = 3/probe) received bilateral injections of retro-AAV-CAG-GFP into the cervical cord and retro-AAV-CAG-tdTomato (catalog #59462-AAVrg, Addgene) into the lumbar cord as described above for a total of 16 spinal injections (100 nl over eight sites in the cervical enlargement, and 100 nl over eight sites in the lumbar enlargement). All animals received postsurgical antibiotics (Ampicillin, 100 mg/kg subcutaneously) and analgesia (Buprenorphine, 0.05 mg/kg subcutaneously) for 2 d postlesion. All animals recovered uneventfully.
Intersectional tracing of supraspinal CST terminals for BrainJ analysis
To label the FL and HL corticospinal tract for projection analysis, two groups of adult C57/blk6 mice received dual injections of pAAV-CAG-FLEX-tdTomato (catalog #28306-AAV1, Addgene) into either the caudal FL motor cortex (n = 4; AP, between 0.0 mm and 0.5 mm rostral to bregma; ML, between 1.0 mm and 1.3 mm lateral to bregma; four injections, 100 nl/site) or the HL motor cortex (n = 4; AP, between −0.5 mm and −1 mm caudal to bregma; ML, between 1.5 mm to 2.0 mm lateral to bregma; four injections, 100 nl/site), and pENN.AAV.hSyn.Cre.WPRE.hGH (catalog #105553-AAVrg, Addgene) into either the cervical (C6/7) or lumbar spinal cord (L4/5), as described above. All animals received postsurgical antibiotics (Ampicillin, 100 mg/kg subcutaneously) and analgesia (Buprenorphine, 0.05 mg/kg subcutaneously) for 2 d postlesion. All animals recovered uneventfully.
Intact adult single-cell isolation for scRNAseq
Adult intact cell isolation was done following our previously published protocol (Golan and Cafferty, 2021). Briefly, mice that received retrograde injections into either the lumbar or cervical cord were anesthetized with isoflurane and transcardially perfused with artificial CSF (aCSF) 10 d after retro-AAV injection. The aCSF contained the following (in mm): 0.5 CaCl2, 25 glucose, 96 HCl, 20 HEPES, 10 MgSO4, 1.25 NaH2PO4, 3 myo-inositol, 12 N-acetylcysteine, 96 NMDG, 2.5 KCl, 25 NaHCO3, 5 sodium l-ascorbate, 3 sodium pyruvate, 0.01 taurine, 2 thiourea, 13.2 trehalose, and was bubbled with carbogen gas (95% O2 and 5% CO2; Tasic et al., 2018). The brain was dissected and submerged in ice-cold carbogenated aCSF for 3 min before being transferred to a brain matrix (Braintree Scientific) where four 500 µm sections through M1 were sliced and transferred to a Petri dish containing ice-cold carbogenated aCSF. Regions of M1 (+1.0 to −0.5 mm relative to bregma) containing labeled CSNs were dissected under a fluorescent microscope and transferred to a 5 ml Eppendorf tube with dissociation buffer with papain.
The dissociation buffer contained the following (in mm): 82 sodium sulfate, 30 potassium sulfate, 10 HEPES, 10 glucose, and 5 magnesium chloride (Saunders et al., 2018). Thirty minutes before dissociation, the dissociation buffer was warmed to 34°C, and 5 ml was added to a vial of lyophilized papain (catalog #LK003178, Worthington Biochemical). The solution was diluted 1:2 with additional dissociation buffer before sample dissociation.
The sample in dissociation buffer with papain was incubated at 34°C for 70 min on a shaker at medium speed. After 70 min, the dissociation solution with papain was replaced with a Stop buffer containing 5 ml of dissociation buffer, 5 mg ovomucoid protease inhibitor (catalog #LK003182, Worthington Biochemical) and 10 mg bovine serum albumin (catalog #AB01088, AmericanBio) for 5 min on ice. The Stop solution was then replaced with 800 µl of dissociation buffer and triturated with Pasteur pipettes with decreasing diameters (600, 300, and 150 µm) and placed on ice. To differentiate between intact neurons and nuclei, a cell-permeable dye was added at a concentration of 1:1000 for 5 min (catalog #L34974, Invitrogen). The sample was then centrifuged at 300 × g for 10 min at 4°C and resuspended in 500 µl of dissociation buffer.
Fluorescence activated cell sorting
In total 44 mice were used for the scRNAseq studies, and processing of tissue was completed in seven batches on different days for the following: three FL-traced cohorts (FL1, n = 8, 4 female/4 male; FL2, n = 6, 3 female/3 male; FL3, n = 5, 2 female/3 male), three HL-traced cohorts (HL1, n = 6, 3 female/3 male; HL2, n = 7, 3 female/4 male; HL3, n = 6, 3 female/3 male), and one naive Rbp4 cre:Ai14 (n = 6 rbp4 cre:Ai14, 3 female/3 male). Fluorescent single cells were isolated from rbp4 cre:Ai14 mice and C57BL/6J mice injected with retro-AAV-CAG-GFP. Sorting was performed on the Sony SH800 Cell Sorter with a 130 µm nozzle, ensuring a sheath pressure and sample pressure of <9 psi. To ensure pyramidal neurons were sorted, size reference beads were used to set a minimum diameter of 10 µm [polystyrene particles (PPS)-5 and PPS-6, Spherotech] for all collected cells. To exclude nonintact cells, the cells were sorted for the presence of the cell-permeable dye. A sample of nonfluorescent tissue from adjacent regions of the cortex was used as negative control to set the gate for the presence of fluorescence (GFP or tdTomato). The sample was then sorted into Eppendorf tubes containing 4°C dissociation buffer and immediately sequenced at the Yale Center for Genome Analysis.
In total 44 animals were used (4 rbp4 cre:Ai14; 16 C57BL/6J or Ai14 injected with retro-AAV-CAG-GFP into C6/7; 24 C57BL/6J or Ai14 injected with retro-AAV-CAG-GFP into L4/5). Dissection of L5 of M1 was guided by fluorescence and was consistent across all animals.
10x Single-cell RNAseq
Construction of 10x genomic single-cell/nuclei 3′ RNA-seq libraries, CITE-seq libraries, and sequencing
The single-cell protocol enables short read sequencing to deliver a scalable microfluidic platform for digital gene expression of 500−10,000 individual cells per sample.
Gel beads-in-emulsions generation and barcoding
Single-cell suspension in Real Time (RT) Master Mix was loaded on the Single Cell A Chip and partitioned with a pool of ∼750,000 barcoded gel beads to form nanoliter-scale gel beads-in-emulsions (GEMs). Each gel bead has primers containing (1) an Illumina R1 sequence [read 1 sequencing primer (R1)], (2) a 16 nt 10x barcode, (3) a 12 nt unique molecular identifier (UMI), and (4) a poly-dT primer sequence (30 nt). After dissolution of the Gel Beads in a GEM, the primers were released and mixed with cell lysate and Master Mix. Incubation of the GEMs then produced barcoded, full-length cDNA from polyadenylated mRNA.
Post-GEM-RT cleanup, cDNA amplification, and library construction
Silane magnetic beads were used to remove leftover biochemical reagents and primers from the post-GEM reaction mixture. Full-length, barcoded cDNA was then amplified by PCR to generate sufficient mass for library construction. Enzymatic fragmentation and size selection were used to optimize the cDNA amplicon size before library construction. R1 was added to the molecules during GEM incubation. Primer 5 (P5), P7, a sample index, and R2 (read 2 primer sequence) were added during library construction via End Repair, A-tailing, Adaptor Ligation, and PCR. The final libraries contain the P5 and P7 primers that were used in Illumina bridge amplification.
Sequencing libraries
The Single Cell 3′ protocol produces Illumina-ready sequencing libraries. A Single Cell 3′ library contains standard Illumina paired-end constructs that begin and end with P5 and P7. The Single Cell 3′ 16 bp 10x barcode and 12 bp UMI are encoded in Read 1, whereas Read 2 is used to sequence the cDNA fragment (91 bp). Minimum sequencing depth is 20,000 read pairs per cell. The RNAseq data is deposited in the National Center for Biotechnology Information Gene Expression Omnibus (GEO) and is accessible through the GEO series.
Single-cell RNAseq analysis
Quality control and normalization
Single-cell analyses were performed with Python version 3.8.5 software and the SciPy (version 1.5.2) and NumPy (version 1.19.2) modules. After extracting the single-cell gene expression matrix, we used the scprep (version 1.0.3) package in Python to eliminate non-neurons and normalize the data (https://github.com/KrishnaswamyLab/scprep/releases). Briefly, we filtered out cells that had a library size of <2500 UMIs as well as cells that had mitochondrial gene expression in the top 10th percentile. To exclude doublets, we used two methods. (1) We used the Python package Scrublet and removed all cells that were detected as doublets (Wolock et al., 2019), and (2) we removed cells that expressed both the male-specific genes (Eif2s3y or Ddx3y) and female-specific genes (Xist or Tsix). After a library size normalization, we used the Python package MAGIC (version 3.0.0) to impute missing data (van Dijk et al., 2018). Using the MAGIC output gene expression matrix, we filtered out contaminating cell types (endothelial cells, microglia, astrocytes, and oligodendrocytes) as well as inhibitory interneurons (via negative expression of vglut2), resulting in an approximate 1% reduction from the starting neuronal population. We then square-root transformed the data for further downstream analysis.
Batch correction
To batch correct collected data, we used the ingest function from the Python package Scanpy (version 1.8.2). Once the datasets had been integrated, a graph was built in the integrated space, and MAGIC was used to impute gene expression based on the integrated graph.
To identify supraspinal connectivity of neurons sequenced in this study, we took a similar approach. Briefly, the dataset collected in Yao et al. (2021) was used as the reference dataset. We then used the ingest function in Scanpy and projected the dataset collected in this study into the manifold calculated from Yao et al. (2021). Last, we used the label-mapping output from Scanpy to identify intratelencephalic (IT), extratelencephalic (ET), near-project (NP), and corticothalamic (CT) cells.
Removing putatively mislabeled neurons from downstream analyses
To conduct accurate differential expression analyses, we first needed to identify putatively mislabeled neurons. Because of inherent inefficiencies in viral labeling and fluorescence-activated cell sorting (FACS), we used a computational approach. We approximated the distribution for each of our cell populations (L5 nonspinal, L2/3, and CSNs) using a Gaussian kernel density estimate (KDE) using the Python package SciPy. To identify putative nonspinal nonfluorescent neurons, we calculated the likelihood that each non-CSN (L5 and L2/3) came from the CSN distribution. We removed putative CSNs (KDE > 50th percentile of CSN distribution) from the non-CSN samples. Additionally, we used a similar approach to exclude a portion of L5 nonspinal neurons as they were putative L2/3.
Differential gene expression analysis
Once we identified a pure nonspinally projecting layer 5 neuron population, we used the Python package diffxpy (https://github.com/theislab/diffxpy) to run a Wilcoxon rank-sum test with a Benjamini–Hochberg correction for multiple comparisons to identify differentially expressed (DE) genes between GFP+ CSNs and nonspinal L5 cells. Using the significantly differentially expressed genes (q < 0.05), we used recursive feature elimination using sklearn.feature_selection.RFE (Pedregosa et al., 2011) to identify the top 10 genes that best differentiate between CSNs and non-CSNs. Differentially expressed genes were then used for QIAGEN Ingenuity Pathway Analysis (IPA; Krämer et al., 2014; https://www.qiagen.com/us/products/discovery-and-translational-research/next-generation-sequencing/informatics-and-data/interpretation-content-databases/ingenuity-pathway-analysis) for identification of commonly regulated pathways and for building support machine vectors.
Linear support vector machine classification
To classify cells as either CSNs or non-CSNs, we trained a support vector machine (SVM) classifier with a linear kernel to predict cell origin based on the top 50 enriched genes in CSNs and the top 50 enriched genes in non-CSNs, reserving half of our cells for cross-validation (Pedregosa et al., 2011). The SVM was balanced for class weight as there were an uneven number of CSNs and non-CSNs. Probability outcomes were predicted by setting probability to True during training. SVM performance was evaluated using the area under the receiver operating characteristics curve (AUC ROC). Classifier performance was statistically compared using a randomized permutation test.
The limb-specific classifier was built in a similar way. Briefly, the top 50 enriched forelimb genes and the top 50 enriched hindlimb genes were used to build a linear SVM, with half of the cells reserved for cross-validation. SVM performance was evaluated using AUC ROC. Classifier performances were statistically compared using a randomized permutation test.
Classification of Rbp4+ neurons
Rbp4+ neurons were run through the same QC quality control pipeline outlined above. They were batch corrected with the sequenced cells in this study using the ingest function in Scanpy. Their expression was batch normalized with MAGIC based on the combined graph from the batch correction. Rbp4+ neurons with batch corrected expression values were then run through the CSN-C, and neurons with a probability of being a CSN > 0.75 were classified as CSNs.
CSN visualization
For low-dimensional visualization of all scRNAseq data, we used the Python package PHATE (potential of heat diffusion for affinity-based transition embedding; version 1.0.7) in three dimensions with the default settings (Moon et al., 2019). All PHATE projections are displayed in 2D; examples of 3D animated versions are in Movies 1, 2, 3 (see Figs. 3C, 6D, 7A). Animated versions of all PHATE projections are available on request.
3D PHATE visualization of CSNs, L2/3, and nontraced L5 neurons in M1. 3D PHATE projection of 5495 intact adult neurons from M1 colored by cell identity (traced CSNs, green; nontraced L5, blue; L2/3 neurons, yellow).
3D PHATE visualization of CSNs and deep layer glutamatergic projection neurons sequenced from Yao et al. (2021) in the Yao et al. (2021) manifold. PHATE projection of CSNs and M1 projection neurons from Yao et al. (2021), colored by BICCN taxonomy, ET (orange), IT (green), CT (purple), and NP (red) projection subtype.
3D PHATE visualization of CSNs traced from the cervical and lumbar enlargements colored by BICCN taxonomy in their native manifold. 3D PHATE visualization of CSNs in their native manifold colored by BICCN taxonomy, CSN-ET (orange), CSN-IT (green), CSN-CT (purple), and CSN-NP (red) projection subtype.
MELD analysis of intermediate FL-HL CSN classification
To investigate the distribution of cervical- and lumbar-labeled CSNs, and investigate the possibility of an intermediate population, we used the Python package MELD (version 1.0.0) to identify the degree to which individual CSNs separated into the two groups (cervical or lumbar) defined by labeling, which was expressed as a sample-associated relative likelihood (Burkhardt et al., 2021). This analysis was performed on unimputed square-root-transformed data to avoid the possibility that imputing over the whole dataset could have softened the distinction between the cervical and lumbar groups.
Immunohistochemistry and single-molecule fluorescent in situ hybridization validation
Single-molecule fluorescent in situ hybridization of adult cortical tissue
Ten days after retro-AAV injections, mice were killed with isoflurane and transcardially perfused with 0.9% NaCl with 10 units/ml heparin (1000 unit/ml; Covetrus) followed by 4% paraformaldehyde (PFA) in PBS. Brains and spinal cords were dissected and postfixed in 4% PFA overnight at 4°C. The next day, 15-µm-thick sections were cut using a vibratome (VT1000 S, Leica), mounted on Superfrost Plus slides and processed for single-molecule fluorescent in situ hybridization (smFISH) according to the ACD RNAScope fluorescent protocol. Briefly, sections were either immediately stored at −80°C until they were ready for processing, or sections were postfixed in 4% PFA and stained with a single probe, Wnt7b (catalog #401131, Advanced Cell Diagnostics), Cacng7 (catalog #556171-C2, Advanced Cell Diagnostics), and Slc16a2 (catalog #545291, Advanced Cell Diagnostics). Each experiment included a positive and negative control probe to ensure validity. Before coverslipping, immunofluorescence using primary antibodies directed against GFP (1:2000; catalog #ab13970, Abcam), and mCherry (1:2000; catalog #ab167453, Abcam) and detected with secondary antibodies Alexa Fluor 488 and 568, (1:500, 1:500; catalog #A11011, Life Technologies) was used to visualize traced CSNs.
smFISH of adult cortical tissue analysis
All analyses were performed by an experimenter blinded to probe identity. For quantification of smFISH results, images were taken at 40× on a Leica SP8 confocal microscope. Three images were taken per hemisphere, three sections were taken per animal, and three animals were used for the analysis of each probe. For analysis, each image was separated into individual channels (probe, forelimb and hindlimb CSN, and DAPI). The images of the traced neurons were then exported to Adobe Photoshop, where the quick selection tool was used to select cell body regions of interest (ROIs) and create masks. This was done manually as CSN axons can be very thick, and automatic ROI selection often cannot distinguish between CSN axons and cell bodies. Images of the probe and the ROI masks are then uploaded to CellProfiler software, and puncta per cell are automatically calculated according to the CellProfiler manual (McQuin et al., 2018).
Immunohistochemistry protocol
Ten days after retro-AAV injections, mice were killed with isoflurane and transcardially perfused with 0.9% NaCl with 10 units/ml heparin (1000 unit/ml; Covetrus) followed by 4% paraformaldehyde in PBS. Brains and spinal cords were dissected and postfixed in 4% PFA overnight at 4°C, and 40-µm-thick sections of brain were then sectioned on a vibratome (Leica Microsystems). The sections were then processed for mCherry with tyramide signal amplification (catalog #NEL700A001KT, PerkinElmer). Immunofluorescence used primary antibodies directed against mCherry (1:500,000; catalog #ab167453, Abcam), Synaptophysin (1:10,000; catalog #MAB329, EMD Millipore), and NeuN (1:200; catalog #24307, Cell Signaling Technology) and detected with secondary antibodies Alexa Fluor 488, 568, and 647 (1:500; catalog #A11039, #A11011, #A21235, Life Technologies).
BrainJ processing
Ten days after AAV injections, mice were killed with isoflurane and transcardially perfused with 0.9% NaCl with 10 units/ml heparin (1000 unit/ml; Covetrus), followed by 4% PFA in PBS. Brains were dissected and postfixed in 4% PFA overnight at 4°C. The next day, 75-µm-thick sections were cut using a vibratome (Leica VT1000 S). Sections underwent an immunohistochemistry protocol to enhance tdTomato expression for microscopy. Briefly, sections received 3× washes of PBS-T, followed by a 10 min incubation in 0.3% H2O2 at room temperature, three more washes of PBS-T, and blocking in 10% normal donkey serum for 1 h at room temperature. Next, sections were incubated overnight in Rb-anti-mCherry (1:1000). The next day, sections were washed 3× in PBS-T and then incubated for 2 h in 546 G-anti-Rb (1:1000). After 3× final PBS-T washes, every other section was mounted on Superfrost Plus slides with DAPI mounting media (1 mg/ml; 1:10,000). Sections were imaged at 4× using an Olympus VS200 Slide Scanner for projection analysis. BrainJ software was used, as previously described (Botta et al., 2020), to perform analysis of supraspinal forelimb and hindlimb ipsilateral and contralateral projections. Briefly, sections were ordered and registered to the Allen Mouse Brain Common Coordinate Framework (CCF) (http://atlas.brain-map.org/) template using Elastix software, version 5.0.1. Next, ilastix, version 1.3.3 was used to identify brainwide soma and projections in the ∼80 mounted sections using pixel-based machine learning. BrainJ output absolute density measurements for 1429 subregions of the brain, as defined by the Kim enhanced and unified mouse brain atlas (Chon et al., 2019). We then classified each subregion as ET, IT, NP, or CT using the Allen Institute for Brain Science framework. A mixed-design ANOVA was performed in R software for each projection subtype to analyze differences between forelimb–hindlimb and ipsilateral–contralateral projections. An alpha level of 0.05 was set a priori, and Bonferroni-adjusted post hoc tests were performed when main effects were found.
Data availability
All analyses were performed using in-house-developed code and implemented in Python. The CSN and limb-specific CSN gene lists along with example classifiers and appropriate sample data and code can be found on GitHub at https://github.com/cafferty-lab/CSN_classification.
Results
A novel pipeline for scRNAseq of retrogradely labeled CSNs
CSNs extend axons from L5b of M1 to innervate every spinal segment, synapsing principally on spinal interneurons to control FL and HL motor output. Spinal retrograde tracing from C6/7 and L4/5 show that CSNs that innervate the FL are spatially and anatomically distinct from those that innervate the HL, mirroring functional dexterity differences and hinting that FL and HL CSNs are also molecularly heterogeneous (Fig. 1A–D'). To explore the molecular heterogeneity among CSN populations, we completed scRNAseq on retrogradely labeled FL and HL CSNs and compared gene expression between these anatomically discrete populations with layer 2/3 and nonspinal projecting L5 neurons in M1 of adult mice. As relatively large pyramidal neurons, CSNs have a lower nucleus to cytoplasm ratio, hence a large portion of transcripts are extranuclear (Bakken et al., 2018). Therefore, it was critical to maintain plasma membrane integrity to achieve maximal quality for whole-cell sequencing. To that end, we created a pipeline for labeling, extracting, dissociating, and sequencing adult neurons with long-distance axonal projections (Golan and Cafferty, 2021). To independently label FL and HL CSNs, we microinfused retro-AAV-CAG-GFP into the gray matter of either C6/7 or L4/5 spinal cord in adult mice (Fig. 1E). After 10 d we macrodissected layers 2/3 and 5 of M1 (+1.0 to −0.5 mm relative to bregma), enzymatically and physically dissociated the tissue to create a single-cell suspension. Dissociated CSNs were then incubated in a cytoplasmic dye to aid whole-cell detection for subsequent FACS (Golan and Cafferty, 2021). Intact CSNs were then collected for scRNAseq using the 10× Chromium system. As not every CSN was fluorescently labeled because of inherent inefficiencies with viral labeling, KDE analysis was used to identify putative CSNs from the nonfluorescent layer 5 samples based on gene expression. Control layer 5 cells with a >50% probability of originating from the CSN distribution were classified as putative CSNs and excluded from downstream analyses (Fig. 1F). To overcome inexact microdissection of layer 2/3 neurons from layer 5 neurons, KDE analysis was used to exclude putative layer 5 neurons from the layer 2/3 sample for downstream analysis (Fig. 1G).
Single-cell RNAseq pipeline for adult CSNs. A, Schematic of spinal tracing via injecting retro-AAV-CAG-GFP at spinal level C6/7 and retro-AAV-CAG-tdTomato at spinal level L4/5. B–D', Photomicrographs from M1 (∼0.0 mm AP bregma, medial cortex) showing FL-CSNs (B), HL-CSNs (C), and an overlay with both single and dual-projecting CSNs (D, D') 14 d after spinal infusion. Schematic overview of the dissection-sequencing procedure. Ei–iv, Mice received injections of retro-AAV-CAG-GFP into either C6/7 or L4/5 (E). After a 2 week incubation, three 500 µm sections through M1 were collected (Ei), layers 2/5 and 5 were macrodissected (Eii) and dissociated, and traced CSNs (green), L2/3 cells (yellow), and L5 nonfluorescent cells (blue) were then sorted using FACS (Eiii) and sequenced using 10x Chromium (Eiv). To exclude putative CSNs from the unlabeled layer 5 population of sequenced neurons, we used KDE analysis to quantify the transcriptional overlap between populations in low-dimensional space. F, Right, untraced layer 5 neurons (blue) with >50% probability (blue, white outline) of originating from the CSN distribution (green) were excluded. To exclude putative layer 5 neurons from the layer 2/3 population of sequenced neurons, we used KDE to quantify the transcriptional overlap between these populations in low-dimensional space. G, Right, Layer 2/3 neurons (orange) with >50% probability (orange, white outline) of originating from the layer 5 distribution (blue) were excluded. Scale bars: B–D, 500 µm; D', 100 µm.
Sequenced neurons had a mean library size of 19,432 UMIs (median, 18,474 UMIs; Fig. 2A) and 5357 genes per cell (median, 5537). Previous studies have detailed molecular differences between and among excitatory and inhibitory neurons spanning the entire cortical depth (Arlotta et al., 2005; Lein et al., 2007; Morris et al., 2010; Molyneaux et al., 2015; Tasic et al., 2016; Saunders et al., 2018; Tasic et al., 2018). Here, we focused specifically on molecular differences between layer 2/3 neurons and L5 neurons in M1, as well as differences among L5 neurons, and comprehensively characterize transcriptional specificity between CSNs. To achieve this, we compared gene expression between retrogradely traced CSNs (n = 1776 cells), nonfluorescent nonspinal L5 neurons (n = 3000 cells), and layer 2/3 neurons (n = 719) in M1, sequencing a total of 5495 adult neurons.
Gene expression confirms sequenced cell identity. A, Sequenced CSNs, L2/3 neurons, and nontraced L5 cells have overlapping distributions of library size; mean library size of 19,432 UMIs; median, 18,474 UMIs. B, Differential gene expression analysis revealed 32 genes were enriched in neurons from male mice and 213 genes enriched in neurons from female mice (Wilcoxon rank sum test with a Benjamini–Hochberg correction, p adjusted < 0.05). C, IPA was used to gain biological insight between the sexes; only three pathways were significantly enriched—the synaptogenesis, protein ubiquitination, and eNOS signaling (right-tailed Fisher's exact test with a Benjamini–Hochberg correction, p < 0.05). D–H, Heat maps representing normalized imputed expression of sex-specific genes (D), intact cell (IC) genes (E), cell-type genes (F), neurotransmitter (NT) genes (G), and cortical layer genes (H), confirming that sequenced cells were sex balanced, intact, neurons, glutamatergic, and from expected layers, respectively.
CSNs are unique among neurons in M1 and are enriched in Wnt7b expression. A, Differential gene expression analysis between traced CSNs and a combination of layer 2/3 and L5 nonfluorescent cells shows that 197 genes were upregulated and 639 downregulated in CSNs, 217 genes were upregulated and 634 downregulated in CSNs compared with layer 5 nontraced neurons, and 467 genes were upregulated and 1533 downregulated in CSNs compared with layer 2/3 neurons (Wilcoxon rank sum test with a Benjamini–Hochberg correction, p adjusted < 0.05). B, QIAGEN IPA was used to identify enriched pathways in each of these analyses, (right-tailed Fisher's exact test with a Benjamini–Hochberg correction, p adjusted < 0.05). C, D, PHATE projection (C, Movie 1) and kernel density estimate plots with hierarchical clustering and a similarity index (D) colored by cell identity illustrate differences between CSNs, L5 nontraced, and L2/3 neurons. E, Wnt7b is significantly enriched in CSNs compared with layer 5 nontraced and L2/3 neurons (Wilcoxon rank sum test with a Benjamini–Hochberg correction; ****p adjusted < 0.001). F, Schematic showing labeling of cervical CSNs with retro-AAV-CAG-GFP (rAAV-GFP, cyan) and lumbar CSNs with retro-AAV-CAG-tdTomato (rAAV-tdT, magenta). G, H, smFISH of retrogradely labeled CSNs shows an enrichment of Wnt7b in FL-CSNs and HL-CSNs compared with nonfluorescent DAPI+ cells in layer 5b (data shown are average number of puncta per cell ± SEM, one-way ANOVA with a Bonferroni correction, ****p < 0.0001). Scale bars: 10 µm. Data tables from the scRNAseq containing the results of the DE analysis comparing CSNs to nontraced L5 and L2/3 neurons are provided in Extended Data Figure 3-1.
Figure 3-1
DE genes for CSNs versus non-CSNs in M1. Data tables from the scRNAseq containing the results of the DE analysis comparing CSNs to nontraced L5 and L2/3 neurons (Wilcoxon rank sum test with a Benjamini–Hochberg correction, p adjusted < 0.05). Download Figure 3-1, XLSX file.
Differential gene expression analysis of neurons from male and female mice revealed 32 DE genes in male mice and 213 DE genes in female mice (Fig. 2B). QIAGEN IPA revealed that the synaptogenesis, protein ubiquitination, and eNOS (endothelial nitric oxide synthase) signaling pathways were differentially enriched between male and female mice (Fig. 2C). To ensure there was no sex bias in this study, we assessed the expression of female-specific genes Xist and Tsix as well as the male-specific genes Ddx3y and Eif3s3y (Fig. 2D), showing female and male representation in all assessed populations. Next, expression of Ywhaz, a gene whose mRNA is enriched in the cytoplasm (Yao et al., 2021), and Neurod2, a nuclear gene, confirms that all sequenced cells contained both a nucleus and cytoplasm (Fig. 2E). To broadly confirm the identity of sequenced cells, we assessed their expression of previously identified cell-type-specific genes. All sequenced cells expressed the neuron-specific gene Snap25 and showed a dearth in expression for other cell-type-specific genes including the endothelial-specific gene Cldn5, the astrocyte-specific gene Gfap, the microglial-specific gene Cx3cr1, and the oligodendrocyte-specific gene Olig2 (Fig. 2F). Neurotransmitter gene expression confirmed that all sequenced neurons were enriched in the glutamate transporter Slc17a7 compared with the GABAergic gene Gad1, the cholinergic gene Chat, or the serotonergic gene Sert (Fig. 2G). Layer 2/3 neurons showed an elevated expression for the previously identified layer 2/3-enriched gene Cux2 (Tasic et al., 2016). CSNs and nonspinal L5 neurons showed a moderate enrichment for established L5 genes Ctip2, Fezf2, and Crym, relative to other cortical layer genes Reln (L1) and Ctgf (L6; Fig. 2H), confirming L5 identity. Although CSNs expressed previously identified L5 genes, none were specific to any L5 population, suggesting that more targeted analyses are necessary to identify unique molecular characteristics of CSNs.
Identifying CSN-enriched genes in M1
To identify CSN-enriched genes within M1, we performed the following DE analyses: (1) traced CSNs versus non-CSNs (nonspinal L5 neurons and layer 2/3 neurons), revealing 197 CSN upregulated genes and 639 downregulated genes; (2) traced CSNs versus nonspinal L5 neurons, revealing 215 CSN upregulated genes and 634 downregulated genes; and (3) traced CSNs versus layer 2/3 neurons, revealing 467 CSN upregulated genes and 1533 downregulated genes (Fig. 3A; Extended Data Fig. 3-1 contains the entire DE lists). To explore the biological import of the transcriptional differences among CSNs and other glutamatergic neurons in M1, we completed QIAGEN IPA (Krämer et al., 2014). Consistently among the top enriched pathways across all comparisons were those associated with axon guidance, signal transduction, and synaptic transmission (Fig. 3B). These pathways are consistent with previous work showing that axon guidance cues including Epha4, Efnb3, Ryk, and Wnt3a have previously been shown to be critical in wiring the CST (Canty and Murphy, 2008).
To explore differences among these populations, we used PHATE for dimensionality reduction and visualization, as unlike the visualization techniques PCA, t-SNE, and UMAP, PHATE maintains both local and global structure and was designed for scRNAseq visualization (van Dijk et al., 2018; Moon et al., 2019). A PHATE projection colored by cell origin (Fig. 3C; Movie 1) illustrates the transcriptional overlap between these populations. To quantify, KDE analysis of each population was plotted, and a similarity index was calculated based on the overlapping regions of the KDE (Fig. 3D). These analyses demonstrate that CSNs are most similar to nonspinal L5 neurons (similarity index = 0.5), whereas they are least similar to layer 2/3 neurons (similarity index = 0.23). Nonspinal L5 neurons show an intermediate level of similarity with layer 2/3 neurons in comparison with CSNs (similarity index = 0.41).
Based on enrichment of guidance cues in CSNs, we selected Wnt7b as a primary candidate to validate our scRNAseq data as it is significantly enriched in all traced CSNs (Fig. 3E). To confirm this finding in vivo, we infused retro-AAV-CAG-GFP into the gray matter of C6/7 spinal cord and retro-AAV-CAG-tdTomato into the gray matter of L4/5 spinal cord in adult wild-type mice, differentially labeling FL and HL CSNs and completing smFISH (Fig. 3F). Photomicrographs of L5b from M1 prepared 2 weeks after spinal infusion show GFP+ FL CSNs and tdTomato+ HL CSNs are significantly enriched in Wnt7b relative to DAPI+ control cells (Fig. 3G,H).
Classifier analysis to predict CSN identity
Based on the high similarity index between CSNs and nonspinal L5 neurons, we reasoned it would be unlikely that any single gene could uniquely identify CSNs from other L5 neurons. To that end, we first explored 12 genes that have previously been used to describe L5 pyramidal neurons (Fig. 4A; Arlotta et al., 2005; Fink et al., 2015; Tasic et al., 2016; Economo et al., 2018; Tasic et al., 2018; Poplawski et al., 2020). Of these genes, many are enriched in CSNs compared with other populations including Ctip2, Cd55, Cdh13, Crym, and Fezf2. However, all are also expressed at moderate levels in nontraced L5 neurons and are only detected in a subset of CSNs. Accordingly, PHATE projections for Ctip2 and Rbp4 show an incomplete and nonspecific expression in CSNs, and although Wnt7b shows greater specificity, coverage remains incomplete (Fig. 4B). To improve on these genes as a means for CSN identification, we sought to distinguish a CSN transcriptional signature. To that end, we trained a linear classifier to predict cell identity based on expression of 100 genes that are most differentially expressed between CSNs and non-CSNs (CSN-C). For comparison, we trained additional CSN classifiers based on the top 10 genes differentially expressed and a random 10 genes from the top 100 that were differentially expressed (Fig. 4C). We compared the performance of these classifiers to one that was trained on sole expression of all 12 previously characterized L5 genes, termed L5 classifier (L5-C), and another classifier trained on a random set of 10 genes. The performance of each classifier was evaluated using the ROC AUC, and the classifiers were compared using a permutation test. This revealed that CSN-C (AUC = 0.96), top 10 genes (AUC = 0.95), random 10 CSN genes (AUC = 0.91) significantly outperform L5-C (AUC = 0.87; Fig. 4C). To visually compare the differences between these classifiers, we plotted CSN-C CSNs and L5-C CSNs in PHATE. CSN-C correctly classified 92% of CSNs, whereas L5-C classified 78% of CSNs (Fig. 4D).
Using linear SVM classifiers to identify CSNs among cells in M1. A, Dot plot of previously characterized L5-enriched genes illustrates enrichment and expression of each gene in CSNs, nontraced L5, and L2/3 neurons. The size of the dot denotes the proportion of each population that expresses a given gene, and the color denotes the normalized mean expression within that subtype. B, PHATE projections colored by normalized gene expression of previously characterized L5 genes including Ctip2 and Rbp4, along with newly identified wnt7b, illustrate expression in CSNs, nontraced L5 neurons, and L2/3 neurons. Top right, Boxed PHATE projection illustrates neuronal cell types in their transcriptional space (CSNs, green; nontraced L5 neurons, blue; L2/3 neurons, orange). C, An ROC curve showing the performance of an SVM) trained to classify CSNs based on L5-C genes; the top 100, top 10, or random 10 novel genes enriched in CSN-C, or 10 random genes. AUC is calculated for each classifier, showing that CSN-C outperforms L5-C, and top 10 and random 10 CSN genes outperform the L5-C (randomized permutation test, p < 0.001). D, PHATE projections showing CSN classification performance with either the CSN-C (top left, forest green) or L5-C (bottom left, purple), Top right, CSNs missed by L5-C (black), and a quantification of the proportion of CSNs that were classified correctly. CSN-C correctly classified 92% of CSNs. L5-C correctly classified 76% of CSNs. E, Rbp4 cre:Ai14 mice received infusions of retro-AAV-CAG-GFP into either the cervical, C6/7, or the lumbar, L4/5, spinal cord. F, F', Low-power (F) and high-power (F') photomicrographs show GFP (cyan) expression in a section through M1 from an Rbp4 cre:Ai14 (magenta) mouse. Cyan cells represent traced CSNs. Magenta cells are Rbp4+ neurons that do not project to the spinal cord, white cells are GFP+/Rbp4+ neurons (F'). G, Quantification of the proportion of Rbp4+ neurons that express GFP in coronal sections through M1 (1.7 to −1.2, relative to Bregma; mean, 7.7%). CSN-C prediction of CSN identity of sequenced Rbp4+ neurons. H, CSN-C predicts 7% of Rbp4 neurons are CSNs with a 75% probability cutoff. Scale bars, F, 500 µm; F', 100 µm.
Classifier analysis to predict CSN identity of Rbp4+ neurons
To validate the use of the CSN-C as a tool for predicting CSN identity, we sequenced neurons from Rbp4 cre:Ai14 mice, a line previously used to comprehensively target and label L5 pyramidal neurons (Alcamo et al., 2008; Tasic et al., 2018). To quantify the proportion of Rbp4+ neurons that project to the spinal cord (Rbp4+ CSNs), we delivered retro-AAV-CAG-GFP into the cervical and lumbar cord of Rbp4 cre:Ai14 mice (Fig. 4E). Photomicrographs show that a minority of Rbp4+ neurons project to the spinal cord (Figs. 4F,F'). Quantification of the proportion of Rbp4+ neurons that expressed GFP, and therefore have terminals in either the cervical or lumbar spinal cord, revealed that on average, 7.7% of Rbp4+ neurons are also CSNs (Fig. 4G). Classification of sequenced Rbp4+ neurons using CSN-C predicted that 7% of Rbp4+ neurons were putative CSNs using a 75% probability cutoff (Fig. 4H). These findings confirm the validity of our CSN gene list to predict CSN identity among L5 neurons in M1.
Molecular identification of limb-specific CSN subtypes
We next wanted to explore potential transcriptional heterogeneity between anatomically discrete FL- and HL-projecting CSNs. PHATE visualization of traced neurons colored by anatomic origin illustrates there is a transcriptional separation between FL and HL traced cells (Figs. 5A). To quantify, a KDE of each population was plotted, and a similarity index was calculated based on the overlapping regions of the KDE. FL and HL CSNs are somewhat similar to each other (similarity index = 0.74). Because of the elevated similarity index relative to the more definitive separation between CSNs versus L5 controls and CSNs versus layer 2/3 (Fig. 3D), we reasoned that a population of CSNs may have a projection to both the cervical and lumbar spinal cord. Indeed, a small percentage of traced CSNs were double positive for GFP and tdTomato when tracers were injected into the cervical and hindlimb respectively (Fig. 1A–D'). To explore whether these cells represented a true population of FL/HL CSNs or whether en passant labeling of HL axons occurred when infusing tracer into the cervical cord, we used MELD to identify CSNs whose transcriptomic characteristics were intermediate between the FL and HL populations as previously identified by cervical and lumbar tracing (Fig. 5B,C). We found 188 such CSNs, 60 of which had been previously identified as FL and 128 as HL (Fig. 5D). These data suggest that as much as 10% of CSNs may have a projection into both the cervical and lumbar spinal cord.
A, Forelimb and hindlimb CSN specialization. PHATE projection colored by cell identity of retrogradely labeled GFP+ cervical CSNs (cyan) and tdTomato+ lumbar CSNs (magenta). B, Histogram of MELD prediction confidence of lumbar identity in all CSNs, showing predicted intermediate cells in gray. C, PHATE projection showing MELD prediction confidence of lumbar identity. D, Comparison of counts for tracing-identified lumbar and cervical CSN populations (top) with MELD-predicted counts for lumbar, cervical, and intermediate CSN populations (bottom). E, DE analysis identified 492 upregulated FL genes and 788 upregulated HL genes (Wilcoxon rank sum test, Benjamini–Hochberg correction, p adjusted < 0.05). F, Violin plots of the top 10 most differentially expressed genes in FL CSNs (cyan) and HL CSNs (magenta; Wilcoxon rank sum test, Benjamini–Hochberg correction, p adjusted < 0.05). G, PHATE projection colored by Cacng7 expression illustrates FL enrichment. H, PHATE projection colored by Slc16a2 expression illustrates HL enrichment. I, J, smFISH of retrogradely labeled GFP+ FL-CSNs and tdTomato+ HL-CSNs shows an enrichment of Cacng7 in FL-CSNs (data shown are average number of puncta per traced CSN ± SEM; unpaired t test, ***p < 0.0001). K, L, smFISH of retrogradely labeled GFP+ FL-CSNs and tdTomato+ HL-CSNs shows an enrichment of Slc16a2 in HL-CSNs (data shown are average number of puncta per traced CSN ± SEM; unpaired t test, ***p < 0.001). M, Dot plot of previously characterized L5-enriched genes illustrates enrichment and expression of each gene in FL and HL CSNs. The size of the dot denotes the proportion of each population that expresses a given gene, and the color denotes the normalized mean expression within that subtype. These L5 genes are not uniquely nor differentially enriched in FL or HL CSNs. N, An ROC curve showing the performance of a linear SVM trained to classify limb-specific CSNs based on previously characterized L5 genes (L5-C) or our novel limb-specific gene list (LS-CSN-C). AUC is calculated for each classifier, showing that LS-CSN-C outperforms C-L5. Scale bars; K, bottom left 30 µm; K, bottom right, 15 µm. Data tables from the scRNAseq containing the results of the DE analysis comparing FL versus HL CSNs are provided in Extended Data Figure 5-1.
Figure 5-1
DE genes for FL versus HL CSNs. Data tables from the scRNAseq containing the results of the DE analysis comparing FL versus HL CSNs (Wilcoxon rank sum test with a Benjamini–Hochberg correction, p adjusted < 0.05). Download Figure 5-1, XLSX file.
Transcriptionally distinct CSN subtypes can be defined by supraspinal terminal fields. A, Schematic based on the BICCN M1 taxonomy showing that projection neurons can be classified into four categories based on their primary terminal field, including ET (areas shaded in orange in the schematic), IT (areas shaded green), CT (areas shaded purple), and NP (areas shaded red). B, PHATE projection showing 10x V3-sequenced M1 neurons from Yao et al. (2021) colored by ET (orange), IT (green), CT (purple), and NP (red) projection subtype. C, PHATE projection of CSNs collected in this study projected onto the same PHATE space as in B from Yao et al. (2021) colored by projection subtype show that CSNs are represented in all four projection classes. D, PHATE projections of CSNs and M1 projection neurons from Yao et al. (2021; Movie 2). E, Violin plot showing gene expression in each projection subtype. F, PHATE projections showing enrichment of established ET projection gene Fam84b and novel ET enriched gene Gprc5, enrichment of established IT gene Slc30a3 and novel IT enriched gene Frzb, enrichment of established CT gene Foxp2 and novel CT enriched gene Syt6, and enrichment in established NP gene Sla2 and novel NP enriched gene Lypd1. G, DE analysis identified 791 upregulated and 181 downregulated genes in CSN-ET cells versus CSN-IT, CSN-CT, and CSN-NP cells, listed with the top six enriched pathways from QIAGEN IPA (Wilcoxon rank sum test, Benjamini–Hochberg correction, p adjusted < 0.05). H, DE analysis identified 275 upregulated and 263 downregulated genes in CSN-IT cells versus CSN-ET, CSN-CT, and CSN-NP cells, listed with the top six enriched pathways from QIAGEN IPA (Wilcoxon rank sum test, Benjamini–Hochberg correction, p adjusted < 0.05). I, DE analysis identified 215 upregulated and 441 downregulated genes in CSN-CT cells versus CSN-ET, CSN-IT, and CSN-NP cells, listed with the top six enriched pathways from QIAGEN IPA (Wilcoxon rank sum test, Benjamini–Hochberg correction, p adjusted < 0.05). J, DE analysis identified 675 upregulated and 1060 downregulated genes in CSN-NT cells versus CSN-ET, CSN-IT, and CSN-CT cells, listed with the top six enriched pathways from QIAGEN IPA (Wilcoxon rank sum test, Benjamini–Hochberg correction, p adjusted < 0.05). To validate CSN innervation of ET, IT, CT, and NP structures, retro-AAV-CAG-tdTomato was injected in the cervical and lumbar cord. K–N', Terminals in supraspinal brain regions were examined including the brainstem (ET, orange, K), the striatum (IT, green, L), motor thalamus (CT, purple, M), and primary motor cortex (NP, red, N). Photomicrographs of CSN terminals in the brainstem (ET, K'), striatum (IT, L'), motor thalamus (CT, M'), and primary motor cortex (NP, N') stained with antibodies against synaptophysin (green) and NeuN (white) show presumptive synaptic connections within these supraspinal terminal fields. Scale bars; K'-N', 100 µm; insets in each, 30 µm. Data tables from the scRNAseq containing the results of the DE analysis comparing each CSN supraspinal subtype versus all the other subtypes are provided in Extended Data Figure 6-1.
Figure 6-1
DE genes among CSN supraspinal subtypes. Data tables from the scRNAseq containing the results of the DE analysis for comparing each CSN supraspinal subtype versus all the other subtypes (Wilcoxon rank sum test with a Benjamini–Hochberg correction, p adjusted < 0.05). Download Figure 6-1, XLSX file.
DE analysis revealed that there were 492 FL-enriched genes and 788 HL-enriched genes (Fig. 5E; Exended Data Fig. 5-1, entire DE lists). Identifying a set of enriched FL and HL CSN genes is central to understanding the potential independent homeostatic biology and injury-induced changes within these discrete populations. To prioritize candidate genes, recursive feature elimination was used to determine the top 10 FL- and HL-specific genes that most robustly distinguished these two populations. The top genes that define the FL population include calcium-related genes Cacng7 and Efcab12, as well as transporters Slc27a1 and Slc4a2 (Fig. 5F). The top genes that define the HL population include ion transporters such as Slc16a2, transcription-related genes such as Zfp281, and receptors such as Ptger1 (Fig. 5F).
PHATE projections show Cacng7 expression is significantly enriched in FL traced neurons (Fig. 5G) and Slc16a2 expression is significantly enriched in HL traced neurons (Fig. 5H). To validate these results in vivo, mice received injections of retro-AAV-CAG-GFP into spinal gray matter at C6/7 and retro-AAV-CAG-tdTomato into spinal gray matter at L4/5. Visualization of Cacng7 using smFISH confirmed differential expression in GFP+ FL–projecting CSNs compared with tdTomato+ HL–projecting CSNs (Fig. 5I,J). A similar analysis was then repeated with the HL enriched gene Slc16a2. In vivo smFISH confirmed HL-enrichment of Slc16a2 in tdTomato+ HL–projecting CSNs relative to GFP+ FL–projecting neurons (Fig. 5K,L).
Classifier analysis to predict FL and HL identity
Although DE analysis revealed FL- and HL-enriched genes, because of the high similarity index between FL and HL CSN populations, none of these genes can independently predict FL or HL identity. To explore, we first examined limb specificity in the previously identified L5-enriched genes (Fig. 5M). Crym is modestly enriched in FL-CSNs (Fink et al., 2015), and Crim1 was modestly enriched in HL-CSNs (Sahni et al., 2021a,b). As the majority of these genes showed no limb specificity, we trained a linear classifier to predict limb specificity based on expression of 100 genes that are most differentially expressed between FL and HL CSNs, termed limb-specific CSN-C (LS-CSN-C). We compared the performance of this classifier with L5-C. We evaluated the performance of each classifier with ROC AUC, and the two classifiers were compared using a permutation test. This revealed that LS-CSN-C (AUC = 0.94) significantly outperforms L5-C (AUC = 0.76; Fig. 5N).
CSN subtypes defined by supraspinal connectivity
Recent studies from the Brain Initiative Cell Consensus Network (BICCN; Callaway et al., 2021) performing scRNAseq on M1 and anterolateral motor cortex (ALM) have classified glutamatergic projection neurons into the following four broad subtypes: IT neurons that have many projection targets including the striatum and cerebral cortex (Lin et al., 2018; Zhang et al., 2021); ET neurons that project from M1 to the brainstem and medulla (Economo et al., 2018; Tasic et al., 2018; Yao et al., 2021); NP neurons that have projections in local areas including M1, M2, and S1 (Tasic et al., 2018); and CT) neurons (Yao et al., 2021; Fig. 6A). The 10×-sequenced deep layer glutamatergic neurons collected in Yao et al. (2021) separate into four discrete populations representing the IT, ET, NP, and CT subtypes when examined in PHATE (Fig. 6B). CSNs collected in this study, when projected into the same PHATE graph, populate all four subtypes, and therefore we classified them as CSN-IT, CSN-ET, CSN-NP, and CSN-CT (Fig. 6C,D; Movie 2). Differential gene expression analysis among these four putative subdivisions of CSNs revealed an enrichment for previously identified projection subtype genes, Fam84b for CSN-ET cells, Slc30a3 in CSN-IT cells, Foxp2 in CSN-CT, and Sla2 in CSN-NP cells (Tasic et al., 2018; Yao et al., 2021; Fig. 6E), and an additional set of subcategory-specific genes (Fig. 6E,F). We completed QIAGEN IPA on genes significantly differentially expressed among these CSN subcategories (Figs. 6G–J). When comparing the CSN-ET population with the other three categories, we found there were 972 DE genes showing an enrichment in pathways associated with synaptic signaling and neurovascular coupling (Fig. 6G). Comparing CSN-IT cells with the other three categories, we found there were 539 DE genes showing an enrichment in pathways associated with axon guidance and CREB signaling (Fig. 6H). Comparing CSN-CT cells with the other categories, we found there were 656 DE genes showing an enrichment in axon guidance and RhoGDI signaling (Fig. 6I). Comparing CSN-NP cells with the other categories we found there were 1735 DE genes showing an enrichment in axon guidance and cAMP signaling and GPCR signaling (Fig. 6J). The entire DE lists are available in Extended Data Figure 6-1.
As CSNs are represented in all four projection subtype categories, we infer that transcriptional subdivisions among CSNs are therefore defined by the spinal termination of their primary axon in the spinal cord and by the location of their supraspinal terminal fields (Fig. 6K–N).
To explore CSN supraspinal connectivity in vivo, we injected retro-AAV-CAG-tdTomato unilaterally into the gray matter of C6/7 and L4/5 spinal cord and surveyed for the presence of tdTomato+ve CSN terminals in the BICCN identified supraspinal structures (Fig. 6A,K–N). We observed td-Tomato+ve CSN terminals in the brainstem (Fig. 6K,K'), the caudoputamen of the striatum (Fig. 6L,L'), the motor thalamus (Fig. 6M,M'), and regions near M1 including S1 and M2 (Fig. 6N,N'), confirming that CSNs project to both the spinal cord and to these supraspinal structures. To further explore whether these collaterals represented putative active synapses, we completed immunohistochemistry, showing that traced CSN terminals express the presynaptic protein synaptophysin and colocalize with NeuN-positive neurons.
Limb-specific specialization among supraspinal CSN subtypes
We found that traced CSNs maintained a supraspinal cluster separation when plotted in their native PHATE space (Fig. 7A; Movie 3), and FL and HL CSNs were represented in each supraspinal subcategory (Fig. 7B). However, the relative distribution of CSN subcategory distribution was different between FL and HL. For FL CSNs 19.6% were ET, 36.5% IT, 40% CT, and 3.8% NP. HL CSNs were 37.6% ET, 31.5% IT, 26.8% CT, and 4% NP (Fig. 7C,D). These data suggest that FL and HL CSNs may differentially innervate supraspinal structures. To explore we used intersectional viral tracing to independently label all CSN terminals whose axons innervate either cervical or lumbar spinal cord. We injected retro-AAV-CAG-Cre into either C5/6 or L4/5 spinal cord and AAV-CAG-FLEX-tdTomato into FL M1 or HL M1, respectively. Using the Java plug-in BrainJ (Botta et al., 2020), we analyzed the density of tdTomato+ve CST projections across the entire CST. We found that the projection density was higher on the ipsilateral side compared with the contralateral side (mixed-model ANOVA, p < 0.001). We also found that the FL-CST had a higher density of projections compared with the HL-CST (mixed-model ANOVA p = 0.008), consistent with prior reports that FL-CSNs are more numerous than HL-CSNs (Wang et al., 2018). Next, we assessed FL and HL-CSN supraspinal projection density. We found an enrichment of FL-CSN-ET projections compared with HL-CSN-ET projections (Fig. 7E,F; mixed-model ANOVA p = 0.003), and FL-CSN-NP projections compared with HL-CSN-NP projections (Fig. 7Q,R; mixed-model ANOVA p < 0.001). In contrast, we found no significant differences between FL and HL for IT (Fig. 7I,J; p = 0.38) and CT (Fig. 7M,N; p = 0.085) projections. Overall, these data demonstrate the supraspinal anatomic diversity of FL-CSNs and HL-CSNs and suggest that the transcriptional differences among these population support unique circuit wiring.
Forelimb- and hindlimb-specific supraspinal CSN innervation and specialization. A, B, CSNs in native PHATE space colored by supraspinal categories (A; Movie 3), and FL and HL identity (B). C, D, Divergent numbers of FL and HL CSNs transcriptionally defined by BICCN supraspinal structures. E, Photomicrographs show representative images of FL (top, cyan outline) and HL (bottom, magenta outline) td-tomato+ve CST projections after intersectional tracing in the brainstem (ET population). F, Projection density of FL-CSN-ET was significantly higher than in HL-CSN-ET (n = 4 FL traced mice, n = 3 HL traced mice, mixed-design ANOVA, **p = 0.003). G, Schematic illustrates FL axons (thick cyan lines) more densely innervating ET structures compared with HL axons (thin magenta lines). H, Number of DE genes comparing FL-CSN-ET with FL-CSN-IT, FL-CSN-CT, and NP combined (orange to black), and FL-CSN-ET versus IT (orange, upregulated, to green, downregulated), FL-CSN-ET versus CT (orange to purple) FL-CSN-NP combined (orange to black). H', These comparisons are repeated for the HL-CSN-ET population. I, Photomicrographs show FL and HL td-tomato+ve CST projections after intersectional tracing in the striatum (IT population). J, There was no significant difference in the projection density between ipsilateral and contralateral FL-CSN-IT and HL-CSN-IT projections (mixed-design ANOVA, p = 0.38). K, Schematic illustrates FL axons and HL axons innervating IT structures (green). L, Number of DE genes comparing FL-CSN-IT with FL-CSN-ET, FL-CSN-CT, and FL-CSN-NP combined (green to black), and FL-CSN-IT versus ET (orange to green), FL-CSN-IT versus CT (purple to green), and FL-CSN-IT vs NP (green to red). L', These comparisons are repeated for the HL-CSN-IT population. M, Photomicrographs show images of FL and HL td-tomato+ve CST projections after intersectional tracing in the thalamus (CT population). VPL, Ventral posterolateral nucleus of thalamus. N, There was no significant difference in the projection density between FL-CSN-CT and HL-CSN-CT projections (mixed-design ANOVA, p = 0.085). O, Schematic illustrates FL axons and HL axons innervating the thalamus. P, Number of DE genes comparing FL-CSN-CT to FL-CSN-ET, FL-CSN-IT, and NP combined (purple to black), and FL-CSN-CT versus ET (purple to orange), FL-CSN-CT vs FL-CSN-IT (purple to green), and FL-CSN-CT versus FL-CSN-NP (purple to red). P', These comparisons are repeated for the HL-CSN-ET population. Q, Photomicrographs show images of FL and HL td-tomato+ve CST projections after intersectional tracing in M2 (NP population). R, Projection density of FL-CSN-NP was significantly higher than in HL-CSN-ET (n = 4 FL traced mice, n = 3 HL traced mice, mixed-design ANOVA, ***p < 0.001). S, Schematic illustrates FL axons more densely innervating NP structures compared with HL axons. T, Number of DE genes comparing FL-CSN-NP to FL-CSN-ET, FL-CSN-IT, and FL-CSN-CT combined (red to black), and FL-CSN-NP versus TT (red to orange), FL-CSN-IT versus NP (green to red), and FL-CSN-NP vs CT (purple to red). T', These comparisons are repeated for the HL-CSN-ET population. Scale bars: E, 500 µm; I, M, Q, 100 µm. Data tables from the scRNAseq containing the results of the DE analysis comparing FL CSN supraspinal subtypes and HL within each CSN supraspinal subtypes are provided in Extended Data Figure 7-1.
Figure 7-1
DE genes for FL and HL supraspinal CSN subtypes. Data tables from the scRNAseq containing the results of the DE analysis comparing FL CSN supraspinal subtypes and HL within each CSN supraspinal subtypes (Wilcoxon rank sum test with a Benjamini–Hochberg correction, p adjusted < 0.05). Download Figure 7-1, XLSX file.
CSN ET specialization
These data reveal that supraspinal structures are differentially innervated by FL-CSNs and HL-CSNs. To explore the transcriptional differences among these populations that may drive and/or support the unique supraspinal terminal distribution within these subdivisions we completed a differential gene expression analysis among FL and HL ET, IT, CT, and NP categories.
There were 964 genes upregulated and 484 genes downregulated in FL-CSN-ET cells compared with FL-CSN-IT, FL-CSN-CT, and FL-CSN-NP combined (Fig. 7H). A more granular comparison among subcategories revealed that between FL-ET and FL-IT, 477 genes were upregulated and 1089 downregulated, and QIAGEN IPA revealed that GPCR (G-protein-coupled receptor) signaling and CREB signaling were among the top pathways dysregulated. Between FL-ET and FL-CT, 465 genes were upregulated, and 1065 genes were downregulated, with GPCR signaling and RHOGDI signaling among the top enriched pathways. Between FL-ET and FL-NP, 630 genes were upregulated, and 1808 genes were downregulated, with CREB signaling Gαi signaling among the dysregulated pathways (Fig. 7H).
For the HL-CSN-ET population, 807 genes were upregulated and 588 genes downregulated in HL-CSN-ET cells compared with HL-CSN-IT, HL-CSN-CT, and HL-CSN-NP combined (Fig. 7H'). Comparisons among the subcategories revealed that between HL-ET and HL-IT, 538 genes were upregulated, and 821 genes were downregulated, with axon guidance signaling and CREB signaling among the top pathways dysregulated. Between HL-ET and HL-CT, 735 genes were upregulated, and 1065 genes were downregulated, with Gai signaling and RHOGDI signaling among the top enriched pathways. Between HL-ET and HL-NP, 1234 genes were upregulated, and 2193 genes were downregulated, with spliceosomal cycle signaling and EIF2 signaling among the top pathways (Fig. 7H').
CSN IT specialization
We found that 282 genes were upregulated and 228 genes downregulated in FL-CSN-IT cells compared with FL-CSN-ET, FL-CSN-CT, and FL-CSN-NP combined (Fig. 7L). A more specific comparison among subcategories revealed that between FL-IT and FL-ET, 1089 genes were upregulated, and 477 genes were downregulated, with valine degradation and isoleucine degradation among the top dysregulated pathways. Between FL-IT and FL-CT, 234 genes were upregulated, and 359 genes were downregulated, with Rho family signaling and SEMA signaling among the top dysregulated pathways. Between FL-IT and FL-NP, 820 genes were upregulated, and 1155 genes were downregulated, with glutamate receptor signaling and putrescine degradation among dysregulated pathways (Fig. 7L).
For the HL-CSN-IT population, 388 genes were upregulated, and 488 genes downregulated in HL-CSN-IT cells compared with HL-CSN-ET, HL-CSN-CT, and HL-CSN-NP combined (Fig. 7L'). Comparisons among the subcategories revealed that between HL-IT and HL-ET, 821 genes were upregulated, and 538 genes were downregulated, with axon guidance signaling and valine degradation among the top dysregulated pathways. Between HL-IT and HL-CT, 523 genes were upregulated, and 749 genes were downregulated, with SEMA signaling and RHOGDI signaling among the top enriched pathways. Between HL-IT and HL-NP, 1377 genes were upregulated, and 1758 genes were downregulated, with zymosterol biosynthesis and putrescine degradation among the top pathways (Fig. 7L').
CSN CT specialization
We found that 201 genes upregulated, and 386 genes downregulated in FL-CSN-CT cells compared with FL-CSN-ET, FL-CSN-IT, and FL-CSN-NP combined (Fig. 7P). A more specific comparison among subcategories revealed that between FL-CT and FL-ET, 1065 genes were upregulated, and 465 genes were downregulated, with circadian rhythm signaling and synaptogenesis signaling among the top pathways dysregulated. Between FL-CT and FL-IT, 359 genes were upregulated, and 234 genes were downregulated, with axon guidance signaling and CREB signaling in neurons among the top dysregulated pathways. Between FL-CT and FL-NP, 757 genes were upregulated, and 1068 genes were downregulated, with glutamate receptor signaling and GABA receptor signaling among the dysregulated pathways (Fig. 7P).
For the HL-CSN-CT population, 425 genes were upregulated and 805 genes downregulated in HL-CSN-IT cells compared with HL-CSN-ET, HL-CSN-CT, and HL-CSN-NP combined (Fig. 7P'). Comparisons among the subcategories revealed that between HL-CT and HL-ET, 1370 genes were upregulated, and 735 genes were downregulated, with synaptogenesis signaling and serotonin receptor signaling among the top pathways dysregulated. Between HL-CT and HL-IT, 749 genes were upregulated, and 523 genes were downregulated, with circadian rhythm signaling and the opioid signaling pathway among the top enriched pathways. Between HL-CT and HL-NP, 1109 genes were upregulated, and 1500 genes were downregulated, with valine degradation and glutamate receptor signaling among the top dysregulated pathways (Fig. 7P').
CSN NP specialization
We found that 709 genes upregulated and 1163 genes downregulated in FL-CSN-NP cells compared with FL-CSN-ET, FL-CSN-IT, and FL-CSN-CT combined (Fig. 7T). A more specific comparison among subcategories revealed that between FL-NP and FL-ET, 1808 genes were upregulated, and 630 genes were downregulated, with GADD45 signaling and senescence signaling among the top pathways dysregulated. Between FL-NP and FL-IT, 1155 genes were upregulated, and 630 genes were downregulated, with cAMP signaling and GPCR signaling among the top dysregulated pathways. Between FL-NP and FL-CT, 1068 genes were upregulated, and 757 genes were downregulated, with cAMP signaling and GPCR signaling among the top dysregulated pathways. (Fig. 7T).
For the HL-CSN-NP, 1096 genes were upregulated, and 1706 genes were downregulated in HL-CSN-NP cells compared with HL-CSN-ET, HL-CSN-IT, and HL-CSN-CT combined (Fig. 7T'). Comparisons among the subcategories revealed that between HL-NP and HL-ET 2193 genes were upregulated, and 1234 genes were downregulated, with molecular mechanisms in cancer signaling and axon guidance signaling among the top pathways dysregulated. Between HL-NP and HL-IT 1758 genes were upregulated, and 1377 genes were downregulated, with axon guidance signaling and the serotonin receptor signaling pathway among the top enriched pathways. Between HL-NP and HL-CT 1500 genes were upregulated, and 1109 genes were downregulated, with GPCR signaling and CREB signaling in neurons among the top dysregulated pathways (Fig. 7T'). The entire DE lists are available in Extended Data Figure 7-1.
Discussion
In this study we used a combination of viral retrograde labeling from the cervical and lumbar spinal cord, FACS, and scRNAseq to annotate a comprehensive transcriptional index of adult CSNs. Using a refined extraction, dissociation, and FACS pipeline, we were able to maintain cytoplasmic integrity of large-diameter pyramidal neurons for whole-cell sequencing (Fig. 1; Golan and Cafferty, 2021) of intact adult CSNs (Figs. 2, 3). Using this approach, we identified 849 genes that were significantly differentially expressed between CSNs and nonspinally projecting L5 neurons. We built a linear classifier to confirm the veracity of this constellation of genes to identify CSNs with improved fidelity over a classifier trained on a set of established L5 genes (Fig. 4C). Furthermore, we found that 1280 genes were differentially expressed between FL and HL CSNs. These data allowed us to train a limb-specific CSN linear classifier based on the top 100 DE genes between FL and HL (Fig. 5N). These gene lists can be used in combination to identify FL and HL CSNs from a mixed population of neurons in any single-cell sequencing dataset from the intact adult mouse brain. Additionally, we found that adult CNSs do not fit exclusively into the BICCN-defined extratelencephalic category (also known as pyramidal tract; Callaway et al., 2021). Rather, they express a profile aligning with all nonspinally projecting categories IT, CT, and NP (Fig. 6A–D). Furthermore, these new CSN subcategories CSN-ET, CSN-IT, CSN-CT, and CSN-NP each display a unique gene expression profile (Fig. 6E–J). Finally, we confirmed that CSNs innervate supraspinal structures consistent with projection categories ET, IT, CT, and NP, highlighting the capacity of the CST to influence supraspinal motor processing.
Molecular specialization among forelimb and hindlimb CSNs
CSN identity has been explored during development by combining retrograde tracing from the pons and upper cervical spinal cord postnatally with temporal bulk microarray (Arlotta et al., 2005) and RNA sequencing (Molyneaux et al., 2015) of FACS-purified traced neurons. Authors identified Ctip2, Crym, and Crim1 as CSN-specific genes up to P14, leaving the CSN transcriptome temporally and spatially incomplete. Building on these transformative studies, we combined spinal retrograde labeling with scRNAseq to comprehensively profile FL and HL CSNs. We were able to identify and anatomically confirm robust CSN-specific and limb-specific genetic signatures that expand beyond previously identified single genes. Our findings revealed a novel constellation of genes that can be used to assign CSN identity with greater fidelity/confidence over previously identified CSN markers.
Indeed, the power of CSN subtype identity was previously highlighted in a study that sought to explore transcriptional divergence among rostral- (FL) and caudal (HL)-projecting CSNs during postnatal patterning. Authors identified Klhl14 and Crim1 as crucial for segment-specific axon targeting by completing bulk microarray sequencing of neurons in the rostral forelimb area (RFA) and comparing them with neurons in caudal forelimb area and hindlimb M1 regions (Sahni et al., 2021a,b). To elaborate on the findings, we explored limb-specific differences in the adult. Neither Klhl14 or Crim1 are significantly enriched in FL or HL CSNs in adulthood, likely reflecting downregulation of patterning machinery in the adult. However, here, CSNs were harvested for scRNAseq from +1.0 to −0.5 mm relative to bregma, potentially undersampling CSNs from the RFA.
Functional specialization among forelimb and hindlimb CSNs
Functional subtypes among M1 CSNs were previously described using slice electrophysiology and morphologic analyses (Oswald et al., 2013). Unbiased clustering found that among cortical neurons, FL- and HL-projecting CSNs fired action potentials at the fastest rate and with the most regular pattern, and that burst firing was exclusive to FL CSNs. This is consistent with FL-enriched expression of Cacng7, a type II transmembrane AMPA receptor regulatory protein that enhances glutamate-evoked currents from glutamate receptors (Kato et al., 2007).
There is also evidence that FL and HL CSNs play separate roles in sensorimotor function. A previous study showed that the lumbar CST functions exclusively to gate sensory input via primary afferent depolarization, and CST-mediated motor commands to the lumbar spinal cord are not direct but rather use polysynaptic circuits in the upper cervical cord (Moreno-Lopez et al., 2021). Our data support this finding as mice deficient in the HL-enriched gene Slc16a2 do not display gross motor abnormalities but do exhibit shorter latency paw withdrawal during a hot plate test, suggestive of hyperalgesia (Wirth et al., 2009; Moreno-Lopez et al., 2021).
This spinal enlargement level transcriptional resolution resolves broad FL and HL-level functional differences; however, there could also be within and between spinal segment specificity. Indeed, it was recently shown that different phases of FL movement are coordinated by CSNs projecting to discrete spinal segments (Ueno et al., 2018).
Plasticity of forelimb and hindlimb CSNs
Data from our laboratory support the notion that FL and HL CSNs require independent interventions to stimulate functional axon growth. Previously, we completed bulk RNAseq of retrogradely traced intact CSNs after contralateral PyX-induced functional sprouting into the denervated side. DE analyses revealed that members of the LPAR1-interacting pathway and the 3-phosphoinositide degradation pathway were enriched in sprouting CSNs. Viral overexpression of candidate genes within these pathways in intact CSNs after PyX resulted in significant sprouting of intact CST axons into the denervated side of the spinal cord (Fink et al., 2017; Kauer et al., 2022). Assessment of fine motor skills showed that the lesioned FL performed better than the lesioned HL, suggesting that our treatments were having a greater impact in FL CSNs. These data support limb-specific therapeutic interventions as the progrowth factors we identified in our original in vivo screen were identified via retrograde tracing from the cervical spinal cord. We hypothesize that tracing from the lumbar cord after contralateral PyX may reveal a unique set of proaxon growth modulators.
CSNs segregate by supraspinal collateralization
In addition to FL and HL specialization, CSNs can also be classified based on their supraspinal projections. We found that FL-CSNs innervate ET and NP structures more densely than HL-CSNs. Furthermore, each of these subdivisions can also be defined transcriptionally, suggesting they each influence fine motor processing in a previously underappreciated fashion. Their expression of a range of axon guidance ligands and receptors, suite of G-protein coupled receptors and synaptogenesis signaling components suggest they may each support a unique response to injury, disease, and aging. In alignment with these findings, previous work using scRNAseq of adult motor regions established a taxonomy of neuronal cell types (Tasic et al., 2016, 2018; Callaway et al., 2021; Yao et al., 2021). To describe the diversity observed among neuronal cell types spanning the cortical depth in M1, cell-type-specific genes describing IT, ET, NP, and CT neurons were identified (Yao et al., 2021). The ET designation was confirmed in part by retrograde tracing from the pons and medulla and scRNAseq of ALM neurons (Economo et al., 2018; Tasic et al., 2018). However, these studies did not include tracing from the spinal cord, thus rendering a full ET atlas incomplete. Previous studies have combined scRNAseq with single-cell reconstructions from sparsely labeled L5 neurons. These studies identified a subtype of L5 ET neurons with projections in the spinal cord, in addition to many supraspinal regions (Callaway et al., 2021; Peng et al., 2021). Furthermore, axonal barcoding of the CST revealed that many CSNs innervate a range of supraspinal structures in addition to their terminal fields in the spinal cord (Hausmann et al., 2022). These data are consistent with our findings that CSNs terminate in multiple supraspinal locations. By aligning sequenced CSNs with previously published M1 sequencing datasets, we found that CSNs are not all characterized solely as ET cells but rather span all four designated subtypes.
There is also evidence that CSN supraspinal connectivity reflects functional differences among CSNs. Studies exploring FL-CSNs with supraspinal collaterals in the striatum found that even among this one CSN subclass, there was anatomic and functional heterogeneity. Using intersectional polysynaptic tracing, authors found that CSNs synapse onto both D1 and D2 neurons in the striatum and GAD2 and calbindin D28K interneurons in the spinal cord. By combining two-photon imaging of FL-CSNs labeled via retro-AAV-GCaMP6f infusion into the cervical cord with a lever-pulling task, authors found that striatal CSNs display three unique activity patterns during pulling onset, pulling offset, and throughout the task (Nelson et al., 2021). This study provides robust evidence that supraspinal CST connectivity is critical to refining motor function, further necessitating a need for a comprehensive evaluation of the functional phenotypes of these molecularly divergent CSN subtypes.
In summary, from these data we conclude that CSNs in M1 are molecularly unique among nonprojection layer 2/3 and layer 5 neurons and can be further subdivided based on the location of their spinal rostrocaudal and supraspinal terminal arbors. Leveraging this molecular heterogeneity provides a powerful opportunity to explore the transcriptional machinery that drives development, maintains homeostasis, and actuates reactive changes after trauma and disease within this critical motor pathway. Comprehensive transcriptional insight in these physiological and maladaptive processes will enable refined therapeutic interventions to repair the damaged CST with enhanced precision.
Footnotes
This work was supported by National Institutes of Health–National Institute of Neurological Disorders and Stroke Grants R01NS09593001 and R21NS108053.
The authors declare no competing financial interests.
- Correspondence should be addressed to William B. Cafferty at william.cafferty{at}yale.edu
