Identification of the neuronal types that form the specialized circuits controlling distinct behaviors has benefited greatly from the simplicity offered by zebrafish. Electrophysiological studies have shown that additional to connectivity, understanding of circuitry requires identification of functional specializations among individual circuit components, such as those that regulate levels of transmitter release and neuronal excitability. In this study we use single cell RNA sequencing (scRNAseq) to identify the molecular bases for functional distinctions between motoneuron types that are causal to their differential roles in swimming. The primary motoneuron (PMn) in particular, expresses high levels of a unique combination of voltage-dependent ion channel types and synaptic proteins termed functional ‘cassettes’. The ion channel types are specialized for promoting high frequency firing of action potentials and augmented transmitter release at the neuromuscular junction, both contributing to greater power generation. Our transcriptional profiling of spinal neurons further assigns expression of this cassette to specific interneuron types also involved in the central circuitry controlling high speed swimming and escape behaviors. Our analysis highlights the utility of scRNAseq in functional characterization of neuronal circuitry, in addition to providing a gene expression resource for studying cell type diversity.
In zebrafish, primary motor neurons (PMNs) control escape movements, and a more heterogeneous population of secondary motor neurons (SMNs) regulate the speed of rhythmic swimming. Using single cell RNA sequencing (scRNAseq), the authors have obtained compelling evidence that PMNs, and two types of interneurons innervating them, express a set of three genes encoding voltage-gated ion channels enabling rapid firing. The PMNs also express high transcript levels of proteins involved in exocytosis, which would be expected to support rapid neurotransmitter release. These results will be important for those working on spinal cord function and zebrafish genomics/transcriptomics.
Functional and anatomical studies of spinal circuitry among the vertebrates have formed the basis of our understanding of control over stereotypic movements (Goulding 2009, Grillner & Jessell 2009). Investigation into movement control continues to benefit from larval zebrafish, which represent a greatly simplified system for resolving the underlying spinal circuitry (Fetcho & Liu 1998, Drapeau et al 2002, Lewis & Eisen 2003, Fetcho et al 2008). Study of circuitry in zebrafish also provides the unique opportunity to trace locomotory circuitry from sensory initiation to the final motor output (Fetcho 1991, Koyama et al 2011, Fidelin & Wyart 2014, Berg et al 2018). Fortuitously, despite the evolutionary distance between fish and mammals, many classes of spinal interneurons involved in movement control are conserved between species, heightening the potential significance of zebrafish circuitry analysis (Grillner 2003, Goulding 2009).
As a new approach towards circuitry analysis, we turned to scRNAseq. For this purpose, we developed a method for isolation and dissociation of spinal cords from 4 days post fertilization (dpf) zebrafish, an age at which much of the anatomy and physiology of swim control has been published (Bhatt et al 2007, McLean et al 2007, Fetcho et al 2008, Liao & Fetcho 2008, McLean et al 2008, Satou et al 2009, Menelaou & McLean 2012, Wang & Brehm 2017, Bello-Rojas et al 2019, Menelaou & McLean 2019, Kishore et al 2020, Satou et al 2020, Wen et al 2020). Our analysis has provided identification of major classes of neuronal and glial types. Of those neuronal types, many are known to be directly involved in locomotory circuitry either in mouse, zebrafish, or both. Using the markers revealed by the transcriptome analysis, we validated a number of key interneuron types, previously shown to be present in zebrafish locomotion circuit. In addition, we identified a new excitatory interneuron type that has a unique transmitter phenotype among interneurons.
Our interest in applying scRNAseq methodology to spinal neurons of zebrafish went beyond identifying transcriptional markers. Rather, we sought to use the neuron-specific transcriptomics to mine for distinctions in the ion channels and synaptic players that are associated with specialized circuit function. Many physiological studies have indicated fundamental differences in excitability and synaptic transmission among neurons controlling escape behavior versus those controlling slow to moderate rhythmic swimming speeds (Bhatt et al 2007, Liao & Fetcho 2008, McLean et al 2008, Satou et al 2009, Menelaou & McLean 2012, Menelaou & McLean 2019, Kishore et al 2020, Satou et al 2020). The signaling molecules causal to these differences among spinal neuron types are largely unknown. To address this outstanding question, we compared different types of motor neurons, where two well studied subtypes serve different roles. The PMns control the strongest contraction that provide for the single powerful bend, initiating escape, and participate in rhythmic swimming at only the highest speeds. By contrast, the SMns collectively regulate the range of slower rhythmic swimming (McLean et al 2007, Gabriel et al 2011, Wang & Brehm 2017). Paired patch clamp recordings, possible only in zebrafish, have characterized these Mn types as functional bookends. The PMns fire action potentials at ultrahigh frequency and the neuromuscular synaptic responses are able to follow with fidelity, whereas the SMns respond with lower frequency action potential firing and synaptic transmission at the neuromuscular junction and are subject to frequent failures (Wang & Brehm 2017, Wen et al 2020).
Our scRNAseq analysis, using Mn subtype specific markers that we validate in this paper, has provided candidates for serving these functional differences. First, the PMns express a trio of unique voltage-dependent ion channels, seen only at very low levels in the SMns, that are tailored for high frequency transmission. The same fast ion channel cassette also enriched in two well-characterized interneuron types that control firing of the PMns and are directly involved in the fast swimming and escape behavior. Second, the PMns also express significantly higher transcript levels of several key proteins involved in exocytosis, collectively termed a synaptic cassette. Thus, scRNAseq offers a new means to interrogate spinal circuitry through assignment of specialized signaling molecules. This application may also prove useful to understanding specialized circuits within the central nervous system (CNS) of mammals.
Transcriptional profiling of the larval zebrafish spinal cord
scRNAseq was performed on 4 dpf larval zebrafish, an age corresponding to many studies of spinal circuity and electrophysiological analysis of neuronal control over swimming behavior (Bhatt et al 2007, Liao & Fetcho 2008, McLean et al 2008, Satou et al 2009, Menelaou & McLean 2012, Wang & Brehm 2017, Menelaou & McLean 2019, Kishore et al 2020, Satou et al 2020, Wen et al 2020). For each of two duplicate experiments, approximately 150 spinal cords were isolated, followed by dissociation into single cells. In each case, ~10,000 to 15,000 spinal cord cells were sequenced to a depth of 40,000 reads per cell and aligned to an improved zebrafish reference genome that is more inclusive of 5’ and 3’ untranslated regions (Lawson et al 2020);see Methods). After applying standard quality control filters, data integration resulted in a total of 11762 spinal cord cells from the combined datasets.
Graph-based, unsupervised clustering of the entire spinal cord transcriptome gave rise to 30 clusters, 27 of which were readily identifiable on the basis of established neuronal or glial markers forming the two broad categories of cell types (Fig. 1A). The neuronal markers were elavl4 and snap25a (Fig. 1B) and glia markers were gfap, slc1a2b, myrf and sox10 (Fig. 1C). The glial cells (46% of total cells) fell into two non-overlapping groups, corresponding to astrocytes/radial glia (gfap+/slc1a2b+), and oligodendrocytes (sox10+/myrf+) (Fig. 1C). Both glial types were composed of multiple clusters, indicating further diversity. The remaining 3 clusters (clusters 6, 11, and 17, corresponding to 11% of total cells) showed mixed expression of neuronal and glial markers. It’s unclear whether this population reflects a true cell type or instead a group of unremoved doublets. Given that these clusters could not be assigned to either neurons or glia with confidence they were excluded from the subsequent analysis. Additionally, while the glia data are available as a resource, they were not analyzed further in this study.
Transcriptional profiling of spinal cord neurons
The neuronal population identified by elavl4+/snap25a+ expression was re-grouped into 33 clusters using Seurat (Fig. 2A). To assign individual neuronal identities to these transcriptome clusters, we used a combination code that relied on the co-expression of neurotransmitter biosynthesis/transporter genes along with differentially expressed marker genes (DEGs; Fig. 2; Table 1; Supplementary Table 1). The first code provided clear separation of neuronal populations into four principal categories; glutamatergic (slc17a6a+/slc17a6b+), glycinergic (slc6a5+), GABAergic (gad1b+/gad2+) or cholinergic (chata+/slc18a3a+) types (Fig 2B, Table 1). The second code relied on not only established markers for both zebrafish and mouse spinal neurons, but also new markers identified in this study. Generating the list of candidate markers was aided by previous studies that profiled the neurotransmitter identity of morphologically distinct neurons in the larval spinal cord (Bernhardt et al 1990, Hale et al 2001, Higashijima et al 2004a, Higashijima et al 2004c). In addition, since many aspects of the transcriptional program that establish the spinal neuronal circuit have been shown to be evolutionarily conserved among vertebrates (Kiehn & Kullander 2004, Goulding 2009, Grillner & Jessell 2009), we cross-referenced our data with recent mouse spinal cord sequencing data to search for homologous marker genes (Delile et al 2019, Blum et al 2021).
Using the combination code, we identified specific classes of sensory neurons, Mns and interneurons (Fig. 2C, Table 1). Two types of sensory neurons, the glutamatergic Rohon-Beard (RB) cells and GABAergic Kolmer-Agduhr (KA; also referred to as cerebrospinal fluid-contacting neurons (CSF-cN) cells) are each represented by two distinct clusters (Fig. 2A & C). Cholinergic clusters formed a prominent group corresponding principally to Mns. Clusters of interneurons corresponded to the inhibitory v1, v2b and v2s, and excitatory v0v, v2a and dl1, and dl2 types (Fig. 2A & C). The relative abundance of the cell types associated with individual clusters (Table 1) were consistent with those published for in vivo labeling experiments (Higashijima et al 2004b, Kimura et al 2006, Ampatzis et al 2014, Gerber et al 2019, Satou et al 2020), suggesting that our preparation protocol and analysis are robust and unbiased in recovering spinal cell populations. Of the 33 clusters in our neuronal dataset we were able to assign identities to 22 clusters with confidence. There are a number of cell types described for in larval zebrafish, such as the excitatory v3 interneuron or the inhibitory v0d interneuron (Satou et al 2020, Bohm et al 2022) that we were unable to identify, likely due to lower expression levels of canonical markers at this developmental stage. It is probable that these interneuron populations are present in unassigned clusters.
Our analysis also resolved subtypes within neuronal classes. For example, the two clusters were assigned as KA neurons, that shared a common set of markers for the cerebrospinal fluid contacting interneuron (pkd1l2a/pkd2l1; Fig. 3A). However, assignment to subtype can be made on the basis of differential expression of sst1.1 versus urp1, previously shown to label KA’ and KA’’ functional groups, respectively (Fig. 3A) (Djenoune et al 2017, Yang et al 2020).
Similarly, the glutamatergic v2a interneurons consist of multiple clusters that represented distinct subtypes. This class of interneurons has previously been divided into two sub-populations, Types I and II, based on morphology, molecular and functional heterogeneity (Bhatt et al 2007, McLean & Fetcho 2009, Ampatzis et al 2014, Menelaou et al 2014). One molecular feature differentiating the two types is the higher expression level of both the vsx2 and shox2 marker genes in Type I compared to Type II v2a (Kimura et al 2006, Menelaou et al 2014, Hayashi et al 2018, Menelaou & McLean 2019). In our dataset both the vsx2 and shox2 expression patterns differed among the v2a clusters (Fig. 3B). The two clusters with strong vsx2 and shox2 expression likely represent Type I v2a neurons, while the third cluster with weak vsx2 expression and no shox2 expression likely represents Type II (Fig. 3B). For this cluster, we further identified gjd2b, the gene encoding the gap junction protein connexin 35.1 δ subunit as an additional marker gene (Fig. 3B). This finding is consistent with previous immunohistochemistry in adult fish demonstrating selective expression of gjd2 in Type II v2a axons (Carlisle & Ribera 2014, Pallucchi et al 2022). Subtype assignment based on expression patterns of these markers was validated using in situ hybridization in larval Tg(vsx2: Kaede) fish (Fig. 3C). Expression of the Kaede fluorescent protein driven by vsx2 promoter labels Type I v2a with strong fluorescence compared to weakly labeled Type II (Fig. 3C). Probes against gjd2b colocalized only with those neurons exhibiting weak fluorescence, confirming its specific expression in the Type II v2a subtype (Fig. 3C). The two subtypes of v2a interneurons both make direct connections with the PMns and are recruited during high speed swimming but differ in terms of connection strength and the type II v2a interneuron is recruited more effectively across the range of high swimming speeds (Kimura et al 2006, Bhatt et al 2007, McLean & Fetcho 2009, Menelaou & McLean 2019). The availability of transcriptomes for these neuronal sub-groups provides an opportunity to further mine for differential molecular features responsible for the functional distinction.
In addition to validating established markers for neuronal types and subtypes, our scRNAseq analysis revealed novel marker genes for identifying interneuron transcriptomes, as shown for the Commissural Local (CoLo) and Dorsal Longitudinal Ascending (DoLA) interneurons.
CoLo interneurons provide the fast contralateral inhibition necessary for a successful escape response through direct contact with the PMns (Fetcho & Faber 1988, Liao & Fetcho 2008, Satou et al 2009, Kishore et al 2020). CoLo sends axons ventrally that turn to the contralateral side and split into short ascending and descending projections (Higashijima et al 2004a, Higashijima et al 2004c, Liao & Fetcho 2008, Satou et al 2009, Kishore et al 2020, Satou et al 2020) (Fig. 4A1). We found a small glycinergic cluster that expressed high levels of chromogranin A (chga) (Fig. 4A2) that was distinct from the chga-enriched cholinergic cluster later assigned to the PMn type (Fig. 5D). Labeling using chga in situ hybridization probes revealed a distinct large neuron located rostrally in each hemi-segment in addition to several large-sized motor neurons (Fig. 4A3, see also Fig. 5D). The stereotypical location and the one per hemi-segment stoichiometry of these chga labeled interneurons is consistent with that of CoLos (Liao & Fetcho 2008, Satou et al 2009, Kishore et al 2020). To validate its identity, we transiently labeled inhibitory interneurons with EGFP under the control of the dmrt3a promotor (Kishore et al 2020, Satou et al 2020). Single GFP labeled CoLos were screened based on morphology (Fig. 4A1) and subsequently shown to co-label with chga (Fig. 4A3).
DoLAs are GABAergic inhibitory interneurons located dorsally in spinal cord with well-described morphology (Bernhardt et al 1990, Higashijima et al 2004c) (Fig. 4B1). Our dataset indicated two GABAergic interneuron clusters, one of which corresponded to DoLA neurons with highly specific enrichment of pnoca expression (Fig. 4B2). In situ hybridization using pnoca probes strongly labeled a small number of distinct dorsal neurons located immediately ventral to the domain occupied by the sensory RB cells (Fig. 4B3). Sparse labelling of spinal neurons with a mCherry reporter showed that the neurons positive for pnoca projected long ascending axons with short ventral projections, as well as with occasional short descending axons, thus matching the morphological characteristics reported for DoLA interneurons at this stage (Fig. 4B4) (Bernhardt et al 1990, Higashijima et al 2004c, Wells et al 2010). The functional roles played by DoLA neurons have remained elusive. pnoca encodes Prepronociceptin, a precursor for several neuropeptides involved in multiple sensory signaling pathways (Martin et al 1998). Its highly specific expression in DoLA suggests that they might function as neuropeptide releasing neurons that modulate sensory functions in larval fish.
Importantly, our analysis also identified a cholinergic/glutamatergic spinal interneuron type not described previously in larval zebrafish. A single cluster in our dataset expressed cholinergic markers that include vesicular acetylcholine transporter vAChT (slc18a3a) and choline acetyltransferase (chata), but lacks all canonical Mn markers (mnx1/mnx2b/isl1) (Fig. 4C1). The cluster also expressed evx1 and evx2, markers associated with interneuron types in the v0 domain (Fig. 2C, Fig. 4C1) (Zagoraiou et al 2009, Juarez-Morales et al 2016). This transcriptional profile suggested that it represented the homologues to the premotor cholinergic v0c interneurons in mouse spinal cord (Zagoraiou et al 2009). Cholinergic interneurons have only recently been shown to be present in adult fish by immunochemistry staining (Bertuzzi & Ampatzis 2018). Similar to their mammalian counterparts, they play roles in modulating Mn excitability (Bertuzzi & Ampatzis 2018). Notably, v0c in larval zebrafish differs from the mouse counterpart on the basis of co-expression of cholinergic and glutamatergic transmitter genes (e.g., slc17a6b, Fig. 4C1).
We capitalized on the unique cholinergic phenotype of v0c among interneurons to provide in vivo labelling in the spinal cord. For this purpose, we injected a tdTomato reporter driven by the vAChT promoter to sparsely label cholinergic neurons in the 4 dpf spinal cord of Tg(mnx1: GFP) fish (Fig. 4C2). In addition to the Mns, we observed mosaic labeling of an interneuron type with distinct position and morphology. The soma was located near the dorsal boundary of the motor column (Fig. 4C2 & 3) and the axonal processes were either bifurcating (20 out of 37 cells) or purely descending (14 out of 37) or ascending (3 out of 37). The descending processes would reach lengths corresponding to multiple segments within the motor column (Fig. 4C2 & C3, average length > 6 segments). There was an overall lack of secondary branches, but enlargements reminiscent of synaptic boutons in close vicinity of Mn soma were observed along the neurites (Fig. 4C3). Multiple neurons of this type were observed in the same segment even with the sparse labeling approach, suggesting that there are likely to be numerous v0cs in each segment. The morphology and anatomic arrangement are consistent with a role in modulating Mn properties, as has been proposed for v0c in both adult zebrafish (Bertuzzi & Ampatzis 2018) and the mammalian homolog (Zagoraiou et al 2009).
Subclustering the Mns based on single-cell transcriptomes
We next focused on transcriptome comparison within Mn populations to examine their molecular heterogeneity. As a first approach, we isolated Mn clusters from the whole spinal cord dataset based on the overlap of two sets of marker genes (Fig. 5A). The first set was based on components of the cholinergic pathway that included slc18a3a, (Fig. 5A) and chata (Fig. 2B). The second set of Mn markers included the transcription factors mnx1, mnx2b and isl1 (Fig. 2C, Fig. 5A) (Appel et al 1995, Hutchinson & Eisen 2006, Zelenchuk & Bruses 2011, Asakawa et al 2012, Seredick et al 2012). Overall, ~27% of the profiled neuronal single cell transcriptomes corresponded to Mns (1354 cells) (Fig. 5A). As a complementary approach, we also generated samples enriched for a larger number of Mns to increase the statistical power of the clustering analysis. This was achieved through fluorescence activated cell sorting (FACS) of spinal cells prepared from the fluorescent transgenic fish line, Tg(mnx1:GFP), that broadly labels Mns (Flanagan-Steet et al 2005, Bello-Rojas et al 2019). Clustering analysis of the sorted data showed that >92% of the FACS sorted cells represent Mns based on the same canonical markers used for whole spinal cord, giving rise to 7790 cells for sub-clustering (Fig. 5B). We integrated the datasets from both the computationally isolated and experimentally purified Mn populations and performed clustering analysis using Seurat. A total of 9144 cells were grouped into 10 clusters (Fig. 5C). A small cluster (3.6% of total) expressing genes gfra1a and tbx3b (Supplemental Fig. 1) was almost entirely sourced from the FACS sorted datasets (>94%). In addition, it shared markers with a population of non-skeletal muscle Mns recently described in mouse spinal cord scRNAseq datasets (Blum et al 2021). Therefore, this cluster was not included in further comparisons among Mns that control skeletal muscle contraction.
The transcriptionally distinct clusters were next linked to previously known Mn types. Two broad types of Mns have been described for larval zebrafish, the PMn and the SMn, that are commonly distinguished by birth date, progenitor lineage and morphological features such as size, location, and periphery innervation pattern (Eisen et al 1986, Myers et al 1986, Menelaou & McLean 2012, Bello-Rojas et al 2019). There are four large-sized PMns (CaP, MiP, vRoP and dRoP) in each hemi-segment of the spinal cord, each innervating approximately one-quarter of axial muscle target field (Bello-Rojas et al 2019, Wen et al 2020). By contrast, 50-70 SMns are in a more ventral location, and display a gradient of sizes and functional properties (Myers et al 1986, Westerfield et al 1986, McLean et al 2007, Asakawa et al 2012, Wang & Brehm 2017, Wen et al 2020). We examined the top DEGs among the clusters (supplemental Table 2), and found two markers, chga and nr2f1a, with non-overlapping expression pattern that could reflect this broad classification (Fig. 5C). chga was enriched in one distinct cluster, while nr2f1a was present in the majority of the remaining cells (Fig. 5C). Taken together, these two mutually exclusive markers represented ~95% of Mn population.
In situ hybridization labeling with probes against chga and nr2f1a revealed spatially segregated Mn groups in the Tg(mnx1:GFP) fish. Specifically, chga probes labeled dorsal Mns that were large in size, consistent with them being the primary group (Fig. 5D). nr2f1a labeling was absent in these cells, but was distributed in larger number of smaller Mns that were more ventrally located in the motor column (Fig. 5D). These segregated patterns of expression strongly suggested that chga+ cluster represented the PMns, while nr2f1a marked the major SMn populations.
For further validation, we used the chga probes for in situ hybridization analysis in Tg(SAIG213A; EGFP) fish, in which a single PMn in each hemi-segment, the dorsal projecting CaP, was fluorescently-labeled among all the Mns (Muto & Kawakami 2011, Wen et al 2020). Strong chga labeling colocalized with the GFP labeled-CaP in each hemi-segment (Fig. 5E). We also labeled the MiP and RoP Mns in the spinal cord using the sparse labeling approach, and observed high level of chga signal in both PMn types using in situ hybridization (Fig. 5F). These results firmly established chga as the marker for PMns, leaving the chga negative clusters representing SMns (Fig. 5G). The number of cells associated with SMn clusters was ~18 fold in excess that of the PMn cluster, consistent with previous cell counts of SMns versus PMns in larval zebrafish (Eisen et al 1986, Myers et al 1986, McLean et al 2007, Menelaou & McLean 2012, Bello-Rojas et al 2019).
Diversity among SMns
In contrast to the single cluster associated with the PMn type, SMns were composed of multiple transcriptionally distinct groups. Approximately 95% of SMns fall into three groups that were distinguished by differential expression of three marker genes, foxb1b, alcamb and bmp16 (Fig. 6A). In situ hybridization using the three probes labeled SMns revealed two closely stratified layers that showed differences in dorsal-ventral positioning within the motor column (Fig. 6A). The foxb1b+ SMns occupied the more dorsal position of our labeled secondaries while both alcamb+ and bmp16+ SMns shared a more ventral location (Fig. 6A). We next tested alcamb+ and foxb1b+ SMns for potential correspondence to the different subsets of SMns labeled with two transgenic lines, Tg(isl1: GFP) and Tg(gata2: GFP) (Appel et al 1995, Meng et al 1997, Higashijima et al 2000). Labeling with in situ probes revealed that expression of alcamb overlapped with GFP+ SMns in the Tg(gata2: GFP) fish (Fig. 6B), while foxb1b co-localized with those in Tg(isl1: GFP) (Fig. 6C). Thus, foxb1b/alcamb expression distinguishes previously identified SMn subsets. The third cluster, corresponding to the bmp expressing class of SMns, were fewer in number and the targets are not known. However, bmp signaling has been linked to specification of slow muscle, raising the possibility that this muscle type receives innervation by bmp+ Mns (Kuroda et al, 2013).
Transcriptome comparisons for the PMn and SMn types
We performed differential expression analysis comparing PMn and SMn transcriptomes in order to identify candidate genes that might account for their functional distinctions previously established using electrophysiology. After applying thresholds to the Mn clusters, based on average levels of gene expression, p-value, and the proportion of cells expressing individual genes, we obtained a list of 508 candidates that showed at least a 30% difference in average expression levels between Mn types, 317 of which were higher in primaries. To guide our identification of genes that play potential roles in governing excitability and synaptic transmission differences between the two groups of Mns (Supplemental Table 2) we used Gene Ontology (GO) enrichment analysis for the PMn specific DEGs (The Gene Ontology Consortium, 2021) (Ashburner et al 2000). We focused on significantly enriched GO terms for DEGs that were associated with 3 broad categories of biological processes, synaptic function, ion channels/transporters and ion homeostasis, and ATP generation (Supplemental Fig. 2), due to their roles in excitability and synaptic transmission. Of the 37 Gene ontology terms identified for the PMn DEGs 24 fell into one of these 3 broad functionally relevant categories (Supplemental Fig. 2). We next manually annotated the molecular function of DEGs by referencing evidence-based database ZFIN (Bradford et al 2022) in order to identify specific functionally relevant genes differentially expressed between PMn and SMn. For the PMn type, DEGs encoding a large number synaptic proteins included those of the core exocytotic machinery such as isoforms of VAMP, SNAP25 and Syntaxin; regulators of exocytosis, including Synaptotagmins, NSF, Complexins and RIM; synaptic vesicle proteins, and synaptic structural proteins such as Synuclein, Bassoon and Piccolo (Fig. 7A). Comparisons of expression levels indicate that no synaptic genes were specifically enriched in SMns, but 29 were enriched in PMns compared to SMns (Fig. 7A). Additionally, 5 synaptic genes were shared at approximately equivalent levels between Mn types (Fig. 7A). Similarly, 8 different voltage dependent ion channel subunits were preferentially expressed in the PMn, with the highest levels corresponding to the P/Q calcium channel, the β4 sodium channel subunit and the Kv3.3 potassium channel. Four other ion channel types were shared by both SMns and PMns. As with the synaptic genes, no ion channel candidates were specifically enriched in the SMn type. The absence of unique ion channel and synaptic gene enrichment in SMns contrasted with the high levels of transcriptional factors and RNA binding proteins, neither of which are likely candidates for conferring functional distinctions between Mn types. Approximately two thirds of the 191 genes enriched in the SMns (118/191) were of these categories. The low levels of candidate ion channels and synaptic genes in the SMns might be expected on the basis of greatly reduced levels of synaptic transmission reflected in quantal content and release probability compared to PMn (Wang & Brehm 2017).
Of particular interest among the genes specific to the PMn type were the trio of highest expressing voltage-dependent ion channel types that have been individually linked in previous studies to augmenting transmitter release and/or AP firing rate. Those channel types include a voltage-dependent P/Q-type calcium channel α subunit (cacna1ab), a Kv3.3 potassium channel α subunit (kcnc3a), and a sodium channel b4 subunit (scn4ba) (Fig. 7B & C) (Eggermann et al 2011, Lewis & Raman 2014, Zhang & Kaczmarek 2016). The identification of the cacna1ab channel as a top DEG in PMns corroborated previous studies firmly establishing the P/Q-type as the presynaptic active zone calcium channel mediating Ca2+-dependent release specifically in the PMns (Wen et al 2013, Wen et al 2020). Enriched expression of the sodium channel b subunit scn4ba in PMns was validated by in situ hybridization (Fig. 7D). scn4ba probes specifically labeled dorsally-located Mns with large sizes in the Tg(mnx1:GFP) fish (Fig. 7D). These represented the PMns, as shown by co-labeling of individually labeled CaP, MiP and RoP Mns (Fig. 7D & E). The zebrafish kcnc3a has been shown to be expressed preferentially in PMn at embryonic ages (Issa et al 2011). We tested its expression at 4 dpf fish using immunohistochemical staining. A kcnc3 sub-type specific antibody efficiently labeled NMJ synaptic terminal marked by α-bungarotoxin (α-Btx, Fig. 7F2). Since the PMn and SMn axons track along with each other and form synapses with the same postsynaptic receptor clusters (Wen et al 2020), further resolution was needed to distinguish the PMn and SMn synaptic terminals. For this purpose, we ablated individual CaPs by transiently expressing the phototoxic KillerRed protein (Formella et al 2018), followed by light inactivation at 2 dpf. By 4 dpf, the ablation of CaP was complete as indicated by the absence of its soma and processes (Fig. 7F1). Kcnc3 antibody staining showed a specific signal reduction in synapses located in the ventral-most musculature (Fig. 7F2), the target field shared by the CaP and numerous other SMns. This result strongly supports the expression specificity of kcnc3a channel in the PMns.
These specific calcium, potassium and sodium channel isoforms co-express specifically among the four types of PMns. We term this collective unit of three different voltage-dependent ion channel types as a ‘channel cassette’. Over 55% of the cells in PMn cluster express the cassette, compared to <3% in the SMns over all (Fig. 7C). Together with other synaptic DEGs revealed by the analysis, our results suggested a molecular logic for the functional specialization in PMn synapses, such as strong release and high frequency firing that are uniquely associated with escape behavior.
Gene ensembles for additional components of escape circuitry
Further insights into the potential contribution the ion channel cassette to fast synaptic function specialization came from examination of their differential expression pattern in other circuitry components involved in high-speed swimming and escape. Those cell types include the CoLo inhibitory and v2a excitatory interneurons, both of which provide direct synaptic connections to the PMn type. Both type I and type II v2a neurons participate in high swim speeds, but the type II is more effectively recruited over the range of higher speeds (Hale et al 2001, Higashijima et al 2004c, Kimura et al 2006, Liao & Fetcho 2008, McLean et al 2008, Satou et al 2009, Ampatzis et al 2014, Menelaou & McLean 2019, Satou et al 2020). Additionally, electrophysiological studies have shown that type II v2a neurons fire faster and more reliably than type I v2a neurons (Menelaou & McLean 2019). As seen in the PMn type, the coexpression of the ion channel cassette is high in both type II v2a and CoLo interneuron types compared to other interneuron types (Fig. 8A). Some of those interneuron types expressed individual components of the cassette, most often kcnc3a or cacna1ab channel types, but not the full cassette. Indeed, the scn4ba sodium channel β subunit is the strongest delineator among the 3 for inclusion into the ion cassette classification (Fig. 8A & 7C).
Top DEGs identified for the PMns that encode synaptic proteins also showed higher levels of expression in Type II v2a and CoLo interneurons compared to other cell types, including vamp1a, vamp1b, syt2a, syt2b, snap25a, stxbp1a and cplx2l, all central players in transmitter release (Fig. 8B). This further strengthens the proposition that a synaptic cassette works collectively with the ion channel cassette to create a strong neuronal circuitry controlling the escape behavior and fast swimming behavior. These results are incorporated into a simplified model that potentially accounts for the differential circuitry that distinguish regulation of slow swimming from the power circuits (Fig. 8C).
Interest in the spinal circuitry controlling muscle movements among vertebrates remains high, especially from the standpoint of understanding inheritable human myasthenic disorders (Walogorksy 2012, Wen 2016a, Ono 2002, Wang 2008, Downes & Granato 2004). As a model for investigation into spinal circuits, zebrafish offers great simplicity when linking spinal circuitry to a specific behavior. Moreover, morphological, physiological and transgenic labeling techniques, have revealed a large repertoire of spinal neuronal types involved in locomotion in zebrafish, many of them sharing homologies with those in mouse (Grillner 2003, Goulding 2009, Grillner & Jessell 2009). In both mouse and larval zebrafish, the understanding of motility circuitry centers on the study of Mns which are the final pathway to movement. Unlike the numerous Mn subtypes in mammals, which are grouped according to different anatomical positions and muscle targets (Stifani 2014), there are two main classes in zebrafish which share the same fast muscle cells as well as the individual neuromuscular synapses (Eisen et al 1986, Myers et al 1986, Menelaou & McLean 2012, Bello-Rojas et al 2019, Wen et al 2020). To aid in the assignment of neuronal types to specific circuits, we performed scRNAseq analysis on larval zebrafish spinal cord. Our study revealed new markers for key components of the spinal circuity that likely contribute to regulation of different behavioral responses, along with identification of a new interneuron type.
Studies on Mn control of movement in zebrafish have focused on two idiosyncratic swimming behaviors; the single powerful tail bend initiating escape and the subsequent, generally less powerful, rhythmic swimming (Budick & O’Malley 2000, Thorsen et al 2004). The spinal circuits mediating these two behaviors have been the subject of many studies, which conclude that the SMns do not play important roles in behavioral escape but instead regulate swimming over a broad range of speeds. Consistent with this role they cover a lower range of AP firing and release much less transmitter per AP compared to the PMn type (McLean et al 2007, Wang & Brehm 2017, Wen et al 2020). Thus, PMns regulate the highest speed swimming and escape behavior. Indeed, the functional distinctions were our impetus to identify individual transcriptomes for the separate Mn types. Using a newly identified marker, chaga, we identified a single small cell cluster corresponding to the four PMns located within each hemisegment. The SMns, by contrast, corresponded to three different clusters based on the markers foxb1b, alcamb and bmp16. These three transcripts were used to localize SMn subtype populations within the spinal cord. The SMn somas containing foxb1b transcripts formed the most dorsal group which are positioned above the alcamb-labeled group in the motor column. The smallest cluster, labeled by bmp16, overlapped with the alcamb label. The relative position of these SMn clusters is consistent with studies showing a correspondence between the dorsal-ventral position within the spinal cord and control of rhythmic swim speed by the SMns. Sequential recruitment of more dorsally located SMns leads to the generation of increased power and faster swim speed (McLean et al 2007, Gabriel et al 2011, Wang & Brehm 2017). This gradient of power is determined through both firing pattern and by amount of transmitter release at the NMJ (Wang & Brehm 2017). Thus, it is plausible that the transcriptomic distinctions among SMns reflect their differential roles in swim speed determination and their functional distinctions in synaptic strength (Wang & Brehm 2017, Wen et al 2020).
Previous electrophysiological studies indicate that PMns fire APs at a high frequency and that transmitter release occurs with a short synaptic delay, high quantal content and a high release probability (Wen et al 2016b). By contrast, SMn synapses are much weaker and variable in firing properties, in keeping with the distinct behavioral roles of PMns and SMns (Wang & Brehm 2017, Wen et al 2020). Consistent with these distinctions, analysis of the transcriptomes for the two Mn types revealed large scale differences. The PMns expressed high levels of transcripts encoding ion channels and exocytotic machinery, all of which are generally not detected in SMns. Overall, the transcriptomic profile comparisons are consistent with the electrophysiological findings of greatly enhanced neuromuscular transmission for the PMns. In particular, the PMns had very high expression of an ion channel cassette formed by three different voltage-dependent ion channel types, each of which had been linked previously to either high release probability of transmitter or very high AP frequency. The cassette transcript with the most restrictive expression pattern encoded the β4 subunit of the voltage-dependent sodium channel NaV1.6. This subunit confers high frequency firing through a fast reversible block of the NaV1.6 channel pore. The blocking kinetics are sufficiently fast to enable the neuron to fire APs at frequencies exceeding the refractory period (Raman & Bean 1997, Grieco et al 2005, Lewis & Raman 2014, Ransdell et al 2017). A second highly enriched voltage-dependent channel in PMns is the Kv3.3 potassium channel. This ion channel type has been associated with high AP firing frequency, as well as with augmented transmitter release in mammalian neurons (Zhang & Kaczmarek 2016, Richardson et al 2022), both due to fast activation kinetics that result in fast repolarization of the AP. The zebrafish kcnc3a gene, encoding the Kv3.3 channel, gives rise to a transient potassium current that has both fast activation and inactivation making it well suited for its proposed role in shortening AP waveform and allowing high frequency firing (Mock et al 2010). The third cassette member, cacna1ab, encoding a P/Q-type calcium channel, was a top DEG in the PMn, in agreement with our previously published finding that the SMn expresses a different calcium channel isoform, most likely N-type based on sensitivity to specific conotoxin isoforms (Wen et al 2020). Mutations in the cacna1ab gene completely abolished AP-evoked release in PMns but left synaptic transmission in SMns intact (Wen et al 2020). As a result, mutant fish are unable to mount a fast escape response, but are still capable of normal fictive swimming (Wen et al 2020). The P/Q type calcium channel has been associated widely with synapses with high release probability (Iwasaki & Takahashi 1998, Wu et al 1999, Stephens et al 2001, Fedchyshyn & Wang 2005, Bucurenciu et al 2010, Eggermann et al 2011, Young & Veeraraghavan 2021), a feature, due in part, to a higher open probability during APs compared to the N-type counterparts (Li et al 2007, Naranjo et al 2015), thereby promoting calcium entry. We hypothesize that the collective actions of this ion channel cassette, with their unique biophysical properties, serve to mediate the escape response and the highest swim speed through ultrafast repetitive firing and maximal release of neurotransmitter. This trio of voltage dependent ion channels is also co-expressed in mammalian neuronal types that are involved in fast signaling, suggesting a conserved role for the cassette in fast behaviors. Examples include the auditory neuron calyx of Held (Iwasaki and Takahashi 1998, Midorikawa et al 2014, Richardson et al 2022, Kim et al 2010) and the pyramidal neurons of the brain (Grieco et al 2005, Akemann and Knopfel 2006, Hillman et al 1991).
A second cassette comprised of a set of genes involved in synaptic function was also revealed by comparing the transcriptomes between the PMns and SMns. The cassette components enriched in PMns were genes encoding isoforms of VAMP, Syntaxin, Synaptotagmin, SNAP25 and Complexin, all components associated with exocytosis and transmitter release. Unlike the ion channel cassette described above, which is strongly linked to synapses with high release probability and/or high firing rate, the actions by which these synaptic genes could differentially support strong and fast synaptic properties remains speculative. It has been suggested, however, at both fly NMJ and mammalian CNS, that the differential synaptic strength among synapses correlates with abundance of proteins involved in transmitter release (Holderith et al 2012, Peled et al 2014, Akbergenova et al 2018).
Further support for the idea that ion channel and synaptic cassettes both play direct roles in formation of specialized circuitry surrounding the PMn was provided by transcriptomic analysis of interneuron types known to interact specifically with the PMn and to be recruited during high speed swimming. Those include the Type II v2a excitatory interneuron and CoLo inhibitory interneuron forming the output pathway for the escape response (Bhatt et al 2007, Liao & Fetcho 2008, Satou et al 2009, Menelaou & McLean 2012, Menelaou & McLean 2019). As with the PMn, both interneuron types coexpress the triple ion channel cassette members at high levels compared to those interneurons that participate less during recruitment during high speed swimming. As shown for distinction among Mn types, the sodium channel β4 appears to be the most restrictive among the three ion channel types in conferring fast firing to the interneurons as well. Unlike the case for the PMn, the AP firing and transmitter release properties of these neuronal types are less well established. However, it is clear that the Type II v2a, in particular, can fire at high frequencies over 600 Hz (Menelaou & McLean 2019). The highly specific enrichment of these two gene cassettes, in neurons involved in escape behavior, suggests that they serve as an integral part of the molecular signature underlying functional specialization. It remains to be seen whether the gene cassettes we identified for zebrafish escape circuit represent a general transcriptional architecture plan to build synapses with great strength and speed in the CNS of higher vertebrates.
Finally, our scRNAseq analyses provides a resource for future identification of gene functions that are causal or associated with human disorders involving Mn dysfunction. In the context of myasthenic disorders in particular, zebrafish has provided a large number of animal models corresponding to human syndromes, including slow channel syndrome, episodic apnea, and rapsyn deficiency (Ono et al 2002, Wang et al 2008, Walogorsky et al 2012a, Walogorsky et al 2012b, Wen et al 2016a). Both myasthenic syndromes and amyotrophic lateral sclerosis (ALS) involve dysfunction at the level of the motor circuits (Ferraiuolo et al 2011). In the case of ALS, it is well known that fast motor neurons are selectively targeted for degeneration (Hadzipasic et al 2014, Nijssen et al 2017). Our analysis comparing transcriptional profiles between fast versus slow motor circuit components offers a new means for probing the transcriptional consequences of neuromuscular disease states.
Materials and Methods
Fish lines and husbandry
The transgenic line Tg(mnx1:GFP) was provided by Dr. David McLean (Northwestern University). Tg(vsx2:Kaede) was provided by Dr. Joseph Fetcho (Cornell University). Tg(SAIG213A;EGFP), Tg(islet1:GFP) and Tg(gata2:GFP) were maintained in the in-house facility. Zebrafish husbandry and procedures were carried out according the standards approved by Institutional Animal Care and Use Committee at Oregon Health & Science University (OHSU). Experiments were performed using larva at 4 dpf. Sex of the larva cannot be determined at this age.
Sparse labeling of spinal neurons
In most cases, we used the Gal4-UAS system to achieve mosaic expression by co-injection of two plasmids into single cell embryos: one containing Gal4 driven by cell type – specific promoters and the other containing fluorescent reporter genes under the control of UAS element. The mnx1 or vAChT promoter drove the expression in Mns (mnx1:Gal4 plasmid provided by Dr. McLean, Northwestern University and vAChT:Gal4 provided by Dr. Joe Fetcho, Cornell University). drmt3a promoter (drmt3a:Gal4 provided by Dr. Shinichi Higashijima, National Institute of Natural Sciences, Japan) was used for expression in glycinergic inhibitory interneurons. mCherry driven by HuC promoter was used to label spinal neurons, including the DoLA interneurons. Injected fish were screened on 3 dpf for sparse fluorescent neurons that could provide detailed morphology.
KillerRed mediated CaP ablation
We transiently expressed the phototoxic KillerRed protein in Tg(SAIG213A:EGFP) fish by injecting a plasmid expressing KillerRed driven by UAS promoter (Addgene plasmid # 115516; a gift from Marco Morsch). Fish with KillerRed expression in CaP were identified. Individual CaPs were ablated by light inactivation at 2 dpf by illuminating for 10 mins with a 560 nm laser set at high power. This completely bleached KillerRed fluorescence, and induced visible blebbing in CaP terminals. Fish were grown to 4 dpf, and the ablation was confirmed by the absence of GFP-labeled soma and neurites.
Whole mount immunocytochemistry
Whole mount immunohistochemistry was performed as described previously (Wen et al 2020). 4 dpf larvae were fixed in 4% paraformaldehyde at 4 °C for 4 hrs. Zebrafish Kcnc3a channel was labeled using a polyclonal antibody originally generated against the human Kcnc3 protein (ThermoFisher Scientific PA5-53714) at a concentration of 2.5 μg/ml. Polyclonal anti-GFP (Abcam ab13970) was used at 1 μg/ml. Alexa Fluor-conjugated secondary antibodies (Thermo Fisher Scientific) were used at 1 μg/ml. To mark the location of synapses, 1 μg/ml CF405s-conjugated α-Btx (Biotium) was included in the secondary antibody incubation to label postsynaptic acetylcholine receptors.
Whole-mount in situ RNA hybridization
Fluorescence in situ RNA hybridization (FISH) was performed on whole-mount 4 dpf larva using the multiplexed hybridization chain reaction RNA–FISH bundle (HCR RNA-FISH) according to the manufacturer’s instructions (Molecular Instruments) (Choi et al 2016). Probe sets for zebrafish alcamb, nr2f1a, chga, scn4ba, foxb1b, bmp16, gjd2b and pnoca were custom-designed based on sequences (Molecular Instruments) and used at 4 nM each. Fluorescent HCR hairpin amplifiers were used at 60 nM each to detect the probes. GFP and mCherry fluorescence in the transgenic lines and transient labeled neurons survived the FISH protocol with signal loss mostly limited to the periphery. Residual fluorescence was sufficient to mark the location of soma in the spinal cord without the need for additional amplification.
After staining, fixed larval were mounted in 1.5% low melting agarose and imaged on a Zeiss 710 laser-scanning microscope equipped with an LD C-Apochromat 40x/1.2 n.a. objective. Z-stacks of confocal images were acquired using Zen (Carl Zeiss) imaging software, and presented as either maximal intensity projection or single focal planes as indicated in the figures (ImageJ, National Institutes of Health).
Single cell suspension from spinal cord for scRNAseq
Single cell suspensions were prepared from Tg(SAIG213A;EGFP) fish for the two full spinal cord datasets, and from Tg(mnx1:GFP) for the two FACS-sorted Mn enrichment datasets. About 150 4 dpf larva were euthanized in 0.02% tricaine, and individually decapitated behind the hindbrain. They were incubated with 20 mg/ml collagenase (Life Sciences) in a buffer containing 134 mM NaCl, 2.9 mM KCl, 1.2 mM MgCl2, 2.1 mM CaCl2, and 10 mM Na-HEPES (pH 7.8) at 28 °C for 2 hr, with intermittent trituration using a p200 pipette aid at 0, 0.5 hr and 1 hr of the incubation. To release spinal cords from remaining tissue, the final triturations were done using fire-polished Pasteur pipettes with decreased opening sizes (300, 200, 100 μm respectively). Intact spinal cords were transferred to L15 media and washed 3 times with fresh media. The spinal cords were incubated with 0.25% trypsin solution (in 1xPBS containing 1 mM EDTA) at 28 °C for 25 min. The digestion was terminated by adding 500 μl stop solution (L15 with 1% fetal bovine serum). The tissue was collected by spinning at 400 g for 3 min at 4°C, washed once with L15 and resuspended in 200 μl of L15 media. Spinal cord cells were dissociated by triturating the digested tissue with fire-polished Pasteur pipettes with 80-100 μm opening. The solution was filtered through a 35 μm strainer into a siliconized collection tube. The suspension was examined on a microscope for cell count, Trypan blue staining based - viability test and proportion of dispersed single cells. Samples with a viability above 70% were used for sequencing.
Single cell suspension prepared from Tg(mnx1:GFP) fish were FAC sorted for EGFP+ cells using a 100 μm nozzle on a BD inFlux cell sorter (Flow Cytometry Shared Resource, OHSU). Cells were collected in 100 μl PBS containing 0.2% bovine serum albumin in a siliconized tube.
Single cell capture, cDNA synthesis, library preparation and sequencing
Single cell capture, cDNA synthesis and library preparation were performed by the Massive Parallel Sequencing Shared Resource at OHSU using the 10x Genomics Chromium v3.0 reagent kit. Single cell suspension for the two full spinal cord replicates targeted 10,000-15,000 cells. For the two FACS sorted samples 4,000-5,000 cells were targeted. Replicate samples were prepared from different clutches of animals. Libraries were sequenced on an Illumina NovaSeq 500 instrument to an average read depth of ~40,000 per cell.
Reference genome generation and alignment
Cellranger v6.1.1 (10X genomics) was used for the reference genome generation and alignment. Our reference genome was generated by modifying a preexisting reference genome, Lawson v4.3.2 (Lawson et al 2020), which we edited to add an EGFP sequence as an artificial chromosome and to correct a selection of gene names and 3’ untranslated region (UTR) annotations. A full account of all changes made to the Lawson reference genome can be found in Supplemental Table 3, with a representative example shown (Supplemental Fig. 3). Alignment to this reference genome was performed using Cellrangers count function with expect-cells set to the targeted number of cells for each replicate.
Preprocessing, normalization, clustering analysis and visualization
Count matrices were processed using the Seurat v4.0 package for R (www.satijalab.org/seurat/) (Hao et al 2021). Genes with expression in less than three cells were excluded from further analysis. Initial quality control was performed on each sample independently. Cells were kept for further analysis if they had a number of unique genes between 400 and 4000, UMI counts between 1500 and 9000, and <5% mitochondrial gene content. No further doublet removal methods were applied. Data was normalized using the Seurat Sctransform v2 package in R, generally following the procedure outlined in the Introduction to SCTransform v2 regularization vignette (Hafemeister & Satija 2019, Choudhary & Satija 2022). Principal components were calculated, followed by nearest neighbor graph calculation using ANNoy implemented through Seurat. Clustering used the Leiden community detection method implemented with the FindClusters function and the Leidenalg Python package (Traag et al 2019). To visualize clustered data sets we used t-distributed stochastic neighbor embedding (t-SNE) (van der Maaten & Hinton 2008) or uniform manifold approximation projection (UMAP) (McInnes & Healy 2018) implemented through Seurat.
Datasets were combined using Seurat’s integration pipeline (Stuart et al 2019), considering 9,000 variable genes. To examine the correspondence between duplicates, exploratory clustering was done in the combined datasets (Supplemental Fig. 4A & B). Intermixing of data points from different samples was inspected both visually, and by plotting the distribution against one another (Supplemental Fig. 4A & B). For the two full spinal duplicates, one cluster of cells exhibited a strong bias towards a single replicate, with over 90% of the cells sourcing from a single replicate (Supplemental Fig. 4A). These cells were not included in the combined dataset for downstream analysis to control for technical and biological variability. The FACS sorted duplicates had no obvious outliners (Supplemental Fig. 4B).
The combined motor neuron dataset was generated by integrating Mns extracted from the full spinal dataset and the FACS sorted dataset, based on expression of canonical Mn markers (Fig. 5A & B). While the relative populations of SMns did vary between the two sample sources, there was no Mn subtypes that could not be identified independently in both methods of Mn isolation (Supplemental Fig. 5A). Differential expression analysis performed independently with either Mn source yielded highly reproducible sets of the DEGs (Supplemental Fig. 5B), further validating the results from the combined dataset.
Annotation of cell clusters
Specific cell types in the neuronal subset of the spinal data were annotated on the basis of significantly differentially expressed genes. Exploratory clustering analysis was first conducted to filter out cells not of spinal cord origin. One small cluster in the whole spine dataset (0.6% of cells), expresses tph2/ucn3l/slc18a2 at high level compared to all the rest of the clusters (Supplemental Fig. 6), marker genes for serotonergic raphe nucleus in the hindbrain (Oikonomou et al 2019, Bradford et al 2022). They reflected a trace amount of hind brain tissue during spinal cord dissection, and were removed from the analysis. Contaminating muscle cells were also removed based on expression of my1pfa/tnnt3b/actc1b myosin/troponin/actin genes (Supplemental Fig. 6). Overall, these contaminants accounted for ~2.1% of the total cell population isolated from the spinal cord. No cells in the FACS sorted data set were identified as originating from contamination from outside the spinal cord.
Data analysis and statistical tests
Differential gene expression was calculated using a Wilcoxon rank sum test implemented with the FindMarkers or FindAllMarkers functions in Seurat. Genes were considered to be enriched in a cluster if they had a log2 fold change > 0.38 (corresponding to > 30% enrichment in expression level), expression in at least 30% of cells in one of the groups being compared, and a Bonferroni adjusted p value of < 10-5. A more conservative adjusted p value of < 10-10 was used for comparison between PMns and SMns.
Gene ontology analysis
DEGs comparing PMns and SMns were used to generate lists of gene ontology (GO) terms enriched in the PMns using the EnrichR package in the aspect of “biological processes” (Chen et al 2013, Kuleshov et al 2016, Xie et al 2021). GO terms were considered significantly enriched with adjusted p value <0.05.
Raw fastq files and unprocessed aligned data can be accessed for free through the Gene Expression Omnibus (accession number GSE232801). All code used in data processing and figure creation has been deposited, and is freely available on github (https://github.com/JimmyKelly-bio/Single-cell-RNA-seq-analysis-of-spinal-locomotor-circuitry-in-larval-zebrafish).
The authors thank Drs. Apiar Saunders and Alex Nechiporuk for advice with the scRNAseq analysis, Massive Parallel Sequencing Shared Resource at OHSU, for single cell capture, cDNA synthesis and library preparation and sequencing, and Flow Cytometry Shared Resource at OHSU for FACS sorting. Kara Grist provided expert zebrafish husbandry. This research was funded through a grant from the NIH (NS105664) to P.B.
Supplemental Table 1. DEGs in spinal neuron clusters.
Supplemental Table 2. DEGs in PMns versus SMns.
Supplemental Table 3. Changes to the Lawson v4.3.2 reference genome.
- Characterization of developmental and molecular factors underlying release heterogeneity at Drosophila synapsesElife 7https://doi.org/10.7554/eLife.38268
- Interaction of Kv3 Potassium Channels and Resurgent Sodium Current Influences the Rate of Spontaneous Firing of Purkinje NeuronsThe Journal of Neuroscience 26https://doi.org/10.1523/JNEUROSCI.5204-05.2006
- Separate microcircuit modules of distinct v2a interneurons and motoneurons control the speed of locomotionNeuron 83:934–943https://doi.org/10.1016/j.neuron.2014.07.018
- Motoneuron fate specification revealed by patterned LIM homeobox gene expression in embryonic zebrafishDevelopment 121:4117–4125https://doi.org/10.1242/dev.121.12.4117
- An mnr2b/hlxb9lb enhancer trap line that labels spinal and abducens motor neurons in zebrafishDev Dyn 241:327–332https://doi.org/10.1002/dvdy.22781
- Gene ontology: tool for the unification of biology. The Gene Ontology ConsortiumNat Genet 25:25–29https://doi.org/10.1038/75556
- Central and peripheral innervation patterns of defined axial motor units in larval zebrafishJ Comp Neurol 527:2557–2572https://doi.org/10.1002/cne.24689
- Principles Governing Locomotion in Vertebrates: Lessons From ZebrafishFront Neural Circuits 12https://doi.org/10.3389/fncir.2018.00073
- Identification of spinal neurons in the embryonic and larval zebrafishJ Comp Neurol 302:603–616https://doi.org/10.1002/cne.903020315
- Spinal cholinergic interneurons differentially control motoneuron excitability and alter the locomotor network operational rangeSci Rep 8https://doi.org/10.1038/s41598-018-20493-z
- Grading movement strength by changes in firing intensity versus recruitment of spinal interneuronsNeuron 53:91–102https://doi.org/10.1016/j.neuron.2006.11.011
- Single-cell transcriptomic analysis of the adult mouse spinal cord reveals molecular diversity of autonomic and skeletal motor neuronsNat Neurosci 24:572–583https://doi.org/10.1038/s41593-020-00795-0
- Voltage imaging identifies spinal circuits that modulate locomotor adaptation in zebrafishNeuron 110:1211–1222https://doi.org/10.1016/j.neuron.2022.01.001
- Zebrafish information network, the knowledgebase for Danio rerio researchGenetics 220https://doi.org/10.1093/genetics/iyac016
- A small number of open Ca2+ channels trigger transmitter release at a central GABAergic synapseNat Neurosci 13:19–21https://doi.org/10.1038/nn.2461
- Locomotor repertoire of the larval zebrafish: swimming, turning and prey captureJ Exp Biol 203:2565–2579https://doi.org/10.1242/jeb.203.17.2565
- Connexin 35b expression in the spinal cord of Danio rerio embryos and larvaeJ Comp Neurol 522:861–875https://doi.org/10.1002/cne.23449
- Enrichr: interactive and collaborative HTML5 gene list enrichment analysis toolBMC Bioinformatics 14https://doi.org/10.1186/1471-2105-14-128
- Mapping a multiplexed zoo of mRNA expressionDevelopment 143:3632–3637https://doi.org/10.1242/dev.140137
- Comparison and evaluation of statistical error models for scRNA-seqGenome Biol 23https://doi.org/10.1186/s13059-021-02584-9
- Single cell transcriptomics reveals spatial and temporal dynamics of gene expression in the developing mouse spinal cordDevelopment 146https://doi.org/10.1242/dev.173807
- The dual developmental origin of spinal cerebrospinal fluid-contacting neurons gives rise to distinct functional subtypesSci Rep 7https://doi.org/10.1038/s41598-017-00350-1
- Acetylcholinesterase function is dispensable for sensory neurite growth but is critical for neuromuscular synapse stabilityDevelopmental Biology 270:232–245https://doi.org/10.1016/j.ydbio.2004.02.027
- Development of the locomotor network in zebrafishProg Neurobiol 68:85–111https://doi.org/10.1016/s0301-0082(02)00075-8
- Nanodomain coupling between Ca(2)(+) channels and sensors of exocytosis at fast mammalian synapsesNat Rev Neurosci 13:7–21https://doi.org/10.1038/nrn3125
- Pathway selection by growth cones of identified motoneurones in live zebra fish embryosNature 320:269–271https://doi.org/10.1038/320269a0
- Developmental transformation of the release modality at the calyx of Held synapseJ Neurosci 25:4131–4140https://doi.org/10.1523/JNEUROSCI.0350-05.2005
- Molecular pathways of motor neuron injury in amyotrophic lateral sclerosisNat Rev Neurol 7:616–630https://doi.org/10.1038/nrneurol.2011.152
- Spinal network of the Mauthner cellBrain Behav Evol 37:298–316https://doi.org/10.1159/000114367
- Identification of motoneurons and interneurons in the spinal network for escapes initiated by the mauthner cell in goldfishThe Journal of Neuroscience 8https://doi.org/10.1523/JNEUROSCI.08-11-04192.1988
- Zebrafish and motor control over the last decadeBrain Res Rev 57:86–93https://doi.org/10.1016/j.brainresrev.2007.06.018
- Zebrafish as a model system for studying neuronal circuits and behaviorAnn N Y Acad Sci 860:333–345https://doi.org/10.1111/j.1749-6632.1998.tb09060.x
- Inhibition and motor control in the developing zebrafish spinal cordCurr Opin Neurobiol 26:103–109https://doi.org/10.1016/j.conb.2013.12.016
- Neuromuscular synapses can form in vivo by incorporation of initially aneural postsynaptic specializationsDevelopment 132:4471–4481https://doi.org/10.1242/dev.02044
- Real-time visualization of oxidative stress-mediated neurodegeneration of individual spinal motor neurons in vivoRedox Biol 19:226–234https://doi.org/10.1016/j.redox.2018.08.011
- Principles governing recruitment of motoneurons during swimming in zebrafishNat Neurosci 14:93–99https://doi.org/10.1038/nn.2704
- The HMG box transcription factors Sox1a and Sox1b specify a new class of glycinergic interneuron in the spinal cord of zebrafish embryosDevelopment 146https://doi.org/10.1242/dev.172510
- Circuits controlling vertebrate locomotion: moving in a new directionNat Rev Neurosci 10:507–518https://doi.org/10.1038/nrn2608
- Open-channel block by the cytoplasmic tail of sodium channel beta4 as a mechanism for resurgent sodium currentNeuron 45:233–244https://doi.org/10.1016/j.neuron.2004.12.035
- The motor infrastructure: from ion channels to neuronal networksNat Rev Neurosci 4:573–586https://doi.org/10.1038/nrn1137
- Measured motion: searching for simplicity in spinal locomotor networksCurr Opin Neurobiol 19:572–586https://doi.org/10.1016/j.conb.2009.10.011
- Selective degeneration of a physiological subtype of spinal motor neuron in mice with SOD1-linked ALSProc Natl Acad Sci U S A 111:16883–16888https://doi.org/10.1073/pnas.1419497111
- Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regressionGenome Biol 20https://doi.org/10.1186/s13059-019-1874-1
- A confocal study of spinal interneurons in living larval zebrafishJ Comp Neurol 437:1–16https://doi.org/10.1002/cne.1266
- Integrated analysis of multimodal single-cell dataCell 184:3573–3587https://doi.org/10.1016/j.cell.2021.04.048
- Graded Arrays of Spinal and Supraspinal V2a Interneuron Subtypes Underlie Forelimb and Hindlimb Motor ControlNeuron 97:869–884https://doi.org/10.1016/j.neuron.2018.01.023
- Visualization of cranial motor neurons in live transgenic zebrafish expressing green fluorescent protein under the control of the islet-1 promoter/enhancerJ Neurosci 20:206–218https://doi.org/10.1523/JNEUROSCI.20-01-00206.2000
- Distribution of prospective glutamatergic, glycinergic, and GABAergic neurons in embryonic and larval zebrafishJ Comp Neurol 480:1–18https://doi.org/10.1002/cne.20278
- Engrailed-1 expression marks a primitive class of inhibitory spinal interneuronJ Neurosci 24:5827–5839https://doi.org/10.1523/JNEUROSCI.5342-03.2004
- Neurotransmitter properties of spinal interneurons in embryonic and larval zebrafishJ Comp Neurol 480:19–37https://doi.org/10.1002/cne.20279
- Localization of P-type calcium channels in the central nervous systemProceedings of the National Academy of Sciences 88:7076–7080https://doi.org/10.1073/pnas.88.16.7076
- Release probability of hippocampal glutamatergic terminals scales with the size of the active zoneNat Neurosci 15:988–997https://doi.org/10.1038/nn.3137
- Islet1 and Islet2 have equivalent abilities to promote motoneuron formation and to specify motoneuron subtype identityDevelopment 133:2137–2147https://doi.org/10.1242/dev.02355
- Spinocerebellar ataxia type 13 mutant potassium channel alters neuronal excitability and causes locomotor deficits in zebrafishJ Neurosci 31:6831–6841https://doi.org/10.1523/JNEUROSCI.6572-10.2011
- Developmental changes in calcium channel types mediating synaptic transmission in rat auditory brainstemJ Physiol 509:419–423https://doi.org/10.1111/j.1469-7793.1998.419bn.x
- Evx1 and Evx2 specify excitatory neurotransmitter fates and suppress inhibitory fates through a Pax2-independent mechanismNeural Dev 11https://doi.org/10.1186/s13064-016-0059-9
- Central pattern generators deciphered by molecular geneticsNeuron 41:317–321https://doi.org/10.1016/s0896-6273(04)00042-x
- Presynaptic Resurgent Na<sup>+</sup> Currents Sculpt the Action Potential Waveform and Increase Firing Reliability at a CNS Nerve TerminalThe Journal of Neuroscience 30https://doi.org/10.1523/JNEUROSCI.3982-10.2010
- alx, a zebrafish homolog of Chx10, marks ipsilateral descending excitatory interneurons that participate in the regulation of spinal locomotor circuitsJ Neurosci 26:5684–5697https://doi.org/10.1523/JNEUROSCI.4993-05.2006
- Orderly compartmental mapping of premotor inhibition in the developing zebrafish spinal cordScience 370:431–436https://doi.org/10.1126/science.abb4608
- Mapping a sensory-motor network onto a structural and functional ground plan in the hindbrainProc Natl Acad Sci U S A 108:1170–1175https://doi.org/10.1073/pnas.1012189108
- Enrichr: a comprehensive gene set enrichment analysis web server 2016 updateNucleic Acids Res 44:W90–97https://doi.org/10.1093/nar/gkw377
- Canonical Wnt signaling induces BMP-4 to specify slow myofibrogenesis of fetal myoblastsSkeletal Muscle 3https://doi.org/10.1186/2044-5040-3-5
- An improved zebrafish transcriptome annotation for sensitive and comprehensive detection of cell type-specific genesElife 9https://doi.org/10.7554/eLife.55792
- Resurgent current of voltage-gated Na(+) channelsJ Physiol 592:4825–4838https://doi.org/10.1113/jphysiol.2014.277582
- From cells to circuits: development of the zebrafish spinal cordProg Neurobiol 69:419–449https://doi.org/10.1016/s0301-0082(03)00052-2
- Differential gating and recruitment of P/Q-, N-, and R-type Ca2+ channels in hippocampal mossy fiber boutonsJ Neurosci 27:13420–13429https://doi.org/10.1523/JNEUROSCI.1709-07.2007
- Shared versus specialized glycinergic spinal interneurons in axial motor circuits of larval zebrafishJ Neurosci 28:12982–12992https://doi.org/10.1523/JNEUROSCI.3330-08.2008
- Pain: nocistatin spells reliefCurr Biol 8:R525–527https://doi.org/10.1016/s0960-9822(07)00338-7
- UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction
- A topographic map of recruitment in spinal cordNature 446:71–75https://doi.org/10.1038/nature05588
- Spinal interneurons differentiate sequentially from those driving the fastest swimming movements in larval zebrafish to those driving the slowest onesJ Neurosci 29:13566–13577https://doi.org/10.1523/JNEUROSCI.3277-09.2009
- Continuous shifts in the active set of spinal interneurons during changes in locomotor speedNat Neurosci 11:1419–1429https://doi.org/10.1038/nn.2225
- A gradient in endogenous rhythmicity and oscillatory drive matches recruitment order in an axial motor poolJ Neurosci 32:10925–10939https://doi.org/10.1523/JNEUROSCI.1809-12.2012
- Hierarchical control of locomotion by distinct types of spinal V2a interneurons in zebrafishNat Commun 10https://doi.org/10.1038/s41467-019-12240-3
- Differences in the morphology of spinal V2a neurons reflect their recruitment order during swimming in larval zebrafishJ Comp Neurol 522:1232–1248https://doi.org/10.1002/cne.23465
- Promoter analysis in living zebrafish embryos identifies a cis-acting motif required for neuronal expression of GATA-2Proc Natl Acad Sci U S A 94:6267–6272https://doi.org/10.1073/pnas.94.12.6267
- Developmental changes in Ca2+ channel subtypes regulating endocytosis at the calyx of HeldThe Journal of Physiology 592:3495–3510https://doi.org/10.1113/jphysiol.2014.273243
- Functional effects of spinocerebellar ataxia type 13 mutations are conserved in zebrafish Kv3.3 channelsBMC Neurosci 11https://doi.org/10.1186/1471-2202-11-99
- Imaging functional neural circuits in zebrafish with a new GCaMP and the Gal4FF-UAS systemCommun Integr Biol 4:566–568https://doi.org/10.4161/cib.4.5.15848
- Development and axonal outgrowth of identified motoneurons in the zebrafishJ Neurosci 6:2278–2289https://doi.org/10.1523/JNEUROSCI.06-08-02278.1986
- Zebrafish CaV2.1 calcium channels are tailored for fast synchronous neuromuscular transmissionBiophys J 108:578–584https://doi.org/10.1016/j.bpj.2014.11.3484
- Motor neuron vulnerability and resistance in amyotrophic lateral sclerosisActa Neuropathol 133:863–885https://doi.org/10.1007/s00401-017-1708-8
- The Serotonergic Raphe Promote Sleep in Zebrafish and MiceNeuron 103:686–701https://doi.org/10.1016/j.neuron.2019.05.038
- The Zebrafish motility mutant twitch once reveals new roles for rapsyn in synaptic functionJ Neurosci 22:6491–6498
- Transformation of an early-established motor circuit during maturation in zebrafishCell Rep 39https://doi.org/10.1016/j.celrep.2022.110654
- Evoked and spontaneous transmission favored by distinct sets of synapsesCurr Biol 24:484–493https://doi.org/10.1016/j.cub.2014.01.022
- Resurgent sodium current and action potential formation in dissociated cerebellar Purkinje neuronsJ Neurosci 17:4517–4526
- Loss of Navbeta4-Mediated Regulation of Sodium Currents in Adult Purkinje Neurons Disrupts Firing and Impairs Motor Coordination and BalanceCell Rep 19:532–544https://doi.org/10.1016/j.celrep.2017.03.068
- Kv3.3 subunits control presynaptic action potential waveform and neurotransmitter release at a central excitatory synapseElife 11https://doi.org/10.7554/eLife.75219
- Functional role of a specialized class of spinal commissural inhibitory neurons during fast escapes in zebrafishJ Neurosci 29:6780–6793https://doi.org/10.1523/JNEUROSCI.0801-09.2009
- Functional Diversity of Glycinergic Commissural Inhibitory Neurons in Larval ZebrafishCell Rep 30:3036–3050https://doi.org/10.1016/j.celrep.2020.02.015
- Zebrafish Mnx proteins specify one motoneuron subtype and suppress acquisition of interneuron characteristicsNeural Dev 7https://doi.org/10.1186/1749-8104-7-35
- The Cav2.1/alpha1A (P/Q-type) voltage-dependent calcium channel mediates inhibitory neurotransmission onto mouse cerebellar Purkinje cellsEur J Neurosci 13:1902–1912https://doi.org/10.1046/j.0953-816x.2001.01566.x
- Motor neurons and the generation of spinal motor neuron diversityFront Cell Neurosci 8https://doi.org/10.3389/fncel.2014.00293
- Comprehensive Integration of Single-Cell DataCell 177:1888–1902https://doi.org/10.1016/j.cell.2019.05.031
- Swimming of larval zebrafish: fin-axis coordination and implications for function and neural controlJ Exp Biol 207:4175–4183https://doi.org/10.1242/jeb.01285
- From Louvain to Leiden: guaranteeing well-connected communitiesSci Rep 9https://doi.org/10.1038/s41598-019-41695-z
- Visualizing Data using t-SNEJournal of Machine Learning Research 9:2579–2605
- Acetylcholine receptor gating in a zebrafish model for slow-channel syndromeJ Neurosci 32:7941–7948https://doi.org/10.1523/JNEUROSCI.0158-12.2012
- Zebrafish model for congenital myasthenic syndrome reveals mechanisms causal to developmental recoveryProc Natl Acad Sci U S A 109:17711–17716https://doi.org/10.1073/pnas.1215858109
- Function of neuromuscular synapses in the zebrafish choline-acetyltransferase mutant bajanJ Neurophysiol 100:1995–2004https://doi.org/10.1152/jn.90517.2008
- A Gradient in Synaptic Strength and Plasticity among Motoneurons Provides a Peripheral Mechanism for Locomotor ControlCurr Biol 27:415–422https://doi.org/10.1016/j.cub.2016.12.010
- Cryptic organisation within an apparently irregular rostrocaudal distribution of interneurons in the embryonic zebrafish spinal cordExp Cell Res 316:3292–3303https://doi.org/10.1016/j.yexcr.2010.06.020
- Primary and secondary motoneurons use different calcium channel types to control escape and swimming behaviors in zebrafishProc Natl Acad Sci U S A 117:26429–26437https://doi.org/10.1073/pnas.2015866117
- Fatigue in Rapsyn-Deficient Zebrafish Reflects Defective Transmitter ReleaseJ Neurosci 36:10870–10882https://doi.org/10.1523/JNEUROSCI.0505-16.2016
- Zebrafish calls for reinterpretation for the roles of P/Q calcium channels in neuromuscular transmissionJ Neurosci 33:7384–7392https://doi.org/10.1523/JNEUROSCI.5839-12.2013
- Nonequivalent release sites govern synaptic depressionProc Natl Acad Sci U S A 113:E378–386https://doi.org/10.1073/pnas.1523671113
- Identified motoneurons and their innervation of axial muscles in the zebrafishJ Neurosci 6:2267–2277https://doi.org/10.1523/JNEUROSCI.06-08-02267.1986
- Calcium channel types with distinct presynaptic localization couple differentially to transmitter release in single calyx-type synapsesJ Neurosci 19:726–736
- Gene Set Knowledge Discovery with EnrichrCurr Protoc 1https://doi.org/10.1002/cpz1.90
- The Genetic Programs Specifying Kolmer-Agduhr InterneuronsFront Neurosci 14https://doi.org/10.3389/fnins.2020.577879
- Presynaptic voltage-gated calcium channels in the auditory brainstemMol Cell Neurosci 112https://doi.org/10.1016/j.mcn.2021.103609
- A cluster of cholinergic premotor interneurons modulates mouse locomotor activityNeuron 64:645–662https://doi.org/10.1016/j.neuron.2009.10.017
- In vivo labeling of zebrafish motor neurons using an mnx1 enhancer and Gal4/UASGenesis 49:546–554https://doi.org/10.1002/dvg.20766
- Kv3.3 potassium channels and spinocerebellar ataxiaJ Physiol 594:4677–4684https://doi.org/10.1113/JP271343