Identification of adult spinal Shox2 neuronal subpopulations based on unbiased computational clustering of electrophysiological properties

Spinal cord neurons integrate sensory and descending information to produce motor output. The expression of transcription factors has been used to dissect out the neuronal components of circuits underlying behaviors. However, most of the canonical populations of interneurons are heterogeneous and require additional criteria to determine functional subpopulations. Neurons expressing the transcription factor Shox2 can be subclassified based on the co-expression of the transcription factor Chx10 and each subpopulation is proposed to have a distinct connectivity and different role in locomotion. Adult Shox2 neurons have recently been shown to be diverse based on their firing properties. Here, in order to subclassify adult mouse Shox2 neurons, we performed multiple analyses of data collected from whole-cell patch clamp recordings of visually-identified Shox2 neurons from lumbar spinal slices. A smaller set of Chx10 neurons was included in the analyses for validation. We performed k-means and hierarchical unbiased clustering approaches, considering electrophysiological variables. Unlike the categorizations by firing type, the clusters displayed electrophysiological properties that could differentiate between clusters of Shox2 neurons. The presence of clusters consisting exclusively of Shox2 neurons in both clustering techniques suggests that it is possible to distinguish Shox2⁺Chx10⁻ neurons from Shox2⁺Chx10⁺ neurons by electrophysiological properties alone. Computational clusters were further validated by immunohistochemistry with accuracy in a small subset of neurons. Thus, unbiased cluster analysis using electrophysiological properties is a tool that can enhance current interneuronal subclassifications and can complement groupings based on transcription factor and molecular expression.

Introduction

The spinal neuronal circuitry participates in the control of a wide variety of movements, ranging from reflexes to highly sophisticated motor skills (Hultborn and Nielsen, 2007; Goulding, 2009; Todd, 2010; Levine et al., 2012; Kiehn, 2016; Grillner and El Manira, 2020). The identification and connectivity of spinal interneuron populations that underlie specific behaviors is essential for a complete understanding of motor function. Spinal neuronal populations are highly diverse and functional populations are anatomically intermingled (Goulding, 2009; Kiehn, 2016; Abraira et al., 2017; Stachowski and Dougherty, 2021). The identification of spinal interneuronal populations based on the expression of transcription factors has defined neurons involved in or essential to rhythm generation, left-right alternation, flexor-extensor alternation, and patterning (Wilson et al., 2005; Gosgnach et al., 2006; Crone et al., 2008; Zhang et al., 2008, 2014; Dyck et al., 2012; Dougherty et al., 2013; Talpalar et al., 2013; Caldeira et al., 2017). However, each of the transcription factor-defined canonical classes of ventral interneurons identified to date is heterogeneous and serves multiple functions (Pierani et al., 2001; Brownstone and Wilson, 2008; Borowska et al., 2013; Bikoff et al., 2016; Shevtsova and Rybak, 2016; Hayashi et al., 2018; Deska-Gauthier et al., 2020; Falgairolle and O'Donovan, 2021).

As with other genetically identified populations, lumbar spinal cord neurons that express the transcription factor Shox2 are heterogeneous in transcription factor expression, connectivity, and function (Dougherty et al., 2013; Ha and Dougherty, 2018). Shox2 neurons overlap with the V2a (Chx10-expressing) population of neurons (Dougherty et al., 2013). Similarly, V2a neurons are also heterogeneous based on molecular expression (Song et al., 2009; Hayashi et al., 2018), firing types (Dougherty and Kiehn, 2010; Zhong et al., 2010; Husch et al., 2015), activity pattern during locomotor deletions (Zhong et al., 2012), and function (Crone et al., 2008, 2009). There is also interrelation between subpopulations. The Shox2 subpopulation lacking Chx10 expression has been linked to the generation of the locomotor rhythm (Dougherty et al., 2013; Dougherty and Ha, 2019), the Chx10⁺ neuronal population has been implicated in the left-right alternation (Crone et al., 2008, 2009), and the Shox2⁺Chx10⁺ neurons are thought to be pre-motoneurons participating in the stabilization of the locomotor burst (Dougherty et al., 2013).

Particularly for transcription factors that are developmentally downregulated or with deletions that are lethal, the relation of population to function comes largely from studies in neonates (Gosgnach et al., 2006; Crone et al., 2008; Zhang et al., 2008; Dougherty et al., 2013; Talpalar et al., 2013), as do experiments relating cellular firing or phasing during locomotor-like activity (Dougherty and Kiehn, 2010; Zhong et al., 2011; Dyck et al., 2012). Other classification schemes used to subdivide potential functional populations in both neonate and adult include location (Harrison et al., 1986; Puskar and Antal, 1997; Petko et al., 2004; Jankowska, 2008), neurotransmitter phenotype (Huang et al., 2000; Goulding et al., 2014; Bikoff et al., 2016; Hughes and Todd, 2020; Stachowski and Dougherty, 2021), morphology (Gobel, 1975; Maxwell et al., 2007; Yasaka et al., 2010), electrophysiological properties (Konnerth et al., 1990; Lopez-Garcia and King, 1994; Grudt and Perl, 2002; Dai et al., 2009; Dyck et al., 2012) and molecular profiles (Goulding, 2009; Kiehn, 2016; Ziskind-Conhaim and Hochman, 2017; Haring et al., 2018; Gatto et al., 2019). Recent RNA sequence profiling of adult spinal neurons has suggested clusters of spinal interneurons where ventral clusters outnumber the canonical V0–V3 developmentally-derived populations and some clusters contain significant numbers of neurons in more than one canonical population (Sathyamurthy et al., 2018). Similarly, mature neurons may be more diverse than those of neonate electrophysiologically and morphologically (Alvarez et al., 2005; Al-Mosawie et al., 2007). Electrophysiological recordings have demonstrated that adult Shox2 neurons display heterogeneous electrophysiological properties (Garcia-Ramirez et al., 2021), suggesting the possibility of more than two subpopulations. Electrophysiological properties establish the ways in which neurons respond to the synaptic inputs and neuromodulation, determine integration capacity, and shape output of individual neurons; therefore, along with connectivity, they govern spinal circuit function (Konnerth et al., 1990; Russo and Hounsgaard, 1999; Lee and Heckman, 2000; Butt et al., 2002; Miles et al., 2005; Smith and Perrier, 2006; Brownstone and Wilson, 2008; Dai et al., 2009; Grillner and El Manira, 2020). Thus, classification of neurons based on their passive and active properties can provide essential characteristics that should correlate with their participation in behavior.

Here, rather than constraining subpopulations of Shox2 neurons based on the expression of a second transcription factor or other molecular marker, we use electrophysiological properties to define subpopulations and then determine how electrophysiolocially-defined subpopulations relate to transcription factor-defined subdivisions. Using neuronal firing type and two clustering approaches, k-means clustering and hierarchical clustering, we identify the electrophysiological parameters that are characteristic to each subpopulation of Shox2 neurons. We found 4 populations of adult Shox2 neurons based on firing type, all of which were in common with Chx10 neurons. Analysis performed on 4 k-means clusters and 6 hierarchical clusters of Shox2 neurons identified clusters which overlap with Chx10 neurons to varying degrees. In both computational analyses, we found clusters that included Shox2 neurons but not Chx10 neurons, suggesting the possibility to identify Shox2⁺Chx10⁻ neurons with a single lineage traced mouse line. Unbiased computational analysis is a useful tool to identify populations of neurons in adult mice where electrophysiological diversity is high and extends beyond current molecular markers and transgenic tools. Further, results could have important implications for the study of neurons in intact preparations where the ability to visually identify neurons is limited.

Materials and methods

Experiments were performed using 51 Shox2::Cre; Rosa26-flox-Stop-flox-tdTomato (Ai9 from Jax mice, #007909) transgenic mice (Madisen et al., 2010; Dougherty et al., 2013) and 6 Chx10-eGFP (from MMRRC, Cat#:011391-UCD) transgenic mice. Electrophysiological properties from 60 of the Shox2 neurons (from 28 of the mice) included in the new analyses were reported in a previous study (Garcia-Ramirez et al., 2021). All experimental procedures followed National Institutes of Health guidelines and were approved by the Institutional Animal Care and Use Committee at Drexel University.

Spinal cord preparations

For terminal electrophysiology experiments, adult mice (8–25 weeks postnatal) were anesthetized with ketamine (150 mg/kg) and xylazine (15 mg/kg), decapitated, and eviscerated. Spinal cords were then removed in ice-cold dissecting solution containing the following (in mM): 222 glycerol, 3 KCl, 11 glucose, 25 NaHCO₃, 1.3 MgSO₄, 1.1 KH₂PO₄, and 2.5 CaCl₂ and gassed with 95% O₂ and 5% CO₂ with pH adjusted to 7.4. The lumbar spinal cord (segments L2–5) was sectioned transversely using a vibrating microtome (Leica Microsystems). Slices were next transferred to an artificial cerebrospinal fluid, containing the following (in mM): 111 NaCl, 3 KCl, 11 glucose, 25 NaHCO₃, 1.3 MgSO₄, 1.1 KH₂PO₄, and 2.5 CaCl₂ at 37°C for 30 min and then passively equilibrated to room temperature for another 30 min before recording. Dissecting and recording solutions were continuously aerated with 95%/5% O₂/CO₂.

Whole-cell patch clamp recordings

All recordings were performed at room temperature. Fluorescently labeled (tdTomato) Shox2 and (eGFP) Chx10 neurons were visualized with a 63x objective lens on a BX51WI scope (Olympus) using LED illumination (Andor Mosaic System or Lumen Dynamics X-Cite 120 LED). Patch electrodes were pulled to tip resistances of 5–8 MΩ using a multistage puller (Sutter Instruments) and were filled with intracellular solution, which contained the following (in mM): 128 K-gluconate, 10 HEPES, 0.0001 CaCl₂, 1 glucose, 4 NaCl, 5 ATP, and 0.3 GTP, with pH adjusted to 7.4. In some experiments, biocytin (2 mg/ml, B4261, Sigma-Aldrich) was included in the patch electrode. Data were collected with a Multiclamp 700B amplifier (Molecular Devices) and Clampex software (pClamp9, Molecular Devices). Signals were digitized at 20 kHz and filtered at 10 kHz.

Measurement and calculation of passive and active membrane properties

We performed whole-cell patch clamp recordings from Shox2 neurons and Chx10 neurons identified by fluorescence expression. Following exclusion of neurons with a resting membrane potential more depolarized than −40 mV and neurons with an action potential peak that did not reach 0 mV, 143 Shox2 neurons and 28 Chx10 neurons were considered for this study. Passive and active membrane properties were recorded or calculated as described briefly here. Current step protocols were run when neurons were at ~-70 mV and in most cases this required bias current. Resting membrane potential (E_m) was recorded shortly after entering whole-cell mode. Input resistance (R_in) was calculated from the current/voltage slope in response to hyperpolarizing voltage steps in which no voltage-gated current activation was evident. Membrane time constant (τ) was calculated as the time to reach 63% of the maximum depolarization in response to a subthreshold depolarizing current step. Membrane capacitance (C_m) was calculated from the time constant and input resistance (C_m = R_in/τ). Rheobase was the minimal current step required to generate an action potential, applied in intervals of 2–3pA. Additionally, we recorded action potential (AP) properties from the first AP observed in the depolarizing current at rheobase (Figure 1A). AP voltage threshold was the membrane potential at the first deflection of the AP. AP half width (AP duration) was determined at the voltage midway from the AP threshold to the peak. The fast afterhyperpolarization (fAHP) amplitude was the difference between the AP threshold and the maximum hyperpolarization. We normalized the fAHP amplitude by dividing by the AP amplitude. The durations of the fAHP and slow afterhyperpolarization (sAHP) were measured by the time elapsed from the AP threshold to the peak fast or slow hyperpolarization. For neurons firing with an initial double or initial burst, the slow and fast AHP values were measured from the AP threshold of the last AP of the double or burst. AP frequency was calculated as the number of action potentials in response to a 1 second depolarizing current injection at 1.5x rheobase. Frequency-current (F/I) curve slope was obtained from the number of APs generated by the administration of depolarizing current pulses from 1 to 49 pA at intervals of 3pA (Figure 1B). The AP frequency at resting membrane potential was calculated from the number of APs recorded in 1 min. We also recorded the voltage onset of the persistent inward currents (PIC_on voltage, Figure 1C), measured as the membrane potential at which a negative deflection began during the application of a slow (28 mV/s) voltage ramp.

FIGURE 1

Figure 1. Measured electrophysiological properties. (A) Representative example of a recording of the membrane potential of a Shox2 neuron during the injection of a 1 s long current step at rheobase (bottom) and magnification of the first action potential generated (top). Arrows point out the voltage threshold and fast and slow AHP. The AP and fAHP amplitude and duration measurement points are shown. (B) Action potential firing frequency in response to the injection of depolarizing currents in two representative Shox2 neurons (squares, cell #40; circles, cell #160). Lines represent the F/I fit slope (arrow). (C) Current response (upper trace) to slow depolarizing ramp injection. Arrows point out the presence of PIC and the corresponding voltage on the ramp to activate the PIC (PIC_on).

Statistical analyses

Data analysis was performed with Clampfit (Molecular Devices) and MATLAB (MathWorks). Statistical tests and post-hoc analyses used are stated for each experiment and performed with GraphPad Prism (GraphPad Software) and MATLAB. All results are presented as mean ± SD. Statistical significance was set at p < 0.05 unless otherwise stated. The distribution of the data was determined by Shapiro–Wilk normality test. The statistical comparisons between Shox2 and Chx10 neurons were performed by Mann-Whitney test or unpaired t-test. Comparisons between clusters were performed by Kruskal-Wallis with Dunn's post-hoc test or repeated-measures one-way ANOVA with Bonferroni post-hoc test. Comparisons between percentages were performed by chi-square test.

Correlations, principal component analysis, and cluster analysis

To determine correlations in between the 12 properties obtained for each cell, we performed a Pearson's linear correlation test and considered values of p < 0.01 to be correlated. Principal component analysis (PCA) and multidimensional cluster analyses were performed on the 6 parameters obtained that were not highly correlated. PCA was used to reduce dimensionality of the number of variables recorded (Figure 2). To visualize clustering results, we show 3D graphs displaying neurons or means of clusters in a 3-dimensional space formed by the 3 first PCAs that explained more than 63% of the variability of the data. We classified Shox2 and Chx10 neurons by 3 different methods: (1) Firing types were determined by responses to depolarizing current steps. (2) For k-means clustering, the number of clusters was determined by the elbow method on a silhouette value graph after iterating the algorithm 1,000 times considering k = 2 to k = 8 clusters. Details can be found in the results section. (3) Hierarchical clusters were determined by the cosine distance between pairs of observations and consideration of an average for the distance between clusters. This combination was used because it resulted in the highest cophenetic correlation coefficient. The number of clusters was determined considering the cutoff below the maximum inconsistency coefficient for each link. For the k-means and hierarchical cluster analyses, we applied MATLAB algorithms, initially standardizing the data to set a mean of 0 and standard deviation of 1 to be able to compare variables with different units.

FIGURE 2

Figure 2. Reducing dimensionality of electrophysiological properties. (A) Correlation matrix of p-values from a Pearson's linear correlation of 12 (Ai) and 6 (Aii) electrophysiological properties from 143 Shox2 and 28 Chx10 neurons. Dark shading corresponds to highly correlated values (p < 0.001). (B) Three-dimensional representation of all 171 neurons included (black dots) after performing principal component analysis (PCA). X axis is the 1^st principal component (PC), Y axis is the 2^nd PC, and Z axis is the 3^rd PC, representing 23.0%, 21.7% and 18.3% of the variability, respectively.

Immunostaining and visualization of biocytin-filled neurons

The slices containing biocytin-filled neurons were fixed overnight (4% paraformaldehyde in 0.1 M PBS; pH 7.4; 4°C), and subsequently placed in PBS at 4°C. To visualize biocytin, the slices were incubated for 2 h at room temperature in DyLight 633 conjugated streptavidin (21844, Thermo Fischer Scientific). Then, the slices were washed in PBS (3 x 10 min), permeated with 1% Bovine Serum Albumin (BSA), 5% Donkey Serum (NGS), 0.1% Fish Gelatin, and 0.2% Triton x-100 for 1 h and incubated in sheep anti-Chx10 antibody 1:100 (AB9016, Chemicon) at room temperature for 48 hrs. Slices were washed in PBS (3 x 10 min) and incubated for 2 h at room temperature with rabbit anti-sheep Dylight 488 1:400 (SA5-10054, Thermo Fischer Scientific) followed by a final wash in PBS (3 x 10 min). Slices were then placed on slides within tissue spacers and coverslipped with a Vectashield mounting medium (Vector labs). Sections were scanned using a laser scanning confocal microscope (Leica True Confocal System SP8) in stacks of 10 optical sections across approximately 50 μm at 20x magnification. Images were condensed into maximum projections using the Leica collection software and brightness and contrast were adjusted in ImageJ.

Results

Correlation of the electrophysiological properties of adult Shox2 and Chx10 neurons

Excitatory spinal neurons expressing the transcription factor Shox2 are heterogeneous in terms of firing properties (Garcia-Ramirez et al., 2021). Since the passive and active properties shape the recruitment of, integration of inputs to, and output of neurons (Wilson et al., 2007; Tazerart et al., 2008; Harris-Warrick, 2010; Dougherty and Ha, 2019; Grillner and El Manira, 2020), we hypothesized that these properties could be used to subclassify Shox2 neurons. Since the Chx10 neuronal population partly overlaps with Shox2 neurons and is also heterogeneous (Dougherty and Kiehn, 2010; Zhong et al., 2010, 2012; Dougherty et al., 2013; Husch et al., 2015) we included a smaller set of Chx10 neurons in our sample. We first considered passive and active membrane properties of Shox2 and Chx10 neurons (Figure 1) and performed a correlation analysis (Figure 2Ai) to eliminate redundancy in the parameters measured. We found high correlations (p < 0.001) between some passive and active cellular properties (Figure 2Ai black squares). We removed the variables that were highly correlated from the analysis, which left 6 variables that were not highly correlated: resting membrane potential (E_m), capacitance (C_m), fast afterhyperpolarization (fAHP) duration, fAHP amplitude, activation voltage of the persistent inward current (PIC_on) and frequency-current (F/I) slope (Figure 2Aii). To visualize the neurons in a three-dimensional plane, we performed a principal component analysis (PCA) based on the six variables (Figure 2B). The first three principal components explained 63% of the variability of the 171 neurons. This analysis shows that dimensionality of the variables recorded can be reduced by removing highly correlated properties and considering only 6 neuronal properties.

Electrophysiological comparison between Shox2 and Chx10 neuronal populations

Graphical display of the three first PCAs shows that Shox2 and Chx10 neurons partially overlap. We compared the initial 12 electrophysiological properties of adult Shox2 neurons (n = 143) with those of Chx10 neurons (n = 28) (Figure 3, Supplementary Table 1). Shox2 neurons had lower input resistance (R_in), larger capacitance (C_m), shorter AP half width, shorter fAHP duration, and more depolarized PIC activation voltages (PIC on) than Chx10 neurons. Resting membrane potential, time constant, rheobase, action potential threshold, the sAHP duration, fAHP amplitude, and F/I slope were similar between the groups. This comparison shows that Shox2 neurons and Chx10 neurons share some electrophysiological characteristics, however, differences are observed between subpopulations.

FIGURE 3

Figure 3. Comparison between Shox2 and Chx10 neurons. (A) Three-dimensional representation of 143 Shox2 neurons (black dots) and 28 Chx10 neurons (red dots) on the three principal components (same as Figure 2, but different angle). Ellipsoids are centered on the mean for the three first PCs and semi-axes are generated by the standard deviation of each of the three first PCs. Black and red ellipsoids represent groups of Shox2 and Chx10 neurons, respectively. (B) Comparison of electrophysiological properties measured from Shox2 neurons (black) and Chx10 neurons (red). Filled circles represent the mean and error bars the standard deviation, unfilled gray circles represent individual neurons. Resting membrane potential (E_m), input resistance (R_in), time constant (τ), calculated capacitance (C_m), rheobase, action potential (AP) threshold, AP width, fast after hyperpolarization (fAHP) duration, slow afterhyperpolarization (sAHP) duration, fast after hyperpolarization (fAHP) amplitude, PIC activation voltage (PIC_on), and slope of the AP frequency/current (F/I) curve. *p < 0.05, unpaired t-test or Mann-Whitney test.

Classification of the neurons based on firing type

To identify differences in active and passive cellular properties based on firing type, we classified the 171 Shox2 and Chx10 neurons based on the response to suprathreshold depolarizing current steps. We found four types of responses (Figure 4A). Neurons firing throughout the current step were most prevalent. Over half of the neurons fired action potentials throughout and are referred to here as tonic firing neurons (Figure 4 cyan, n = 88, 51.5%). Initial doublet neurons (blue, n = 49, 28.7%) fire a doublet of action potentials at the start of the step (initial interspike interval <40 ms) and continue firing throughout the step but at a lower frequency (steady frequency of 3.9 ± 5.1 Hz). Neurons with burst of action potentials at the start of the current step are called initial burst firing neurons (red, n = 24, 14%). These neurons displayed three or more initial spikes at high frequency (>25 Hz) and were either silent or fired action potentials later in the step but these action potentials were not a regular frequency like the tonic neurons. Lastly, a small number of neurons, delay neurons (green, n = 10, 5.9%) fired action potentials after a delay from the beginning of the current step. We next looked at the prevalence of firing types in Shox2 and Chx10 neuronal populations. Of the 143 Shox2 neurons, tonic firing neurons were most common (n = 78, 54.6%), followed by initial doublet (n = 36, 25.2%) and initial burst (n = 21, 14.7%) firing, with delayed firing neurons (n = 8, 5.6%) being rare (Figure 4Bii, left side). Of the 28 Chx10 neurons, initial doublet firing neurons were most prevalent (n = 13, 46.4%), over one third were tonic firing (n = 10, 35.7%), and few neurons fired with an initial burst (n = 3, 10.7%) or were delayed firing (n = 2, 7.1%, Figure 4Bii, right side). The difference in the distribution of the 4 firing types in Shox2 and Chx10 neuronal populations was statistically significant (${\text{χ}}_{(3)}^2$ = 30.71, p < 0.0001). Chx10 neurons preferentially displayed initial doublet firing and Shox2 neurons were mostly tonic firing.

FIGURE 4

Figure 4. Shox2 and Chx10 neurons sorted based on type of firing. (A) Representative traces of the four action potential firing types in response to the injection of a 1-s depolarizing current step. (Bi) Pie chart and incidence of the four types of firing responses observed in all 171 neurons. 83.6% of the neurons were Shox2 neurons (solid) and 16.4% were Chx10 neurons (hatched). (Bii) Incidence of the four firing types observed in 143 Shox2 neurons (left) and in 28 Chx10 neurons (right). (C) Three-dimensional representation of all neurons (Shox2 filled dots, Chx10 unfilled dots) on the three principal component axes (same as Figure 2B, but different angle displayed). Ellipsoids are centered on the mean for the three first PCs and semi-axes are generated by the standard deviation of each of the three first PCs. Cyan, blue, red, and green ellipsoids represent tonic, initial doublet, initial burst, and delay firing subgroups, respectively. (D) Comparisons of electrophysiological properties of tonic, initial doublet, initial burst and delay firing neurons. Colors as in other panels. Filled circles represent the mean and error bars the standard deviation, unfilled gray circles represent individual neurons. Properties as in Figure 3. *p < 0.05, one-way ANOVA with Tukey post-hoc test or Kruskal-Wallis with Dunn's post-hoc test. * significantly different from the corresponding color coded group.

To identify the electrophysiological characteristics that corresponded to the different firing types, we performed statistical analyses on both groups considered as a single population (Shox2 and Chx10 neurons, Figure 4D and Supplementary Table 2). Neurons with delayed firing had lower input resistances than initial doublet firing neurons, longer action potential durations than tonic neurons, shorter slow AHP durations than burst firing neurons, and more depolarized PIC activations than tonic and initial doublet neurons. Initial burst firing neurons had longer slow AHPs and more hyperpolarized PIC activation thresholds than both delay and tonic firing neurons. Initial doublet neurons had longer slow AHPs than tonic firing neurons. Thus, there was no property that distinguished one firing type from all others but there were differences in properties, particularly those related to action potential duration and AHP.

We also compared properties of Shox2 and Chx10 neurons within each firing group (Table 1). Considering only tonic firing neurons, Shox2 neurons had longer time constants (Mann-Whitney test, U = 236, p = 0.04), higher capacitance (Mann-Whitney test, U = 199, p = 0.01), shorter AP half width (Mann-Whitney test, U = 195.5, p = 0.009), and shorter fast (Mann-Whitney test, U = 181, p = 0.005) but longer slow (Mann-Whitney test, U = 222, p = 0.02) AHPs. We did not find differences between Shox2 and Chx10 neurons with initial doublet firing patterns. The numbers of delayed and initial burst firing Chx10 neurons were low and precluded from statistical analysis. The electrophysiological properties between neurons classified by type of firing were expected to be different since action potential features and spike frequency rely on these properties. Although firing types were differentially distributed in Shox2 and Chx10 populations, it is not possible to separate Shox2 and Chx10 neuronal populations by firing type since the types are common to both.

TABLE 1

Table 1. Comparison of Shox2 and Chx10 neurons within each firing type.

k-means cluster analysis to classify Shox2 neuronal subtypes

Firing type classification is based on a subjective analysis which could bias the results. To classify Shox2 neurons based on active and passive cellular properties in an unbiased way, we performed a k-means cluster analysis (see Methods) on the set of 171 neurons, which included mainly Shox2 neurons and a small sampling of Chx10 neurons to serve as comparison with a known subdivision with the Shox2 population. The cluster analysis considered the six electrophysiological parameters that were not highly correlated, as described above (Figure 2A). To select the number of clusters (k) that best fit the data, we ran the k-means algorithm one thousand times from k = 2–8 and plotted the highest silhouette value in each condition (Figure 5A). We found that that the silhouette value reached a plateau from k = 4 to k = 8 and, therefore, we selected k=4 for the analysis.

FIGURE 5

Figure 5. k-means clustering. (A) Plot of the highest silhouette values for each number of clusters (k, x axis). (Bi) Incidence of the four k clusters in 171 neurons, with 83.6% of the neurons from the Shox2 population (solid) and 16.4% were Chx10 neurons (hatched). k1 cluster (purple, n = 23), k2 cluster (cyan, n = 11), k3 cluster (yellow, n = 68), and k4 cluster (gray, n = 69). The same color code is used throughout the figure. (Bii) Incidence of k clusters observed in 143 Shox2 neurons (left): k1 cluster (n = 23), k2 cluster (n = 7), k3 cluster (n = 56) and k4 cluster (n = 57). Incidence of the k clusters observed in 28 Chx10 neurons (right): k1 (n = 0), k2 cluster (n = 4), k3 cluster (n = 12), and k4 cluster (n = 12). (C) Plot of the cluster silhouette values to demonstrate how similar the individual neurons are to the neurons of their own clusters. Each bar represents a single neuron. (D) Three-dimensional representation of all neurons (Shox2 filled dots, Chx10 unfilled dots) on the three principal components (same as Figure 2B but different angle displayed). Ellipsoids are centered on the mean for the three first PC and semi-axes are generated by the standard deviation of each of the three first PCs. (E) Comparisons of electrophysiological properties of the neurons from k1, k2, k3 and k4 clusters. Filled circle represents the mean and error bars the standard deviation, unfilled gray circles represent individual neurons. *p < 0.05, one-way ANOVA with Tukey post-hoc test or Kruskal-Wallis with Dunn's post-hoc test. * significantly different from the corresponding color coded group.

The k-mean cluster algorithm generated four clusters (Figure 5Bi). Two of the clusters contained the majority of the neurons (k3 and k4), with each containing ~40% of the population. Since Shox2 neurons were 83% of the neurons included in the analysis, the distribution of the 143 Shox2 neurons was largely similar to that of the entire population (Figure 5Bii, left side) with 23, 7, 56, and 57 neurons in clusters were k1, k2, k3, and k4, respectively. Of the 28 Chx10 neurons (Figure 5Bii, right side), no neurons were in k1, 4 were in k2, and k3 and k4 each included 12 Chx10 neurons. The difference in the distribution of the 4 k-clusters was statistically significant comparing Shox2 and Chx10 populations (${\text{χ}}_{(3)}^2$ = 22.85, p<0.0001). k1 cluster is composed exclusively of Shox2 neurons (binomial test p = 0.01). Chx10 neurons appear more prevalent in k2 cluster than Shox2 neurons, with 14% of Chx10 neurons but only <5% of Shox2 neurons belonging to that cluster, but this is not significant (binomial test p = 0.09). The proportion of Chx10 neurons in each of the k3 and k4 clusters (39.2% and 39.9%) is similar to that of Shox2 neurons (42.9% each, binomial test, p = 0.7 and p = 0.8, for k3 and k4, respectively). The silhouette values for each of the 171 neurons in the k-clusters (Figure 5C) show that most of the neurons have positive values, indicating that the neuron is close to its cluster centroid and far from the centroids of the other clusters. Few neurons (n = 8) have negative values which indicate that the neuron is close to its cluster centroid but also close to other cluster centroids. Together, this shows that the k-means analysis (Figure 5D) associates the majority of neurons suitably and that k3 and k4 are comprised of Shox2 and Chx10 neurons rather equally and k1 is exclusively constituted of Shox2 neurons.

To identify the electrophysiological characteristics of each k-cluster, we performed statistical analysis on the 4 k-clusters considering the 171 neurons (Figure 5, Supplementary Table 3). Unlike when separated by firing type, here, there were more differences between clusters in both the 6 properties used for the analysis and the 6 other properties that were found to be highly correlated (Figure 5E, compared to Figure 4D). Neurons in the k1 cluster are distinguished from all other clusters in that they have the lowest input resistance and highest capacitance, suggesting that these are larger neurons. They also had the longest time constant and most depolarized voltage thresholds, although these were not significantly different from the neurons in the k2 cluster. Neurons in the k2 cluster can be distinguished from the other clusters by their longer AP and fast AHP durations. Neurons in the k3 cluster have the largest fast AHPs and have the most depolarized membrane potentials, although the latter is not different than neurons in the k2 cluster. Neurons in the k4 cluster have the most hyperpolarized resting membrane potentials and narrow action potential, although neither of these are significantly different from k1 neurons. Neurons in the k4 cluster have more hyperpolarized PIC activation voltages than neurons in k1. Taken together, this suggests that k1 cluster properties may be characteristic of Shox2⁺Chx10⁻ neurons while properties of k3 and k4 clusters, and k2 cluster to a lesser extent, may be characteristic of Shox2⁺Chx10⁺ neurons.

We next performed statistical analysis of the k clusters considering Shox2 and Chx10 neurons separately (Table 2). In k2, which had a higher proportion of Chx10 neurons and few Shox2 neurons, the only significant difference was that C_m of k2 Shox2 neurons was higher than of k2 Chx10 neurons (Mann-Whitney test, U = 2, p = 0.02). We did not find differences in k3 cluster (with equal proportions of Shox2 and Chx10 neurons) between Shox2 and Chx10 populations. The k4 also had equal proportions of Shox2 and Chx10 neurons and the only differences between the Shox2 and Chx10 neurons in that cluster were in AP halfwidth (Mann-Whitney test, U = 168.5, p = 0.005) and fast AHP amplitude (unpaired t test, t_(2.04) = 67, p = 0.04). Notably, there were less differences between Shox2 and Chx10 neurons from individual k clusters than there were differences in the whole populations of Shox2 and Chx10 neurons. This demonstrates similarity within clusters regardless of the transcription factor expressed.

TABLE 2

Table 2. Comparison of properties of Shox2 and Chx10 neurons in each k cluster.

Hierarchical cluster analysis to classify subpopulations of Shox2 and Chx10 neurons

The number of clusters in the k-means analysis here was defined by the elbow rule on a graph of the silhouette values (Allen et al., 2014; Vergara et al., 2020) which can be considered arbitrary. Hierarchical clustering produces a dendrogram that facilitates the visualization of natural clustering, and, therefore, represents a more unbiased method of classification (Martinez et al., 2017; Di Miceli et al., 2020). Thus, we performed an unbiased hierarchical clustering (see Methods) considering the same set of data used for k-means algorithm (171 neurons, 6 variables in Figure 2Aii). The algorithm established 6 clusters generated by natural divisions with cutoffs below the value of the maximum inconsistency coefficient. The component neurons of each hierarchical cluster (H cluster) can be visualized in a dendrogram (Figure 6A).

FIGURE 6

Figure 6. Hierarchical clustering. (A) Dendrogram plot of the hierarchical binary cluster tree of 171 neurons. Clusters are generated by natural divisions with the cutoff below the value of the maximum inconsistency coefficient. Clusters are color coded: blue for H1 cluster, orange for H2 cluster, yellow for H3 cluster, purple for H4 cluster, green for H5 cluster, and cyan for H6 cluster. The same color code is used for the other panels. The lines below correspond to the type of neurons with Shox2⁺ neurons in black and Chx10⁺ neurons in red. (Bi) Incidence of the six H clusters in all 171 Shox2 neurons (solid) and Chx10 neurons (hatched). H1 cluster (n = 15), H2 cluster (n = 20), H3 cluster (n = 32), H4 cluster (n = 25), H5 cluster (n = 27) and H6 cluster (n = 52). (Bii) Incidence of the H clusters observed in 143 Shox2 neurons (left): H1 (n = 15), H2 (n = 20), H3 (n = 24), H4 (n = 20), H5 (n = 25) and H6 cluster (n = 39). Incidence of the H clusters observed in 28 Chx10 neurons (right): H1 and H2 (n = 0), H3 (n = 8), H4 (n = 5), H5 (n = 2) and H6 cluster (n = 13). (C) Three-dimensional representation of all neurons (Shox2 filled, Chx10 open) on the three principal components (same as Figure 2B but a different angle). Ellipsoids are centered on the mean for the three first PCs and semi-axes are generated by the standard deviation of each three first PC. (D) Plot of the cluster silhouettes values. (E) Comparisons of electrophysiological properties of the neurons from H1, H2, H3, H4, H5 and H6 clusters. Filled circles and error bars represent the mean and standard deviation, open gray circles represent individual neurons. *p < 0.05, one-way ANOVA with Tukey post-hoc test or Kruskal-Wallis with Dunn's post-hoc test. *significantly different from the corresponding color coded group.

Each of the six H clusters (Figure 6Bi) contained 8.8–30.4% of the input population. Based on the silhouette values and the ellipsoids on the PCA graphs (Figures 6C,D), the first 5 clusters described their constituent neurons well with mainly positive values, where the largest cluster was H6 and 29% of neurons had negative values. When considering the 143 Shox2 neurons (Figure 6Bii, left) separately, neurons were relatively evenly distributed between H1–H5 (n = 15–25 in each) with more neurons in H6 (n = 39). The 28 Chx10 neurons (Figure 6Bii, right) were distributed between H3–H6 clusters (n = 2–13) but there are no Chx10 neurons in H1 and H2 clusters. The difference in the distribution of the 6 H-clusters was statistically significant when comparing Shox2 and Chx10 populations (${\text{χ}}_{(5)}^2$ = 14.95, p = 0.01). Note that H1 and H2 clusters were composed of Shox2 neurons exclusively and these two clusters were most similar to each other, as seen by the height of the next branch point in the analysis. Although clusters H3-H6 contained both Shox2 and Chx10 neurons, the uneven numbers of Shox2 and Chx10 input neurons (143 vs. 28) should be considered (Figure 6Bii). H3 and H6 clusters were enriched in Chx10 neurons, H5 was enriched in Shox2 neurons, and the proportions of Shox2 and Chx10 neurons in H4 were similar.

To identify the electrophysiological characteristics of each H-cluster, we performed statistical analysis on these clusters considering all 171 neurons (Figure 6E, Supplementary Table 4). Neurons in the H1 cluster displayed the smallest fast AHP amplitudes and highest firing frequencies in response to depolarizing current steps but these were not significantly different from all other clusters. The neurons in the H2 cluster had the lowest input resistance (significantly different from all but H5), highest capacitance, and most depolarized AP threshold (significantly different from all but H1 and H5). Neurons in the H3 cluster had the most hyperpolarized resting membrane potentials, AP thresholds more hyperpolarized than neurons in H2 and H5, and the most hyperpolarized PIC activations (significantly different from all but H4). Neurons in the H4 cluster had the highest input resistances, lowest capacitance values, and longest fast AHP durations, although none of these properties were significantly different from neurons in the H6 cluster. Neurons in the H5 cluster had the highest fast AHP amplitude (significantly different from all but H4). Lastly, neurons in the H6 cluster had the lowest firing frequencies in response to depolarizing current steps (F-I slope, significantly different from all but H3). Thus, hierarchical clustering generated six distinct clusters, two of which were composed only of Shox2 neurons, suggesting that larger neurons (low R_in, high C_m) with depolarized voltage thresholds are characteristic of Shox2⁺Chx10⁻ neurons.

We also performed statistical analysis of the H clusters considering Shox2 and Chx10 neurons separately (Table 3). There were no Chx10 neurons in H1 and H2 clusters for analysis and H5 contained only 2 Chx10 neurons which precluded statistical testing. There were no differences between Shox2 and Chx10 neurons in the H6 cluster. There were only a few differences between Shox2 and Chx10 neurons in H3 and H4 clusters. Shox2 neurons in H3 had more depolarized AP thresholds and shorter AP half widths than Chx10 neurons in H3. In H4, the only significant difference was that Shox2 neurons had higher capacitance than the Chx10 neurons. Similar to the k-means clusters, there were few differences between Shox2 and Chx10 neurons within any of the hierarchical clusters. These were outnumbered by differences between the clusters, suggesting relatively homogeneous populations within each cluster.

TABLE 3

Table 3. Comparison of properties of Shox2 and Chx10 neurons in each of the hierarchical clusters.

Comparison of k-mean and hierarchical clusters

In order to validate clustering results and determine the differences between them, we paired the clusters generated by k-means and hierarchical algorithms (Supplementary Figure 1). k-means and H clusters overlapped in 13 combinations (Figure 7A, left) and showed high correspondence. Considering the H clusters with respect to the k clusters, all 15 of the neurons in H1 corresponded to k3, all (32) in H3 were in k4, and nearly all (19 of 20) of the neurons in H2 were in k1. Further, the majority of H5 and H6 neurons were in k4 and k3, respectively, where neurons in H4 were distributed but absent from k1. The correspondence can also be considered the other way around. Most neurons in the k1 cluster (19 of 23) are in H2, most neurons in k2 (9 of 11) correspond to H4 neurons, neurons belonging to k3 were mainly distributed on H1 and H6, and neurons in k4 neurons were dispersed but mainly in H3 and H5 (Figure 7B). Furthermore, there were 6 combinations with high overlap (>80% of maximum possibilities of pair appearance, Figure 7A, right). These 6 combinations accounted for 86% of the 171 input neurons (143 Shox2 neurons and 25 Chx10 neurons).

FIGURE 7

Figure 7. Comparison between k and H clusters. (A) Thirteen possible combinations between k-cluster and H-clusters (k & H pairing, left) with the corresponding maximum percentage of appearance (right). The pairings with percentages higher than 80% are marked with red arrows. k-clusters (left boxes), k1 (purple), k2 (blue), k3 (yellow), k4 (gray), H-clusters (right boxes), H1 (blue), H2 (orange), H3 (yellow), H4 (purple), H5 (green) and H6 (cyan). (B) Table shows the k & H pairing. Numbers indicate the neurons within the corresponding cluster, k in columns and H in rows. Red numbers are the neurons with pairing with percentages higher than 80%.

Further validation of clustering with immunohistochemistry identified cells

A subset of the Shox2 neurons used for this analysis were recorded with biocytin in the electrode and were recovered for labeling with an antibody to Chx10. In total, we stained 8 neurons to determine Shox2⁺Chx10⁺ or Shox2⁺Chx10⁻ identity (Figure 8A). Of those stained, 2 Shox2 neurons were Chx10⁻ (Shox2⁺Chx10⁻) and 6 Shox2 neurons expressed Chx10 (Shox2⁺Chx10⁺). We next determined which clusters each of the 8 neurons belonged to. Considering k clusters (Figure 8Bi), we found that Shox2⁺Chx10⁻ were classified in k1 (1 neuron) and k3 (1 neuron) clusters, while Shox2⁺Chx10⁺ neurons were classified in k2 (1 neuron), k3 (2 neurons), and k4 (1 neuron) clusters. This matches rather well with predictions from the clustering results (Figure 5) because k1 contained only Shox2 neurons and k2–k4 contained both Shox2 and Chx10 neurons. Considering H clusters (Figure 8Bii), we found that Shox2⁺Chx10⁻ where classified in H1 (1 neuron) and H2 (1 neuron), the clusters devoid of Chx10 neurons. The Shox2⁺Chx10⁺ neurons were in H3 (1 neuron), H4 (2 neurons), H5 (1 neuron) and H6 (2 neurons), all clusters which contain both Shox2 and Chx10 neurons. Thus, these results further support the ability of the H clustering method to distinguish populations of Shox2 neurons.

FIGURE 8

Figure 8. Validation of clustering analyses with Chx10 staining of Shox2 neurons. (A) Examples of a Shox2⁺Chx10⁻ (Ai) and a Shox2⁺Chx10⁺ (Aii) neuron filled with biocytin (cyan) during recording of an identified Shox2 neuron (tdTomato, red) and stained with a Chx10 antibody. (B) Three-dimensional representations of 8 Shox2 neurons processed for Chx10 identity (Shox2⁺Chx10⁻, diamonds; Shox2⁺Chx10⁺, squares) on the three principal components (same as Figure 2B, but different angle) for k-mean clusters (Bi) and hierarchical clustering (Bii). Ellipsoids are centered on the mean for the three first PC and sem-axes are generated by the standard deviation of each of the three first PCs.

Discussion

In the present study, electrophysiological properties were used to define populations of adult spinal Shox2 interneurons. Since Shox2 neurons overlap partly with Chx10-expressing neurons (Dougherty et al., 2013), we included a group of Chx10 neurons in our analysis. We classified the neurons based on their type of firing into 4 groups, and there were few significant differences between the groups in terms of electrophysiological properties. On the contrary, computational clustering methods delineated groups which could be readily distinguished by active and passive membrane properties. The k-means algorithm was run for 4 populations of neurons while hierarchical clustering analysis defined 6 populations. Interestingly, in each of the two clustering analyses, we found clusters containing exclusively Shox2 neurons, suggesting possible definition of Shox2⁺Chx10⁻ populations by electrophysiological properties. Finally, as preliminary validation, Chx10 antibody staining of a small subset of biocytin-filled and recovered Shox2 neurons showed that all (8/8) of the filled and post-processed neurons were appropriately found in the expected hierarchical (H) clusters based on Chx10 presence/absence. The k-means clustering corresponded to the expected in 7/8 cases. Taken together, the data demonstrate that it is possible to classify neurons expressing Shox2 in at least six different subpopulations based on active and passive membrane properties. These subpopulations correspond well with Chx10 absence/presence and may provide further divisions between neurons currently defined with intersectional genetics.

Classification by firing type

The firing type of neurons, defined by the response to current steps, is used to classify interneurons throughout the CNS (Lopez-Garcia and King, 1994; Prescott and De Koninck, 2002; Ruscheweyh and Sandkuhler, 2002; Dougherty and Hochman, 2008; Dai et al., 2009; Dougherty and Kiehn, 2010; Yasaka et al., 2010; Zhong et al., 2010; Dougherty and Chen, 2016). Prior studies have noted that Shox2 and Chx10 neurons are not entirely homogeneous (Dougherty and Kiehn, 2010; Zhong et al., 2010; Dougherty et al., 2013; Ampatzis et al., 2014; Garcia-Ramirez et al., 2021). There were significant differences in the distribution of neurons classified based on the type of firing. The majority of Shox2 neurons displayed a tonic response while the majority of Chx10 neurons displayed initial doublet firing. Adult Chx10 neurons have been shown to be predominantly tonic firing (Husch et al., 2015) which aligns with 83% of the Chx10 neurons in the present study firing throughout the current step. However, over half of the Chx10 neurons firing throughout a current step here had an initial doublet. The discrepancy may be due to the level of description, type of recording (whole-cell vs. perforated patch), or differences in neurons maintaining fluorescence expression in the mouse lines used (Chx10-CFP vs Chx10GFP). Initial burst (called phasic in some studies) and delayed firing types have been reported in neonatal Chx10 neurons (Dougherty and Kiehn, 2010; Zhong et al., 2010). The types of firing patterns observed are largely consistent between neonatal and adult but shifts in the distributions with maturation have been noted previously (Husch et al., 2015).

Unbiased cluster analysis to classify Shox2 and Chx10 neurons

The differences in the groups generated based on the type of firing was low. Additionally, a possible subclassification that merged Shox2 and Chx10 populations was not evident. Furthermore, subjectivity is another disadvantage for this type of classification. Unsupervised computational algorithms, have been used to unbiasedly classify neuronal populations based on physiological, molecular, and anatomical features (Dombeck et al., 2009; Karagiannis et al., 2009; Helm et al., 2013; Li et al., 2017; Martinez et al., 2017; Sathyamurthy et al., 2018; Mickelsen et al., 2019; Di Miceli et al., 2020). Two of the most prevalent algorithms implemented in such analyses are k-means and hierarchical clustering (Armananzas and Ascoli, 2015; Zeng and Sanes, 2017). Both analyses group individual objects (in this case neurons) based on similarities in input data (i.e. electrophysiological characteristics). The k-means algorithm requires a pre-established number of clusters (k value) set by the operator. We tested the k-mean algorithm in a range from 2 to 8 clusters, based on the idea of at least two populations of Shox2 neurons depending on the expression of Chx10 (Dougherty et al., 2013) and at least four groups based on the type of firing. We chose to proceed with 4 clusters because that resulted in the highest mean silhouette value (Allen et al., 2014; Vergara et al., 2020). However, we do not exclude the possibility of finding higher silhouette values with more than 8 clusters. The hierarchical cluster algorithm does not require the specification of the number of clusters prior to initial analysis. We found 6 natural divisions in the data based on a cutoff threshold below the value of the maximum inconsistency. Although the analyses resulted in a different number of clusters there was a high overlap between the two methods when clusters were directly compared.

A main motivation for the analyses was to determine if there was a way to identify subpopulations within a class of interneurons using electrophysiological properties rather than combinatorial genetics. The neurons expressing Shox2 can be divided into 2 subpopulations based on the expression of Chx10, Shox2⁺Chx10⁻ and Shox2⁺Chx10⁺. When considering both transcription factors, it should be noted that there are also Chx10 neurons which do not express Shox2 (Shox2⁻Chx10⁺). One may expect that electrophysiological properties are related to transcription factor expression and inferred function. Thus, there would be clusters that contain each of the possible combinations. If this is the case, the distribution of Shox2 and Chx10 neurons in the clusters identified by the k-means and hierarchical algorithms could potentially be used to predict the type of neurons (Shox2⁺Chx10⁻ or Shox2⁺Chx10⁺) that correspond to each cluster. For example, the lack of Chx10 neurons in k1, H1 and H2 clusters suggest that these clusters are composed of Shox2⁺Chx10⁻ neurons. H5 cluster could also be considered as a putative Shox2⁺Chx10⁻ cluster, since 93% of H5 cluster is comprised by Shox2 neurons, suggesting higher influence of Shox2 features than Chx10. As k3, k4, H3, H4 and H6 clusters have similar distributions of Shox2 neurons (82, 83, 75, 80 and 75%, respectively) as the total proportion of Shox2 neurons sampled (83%), we suggest that these clusters are composed of Shox2⁺Chx10⁺ neurons. In contrast, Shox2 neurons constitute 64% of the neurons in k2. Even though this percentage is lower than the total distribution of Shox2, this is not significantly different (binomial test, p = 0.09) and therefore we cannot classify them as Shox2⁻Chx10⁺ neurons.

We do not have any clusters of Shox2⁻Chx10⁺ neurons. This is most likely because there is a low number of Chx10 neurons in our analysis which resulted in an underrepresentation of the population of Shox2⁻Chx10⁺ neurons. Our primary focus here was the division of the Shox2 population and how it would match up with the subgroupings by transcription factor (Chx10) expression. A future analysis may include a more balanced sample of Chx10 neurons, and we expect that would generate clusters consisting of only Chx10 neurons.

In order to validate clusters obtained from both hierarchical and k-mean algorithms, previous studies have compared the output of both clusters and then corrected one or both algorithms (Karagiannis et al., 2009; McGarry et al., 2010; Perrenoud et al., 2012; Helm et al., 2013; Pohlkamp et al., 2014; Martinez et al., 2017). Here, we compared k-mean and hierarchical clusters and found that in 86% of the neurons were found in corresponding k-means and H-cluster pairs. Based on the distribution of Shox2 and Chx10 neurons in the k and H clusters, neurons can be arranged in a gradient from most Shox2 neuron-like to a most Chx10 neuron-like (Figure 9). The fact that k1, H1 and H2 clusters are composed only of Shox2 neurons suggests that these clusters have the features of Shox2⁺Chx10⁻ neurons. k3, k4, H3, H4 and H6 are composed of Shox2⁺Chx10⁺ neurons and the k2 cluster has an inclination toward a Chx10 phenotype. The arrangement displayed here could be considered as distinct populations of neurons or a phenotypic gradient of Shox2 to Chx10 neuronal features.

FIGURE 9

Figure 9. Arrangement of clusters in suggested order of gradient of Shox2-like to Chx10-like properties. k-means (above) and hierarchical (below) clusters arranged based on their similitudes with Shox2 (left and black portion of the middle bar) and Chx10 (right and red portion of the middle bar) neuronal characteristics.

Electrophysiological characteristics of Shox2 and Chx10 neurons

The passive and active electrophysiological properties of the neurons delineate firing responses of the cell to synaptic inputs (Russo and Hounsgaard, 1999; Butt et al., 2002; Smith and Perrier, 2006; Grillner and El Manira, 2020). Shox2 and Chx10 neurons are involved on the generation of the locomotor rhythm and pattern (Crone et al., 2008; Dougherty and Kiehn, 2010; Zhong et al., 2010, 2011; Dougherty et al., 2013; Dougherty and Ha, 2019). Therefore, the electrophysiological passive and active properties of the neuronal clusters and the differences between them should provide insights into the neuronal characteristics required for unique functional populations. We found that Shox2 neurons have lower input resistances and higher capacitance than Chx10 neurons, suggesting that Shox2 neurons are larger than Chx10 neurons. In fact, k1 and H2 clusters, that we consider to be Shox2⁺Chx10⁻ neurons also displayed low input resistances and high capacitances. On the other hand, the H4 cluster linked with a more Shox2⁺Chx10⁺-like phenotype had high input resistances and low capacitances.

One characteristic of motoneurons is their spike frequency adaptation in response to depolarizing current steps. Spike frequency adaptation is related to the inactivation of Na⁺ channels and functionally related to the initiation of muscle contraction (Miles et al., 2005; Brownstone, 2006). The initial doublet and tonic firing types observed here differed in spike frequency adaptation displayed. Spike frequency adaptation is evident in the low slope values from initial double compared with tonic neurons (Table 2). The significant differences observed in the percentage of type of firing between Shox2 and Chx10 are due to mainly to the larger proportion of Chx10 neurons with the initial doublet firing and higher incidence of the tonic type of firing for Shox2 neurons. Chx10 neurons also displayed longer action potentials and fAHPs in comparison to Shox2 (Table 2). In motoneurons, the fAHP has been linked to the activation of high threshold transient A-type K⁺ channels (Hess and El Manira, 2001; Harris-Warrick, 2002). fAHP contributes to the adaptation observed in initial doublet firing neurons (Baldissera and Gustafsson, 1974; Madison and Nicoll, 1984; Mrowczynski et al., 2015). Here, we found reduced fAHP amplitudes in neurons in the H1 cluster (Shox2⁺Chx10⁻) and longer fAHP durations in neurons in the k2 cluster (Shox2⁻Chx10⁺ enriched), which corresponds well with the proportion of initial doublet firing neurons in each group. Taken together, the computational clustering analysis performed here exposed potential novel differences between Shox2 populations that allow for hypotheses to be made regarding the characteristics of these populations. It is important to note that the electrophysiological recordings were performed at room temperature, as with most recordings of adult locomotor-related neurons in vitro (Husch et al., 2012; Mitra and Brownstone, 2012; Hadzipasic et al., 2014; Smith and Brownstone, 2020). The dynamics and activity of ion channels may differ at physiological temperatures (Russell et al., 1994; Li and Burke, 2001; Robertson and Money, 2012; O'Leary and Marder, 2016). Although it is likely to affect values used here for classification, diversity in properties is seen in spinal neurons in vivo (Graham et al., 2004) and similar clustering methodologies should be applicable.

Computational clustering impact

Recent advances in computational technologies have increased our capabilities to cluster spinal neurons, based on molecular profiles from RNA sequence analysis (Li et al., 2017; Menon, 2018; Sathyamurthy et al., 2018; Mickelsen et al., 2019). RNA sequence analysis provides us information about the possibilities of neurons to express proteins that in turn will form ion channels or define neurotransmitter machinery. But, the ways by which these proteins actively participate during behavior is dynamic (Hager et al., 2009; Tanay and Regev, 2017). Classification based on electrophysiological properties, should be considered as a complement to the biomolecular identification of neurons. Clustering methods with electrophysiological properties have been used to identify neurons from other areas of the central nervous system (Karagiannis et al., 2009; Perrenoud et al., 2012; Helm et al., 2013; Martinez et al., 2017; Di Miceli et al., 2020). Here, in addition to computational clustering to classify Shox2 and Chx10 neurons, we were able to test the predicted clusters with immunocytochemistry approaches with high accuracy. The results observed and further validated allow some of the electrophysiological properties that potentially conferred the capacity to a group of neurons to perform specific locomotor task to be inferred.

Limitations of computational clustering techniques include the necessity of a large sample of neurons for the analysis (Armananzas and Ascoli, 2015). Additionally, a large number of parameters should be measured for each neuron in the sample (Armananzas and Ascoli, 2015). Although here we demonstrate that it is possible with a more limited number of parameters, the larger number was collapsed by determining those that were highly correlated. Another limitation is that most clustering algorithms (including the two used here) require the number of clusters generated to be pre-determined. In order to minimize subjectivity in the number of clusters, these can be determined based on mathematical criteria, as we did here. For k-means, the elbow method considering silhouette values can be used (Allen et al., 2014; Vergara et al., 2020) and a cutoff can be established considering inconsistency coefficients for hierarchical clustering (McGarry et al., 2010). However, computational clustering overcomes limitations of other methods by being unbiased. For example, the classification of neurons based on type of firing is subjective and, therefore, susceptible to discrepancies. Both k-means and hierarchical clustering are easy to implement, the computations are fast, and they consider all of the data including apparent outliers (Karagiannis et al., 2009; Armananzas and Ascoli, 2015).

The computational unsupervised clustering technique used here would be useful to classify subpopulations of spinal interneurons in multiple applications. It could be used when recording from neurons in in vitro preparations to avoid the necessity of combinatorial transgenic mice. Additionally, in vivo recordings often performed blindly (Tao et al., 2015; Lee and Lee, 2017) but more intact preparations are important for the understanding of neuronal spinal physiology in awake animal models. Furthermore, electrophysiological methods to classify neurons could be combined with other types of classifications, such as large scale snRNA-sequencing (Sathyamurthy et al., 2018; Patterson-Cross et al., 2021) or anatomical criteria (Karagiannis et al., 2009; McGarry et al., 2010), to obtain a multidimensional identification of spinal neuronal populations. Finally, the cluster analysis performed here could eventually be converted to a semi-supervised machine learning model by training the algorithm to classify a larger scale of spinal neurons (Buccino et al., 2018; Paninski and Cunningham, 2018) that is capable of considering additional factors including behavioral responses, connectivity, and/or neuromodulation, thereby integrating cellular physiology to behavior.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Ethics statement

The animal study was reviewed and approved by Institutional Animal Care and Use Committee at Drexel University.

Author contributions

DG-R and KD designed the study, drafted the manuscript, and wrote the manuscript. DG-R, SS, JM, and NH performed the electrophysiological recordings and individual neuron analyses. SS performed the immunohistochemistry experiments. DG-R ran the cluster analyses and performed all statistical analyses. All authors read, edited, and approved the final manuscript.

Funding

This work was supported by NIH NINDS R01 NS104194 and R21 NS118226 (KD), the Edward Jekkal Muscular Dystrophy Association Fellowship (DG-R), and NIH T32 NS121768 (SS).

Acknowledgments

The authors are grateful to Lihua Yao for technical assistance, Nicholas Stachowski for comments on the manuscript, and Simon Danner and the Marion Murray Spinal Cord Research Center for discussions related to the results and analyses.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncir.2022.957084/full#supplementary-material

Supplementary Figure 1. Correspondence of k and H clusters. Correspondence of k-clusters (left squares) and H-clusters (right squares) for each of the 171 Shox2 neurons (cell identification number on the left, Shox2 in black, Chx10 in red). k-clusters (left boxes), k1 (purple), k2 (blue), k3 (yellow), k4 (gray), H-clusters (right boxes), H1 (blue), H2 (orange), H3 (yellow), H4 (purple), H5 (green), and H6 (cyan).

Supplementary Table 1. Statistical comparisons of Shox2 and Chx10 neuronal properties presented in Figure 3.

Supplementary Table 2. Statistical comparisons of membrane properties by firing type of Shox2 and Chx10 neurons presented in Figure 4.

Supplementary Table 3. Comparison of cellular properties in neurons separated by k-means clustering presented in Figure 5.

Supplementary Table 4. Comparison of cellular properties in neurons separated by hierarchical clustering presented in Figure 6.

References

Abraira, V. E., Kuehn, E. D., Chirila, A. M., Springel, M. W., Toliver, A. A., Zimmerman, A. L., et al. (2017). The cellular and synaptic architecture of the mechanosensory dorsal horn. Cell 168, 295–310.e219. doi: 10.1016/j.cell.2016.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, E. A., Damaraju, E., Plis, S. M., Erhardt, E. B., Eichele, T., and Calhoun, V. D. (2014). Tracking whole-brain connectivity dynamics in the resting state. Cere. Cortex 24, 663–676. doi: 10.1093/cercor/bhs352

PubMed Abstract | CrossRef Full Text | Google Scholar

Al-Mosawie, A., Wilson, J. M., and Brownstone, R. M. (2007). Heterogeneity of V2-derived interneurons in the adult mouse spinal cord. Eur. J. Neurosci. 26, 3003–3015. doi: 10.1111/j.1460-9568.2007.05907.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Alvarez, F. J., Jonas, P. C., Sapir, T., Hartley, R., Berrocal, M. C., Geiman, E. J., et al. (2005). Postnatal phenotype and localization of spinal cord V1 derived interneurons. J. Comp. Neurol. 493, 177–192. doi: 10.1002/cne.20711

PubMed Abstract | CrossRef Full Text | Google Scholar

Ampatzis, K., Song, J., Ausborn, J., and El Manira, A. (2014). Separate microcircuit modules of distinct v2a interneurons and motoneurons control the speed of locomotion. Neuron 83, 934–943. doi: 10.1016/j.neuron.2014.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Armananzas, R., and Ascoli, G. A. (2015). Towards the automatic classification of neurons. Trends Neurosci. 38, 307–318. doi: 10.1016/j.tins.2015.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Baldissera, F., and Gustafsson, B. (1974). Firing behaviour of a neurone model based on the afterhyperpolarization conductance time course and algebraical summation. Adaptation and steady state firing. Acta. Physiol. Scand. 92, 27–47. doi: 10.1111/j.1748-1716.1974.tb05720.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bikoff, J. B., Gabitto, M. I., Rivard, A. F., Drobac, E., Machado, T. A., Miri, A., et al. (2016). Spinal inhibitory interneuron diversity delineates variant motor microcircuits. Cell 165, 207–219. doi: 10.1016/j.cell.2016.01.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Borowska, J., Jones, C. T., Zhang, H., Blacklaws, J., Goulding, M., and Zhang, Y. (2013). Functional subpopulations of V3 interneurons in the mature mouse spinal cord. J. Neurosci. 33, 18553–18565. doi: 10.1523/JNEUROSCI.2005-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Brownstone, R. M. (2006). Beginning at the end: repetitive firing properties in the final common pathway. Prog. Neurobiol. 78, 156–172. doi: 10.1016/j.pneurobio.2006.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Brownstone, R. M., and Wilson, J. M. (2008). Strategies for delineating spinal locomotor rhythm-generating networks and the possible role of Hb9 interneurones in rhythmogenesis. Brain Res. Rev. 57, 64–76. doi: 10.1016/j.brainresrev.2007.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Buccino, A. P., Kordovan, M., Ness, T. V., Merkt, B., Hafliger, P. D., Fyhn, M., et al. (2018). Combining biophysical modeling and deep learning for multielectrode array neuron localization and classification. J. Neurophysiol. 120, 1212–1232. doi: 10.1152/jn.00210.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Butt, S. J., Harris-Warrick, R. M., and Kiehn, O. (2002). Firing properties of identified interneuron populations in the mammalian hindlimb central pattern generator. J. Neurosci. 22, 9961–9971. doi: 10.1523/JNEUROSCI.22-22-09961.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Caldeira, V., Dougherty, K. J., Borgius, L., and Kiehn, O. (2017). Spinal Hb9::Cre-derived excitatory interneurons contribute to rhythm generation in the mouse. Sci. Rep. 7, 41369. doi: 10.1038/srep41369

PubMed Abstract | CrossRef Full Text | Google Scholar

Crone, S. A., Quinlan, K. A., Zagoraiou, L., Droho, S., Restrepo, C. E., Lundfald, L., et al. (2008). Genetic ablation of V2a ipsilateral interneurons disrupts left-right locomotor coordination in mammalian spinal cord. Neuron 60, 70–83. doi: 10.1016/j.neuron.2008.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Crone, S. A., Zhong, G., Harris-Warrick, R., and Sharma, K. (2009). In mice lacking V2a interneurons, gait depends on speed of locomotion. J. Neurosci. 29, 7098–7109. doi: 10.1523/JNEUROSCI.1206-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, Y., Carlin, K. P., Li, Z., McMahon, D. G., Brownstone, R. M., and Jordan, L. M. (2009). Electrophysiological and pharmacological properties of locomotor activity-related neurons in cfos-EGFP mice. J. Neurophysiol. 102, 3365–3383. doi: 10.1152/jn.00265.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Deska-Gauthier, D., Borowska-Fielding, J., Jones, C. T., and Zhang, Y. (2020). The Temporal neurogenesis patterning of spinal p3-V3 interneurons into divergent subpopulation assemblies. J. Neurosci. 40, 1440–1452. doi: 10.1523/JNEUROSCI.1518-19.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Miceli, M., Husson, Z., Ruel, P., Laye, S., Cota, D., Fioramonti, X., et al. (2020). In silico hierarchical clustering of neuronal populations in the rat ventral tegmental area based on extracellular electrophysiological properties. Front. Neural. Circuits 14, 51. doi: 10.3389/fncir.2020.00051

PubMed Abstract | CrossRef Full Text | Google Scholar

Dombeck, D. A., Graziano, M. S., and Tank, D. W. (2009). Functional clustering of neurons in motor cortex determined by cellular resolution imaging in awake behaving mice. J. Neurosci. 29, 13751–13760. doi: 10.1523/JNEUROSCI.2985-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Dougherty, K. J., and Ha, N. T. (2019). The rhythm section: an update on spinal interneurons setting the beat for mammalian locomotion. Curr. Opin. Physiol. 8, 84–93. doi: 10.1016/j.cophys.2019.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Dougherty, K. J., and Hochman, S. (2008). Spinal cord injury causes plasticity in a subpopulation of lamina I GABAergic interneurons. J. Neurophysiol. 100, 212–223. doi: 10.1152/jn.01104.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Dougherty, K. J., and Kiehn, O. (2010). Firing and cellular properties of V2a interneurons in the rodent spinal cord. J. Neurosci. 30, 24–37. doi: 10.1523/JNEUROSCI.4821-09.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Dougherty, K. J., Zagoraiou, L., Satoh, D., Rozani, I., Doobar, S., Arber, S., et al. (2013). Locomotor rhythm generation linked to the output of spinal shox2 excitatory interneurons. Neuron 80, 920–933. doi: 10.1016/j.neuron.2013.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Dougherty, P. M., and Chen, J. (2016). Relationship of membrane properties, spike burst responses, laminar location, and functional class of dorsal horn neurons recorded in vitro. J. Neurophysiol. 116, 1137–1151. doi: 10.1152/jn.00187.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Dyck, J., Lanuza, G. M., and Gosgnach, S. (2012). Functional characterization of dI6 interneurons in the neonatal mouse spinal cord. J. Neurophysiol. 107, 3256–3266. doi: 10.1152/jn.01132.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Falgairolle, M., and O'Donovan, M. J. (2021). Optogenetic activation of V1 interneurons reveals the multimodality of spinal locomotor networks in the neonatal mouse. J. Neurosci. 41, 8545–8561. doi: 10.1523/JNEUROSCI.0875-21.2021

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia-Ramirez, D. L., Ha, N. T. B., Bibu, S., Stachowski, N. J., and Dougherty, K. J. (2021). Spinal cord injury alters spinal Shox2 interneurons by enhancing excitatory synaptic input and serotonergic modulation while maintaining intrinsic properties in mouse. J. Neurosci. 41, 5833–5848. doi: 10.1523/JNEUROSCI.1576-20.2021

PubMed Abstract | CrossRef Full Text | Google Scholar

Gatto, G., Smith, K. M., Ross, S. E., and Goulding, M. (2019). Neuronal diversity in the somatosensory system: bridging the gap between cell type and function. Curr. Opin. Neurobiol. 56, 167–174. doi: 10.1016/j.conb.2019.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Gobel, S. (1975). Golgi studies in the substantia gelatinosa neurons in the spinal trigeminal nucleus. J. Comp. Neurol. 162, 397–415. doi: 10.1002/cne.901620308

PubMed Abstract | CrossRef Full Text | Google Scholar

Gosgnach, S., Lanuza, G. M., Butt, S. J., Saueressig, H., Zhang, Y., Velasquez, T., et al. (2006). V1 spinal neurons regulate the speed of vertebrate locomotor outputs. Nature 440, 215–219. doi: 10.1038/nature04545

PubMed Abstract | CrossRef Full Text | Google Scholar

Goulding, M. (2009). Circuits controlling vertebrate locomotion: moving in a new direction. Nat. Rev. Neurosci. 10, 507–518. doi: 10.1038/nrn2608

PubMed Abstract | CrossRef Full Text | Google Scholar

Goulding, M., Bourane, S., Garcia-Campmany, L., Dalet, A., and Koch, S. (2014). Inhibition downunder: an update from the spinal cord. Curr. Opin. Neurobiol. 26, 161–166. doi: 10.1016/j.conb.2014.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Graham, B. A., Brichta, A. M., and Callister, R. J. (2004). In vivo responses of mouse superficial dorsal horn neurones to both current injection and peripheral cutaneous stimulation. J. Physiol. 561, 749–763. doi: 10.1113/jphysiol.2004.072645

PubMed Abstract | CrossRef Full Text | Google Scholar

Grillner, S., and El Manira, A. (2020). Current principles of motor control, with special reference to vertebrate locomotion. Physiol. Rev. 100, 271–320. doi: 10.1152/physrev.00015.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Grudt, T. J., and Perl, E. R. (2002). Correlations between neuronal morphology and electrophysiological features in the rodent superficial dorsal horn. J. Physiol. 540, 189–207. doi: 10.1113/jphysiol.2001.012890

PubMed Abstract | CrossRef Full Text | Google Scholar

Ha, N. T., and Dougherty, K. J. (2018). Spinal Shox2 interneuron interconnectivity related to function and development. Elife. 7:e42519. doi: 10.7554/eLife.42519.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Hadzipasic, M., Tahvildari, B., Nagy, M., Bian, M., Horwich, A. L., and McCormick, D. A. (2014). Selective degeneration of a physiological subtype of spinal motor neuron in mice with SOD1-linked ALS. Proc. Natl. Acad. Sci. U.S.A. 111, 16883–16888. doi: 10.1073/pnas.1419497111

PubMed Abstract | CrossRef Full Text | Google Scholar

Hager, G. L., McNally, J. G., and Misteli, T. (2009). Transcription dynamics. Mol. Cell. 35, 741–753. doi: 10.1016/j.molcel.2009.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Haring, M., Zeisel, A., Hochgerner, H., Rinwa, P., Jakobsson, J. E. T., Lonnerberg, P., et al. (2018). Neuronal atlas of the dorsal horn defines its architecture and links sensory input to transcriptional cell types. Nat. Neurosci. 21, 869–880. doi: 10.1038/s41593-018-0141-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, P. J., Jankowska, E., and Zytnicki, D. (1986). Lamina VIII interneurones interposed in crossed reflex pathways in the cat. J. Physiol. 371, 147–166. doi: 10.1113/jphysiol.1986.sp015965

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris-Warrick, R. M. (2002). Voltage-sensitive ion channels in rhythmic motor systems. Curr. Opin. Neurobiol. 12, 646–651. doi: 10.1016/S0959-4388(02)00377-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris-Warrick, R. M. (2010). General principles of rhythmogenesis in central pattern generator networks. Prog. Brain Res. 187, 213–222. doi: 10.1016/B978-0-444-53613-6.00014-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayashi, M., Hinckley, C. A., Driscoll, S. P., Moore, N. J., Levine, A. J., Hilde, K. L., et al. (2018). Graded arrays of spinal and supraspinal V2a interneuron subtypes underlie forelimb and hindlimb motor control. Neuron 97, 869–884.e865. doi: 10.1016/j.neuron.2018.01.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Helm, J., Akgul, G., and Wollmuth, L. P. (2013). Subgroups of parvalbumin-expressing interneurons in layers 2/3 of the visual cortex. J. Neurophysiol. 109, 1600–1613. doi: 10.1152/jn.00782.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Hess, D., and El Manira, A. (2001). Characterization of a high-voltage-activated IA current with a role in spike timing and locomotor pattern generation. Proc. Natl. Acad. Sci. U.S.A. 98, 5276–5281. doi: 10.1073/pnas.091096198

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, A., Noga, B. R., Carr, P. A., Fedirchuk, B., and Jordan, L. M. (2000). Spinal cholinergic neurons activated during locomotion: localization and electrophysiological characterization. J. Neurophysiol. 83, 3537–3547. doi: 10.1152/jn.2000.83.6.3537

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, D. I., and Todd, A. J. (2020). Central nervous system targets: inhibitory interneurons in the spinal cord. Neurotherapeutics 17, 874–885. doi: 10.1007/s13311-020-00936-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Hultborn, H., and Nielsen, J. B. (2007). Spinal control of locomotion–from cat to man. Acta. Physiol. 189, 111–121. doi: 10.1111/j.1748-1716.2006.01651.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Husch, A., Dietz, S. B., Hong, D. N., and Harris-Warrick, R. M. (2015). Adult spinal V2a interneurons show increased excitability and serotonin-dependent bistability. J. Neurophysiol. 113, 1124–1134. doi: 10.1152/jn.00741.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Husch, A., Van Patten, G. N., Hong, D. N., Scaperotti, M. M., Cramer, N., and Harris-Warrick, R. M. (2012). Spinal cord injury induces serotonin supersensitivity without increasing intrinsic excitability of mouse V2a interneurons. J. Neurosci. 32, 13145–13154. doi: 10.1523/JNEUROSCI.2995-12.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Jankowska, E. (2008). Spinal interneuronal networks in the cat: elementary components. Brain Res. Rev. 57, 46–55. doi: 10.1016/j.brainresrev.2007.06.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Karagiannis, A., Gallopin, T., David, C., Battaglia, D., Geoffroy, H., Rossier, J., et al. (2009). Classification of NPY-expressing neocortical interneurons. J. Neurosci. 29, 3642–3659. doi: 10.1523/JNEUROSCI.0058-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiehn, O. (2016). Decoding the organization of spinal circuits that control locomotion. Nat. Rev. Neurosci. 17, 224–238. doi: 10.1038/nrn.2016.9

PubMed Abstract | CrossRef Full Text | Google Scholar

Konnerth, A., Keller, B. U., and Lev-Tov, A. (1990). Patch clamp analysis of excitatory synapses in mammalian spinal cord slices. Pflugers Arch. 417, 285–290. doi: 10.1007/BF00370994

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, D., and Lee, A. K. (2017). In vivo patch-clamp recording in awake head-fixed rodents. Cold Spring Harb. Protoc. 2017, prot095802. doi: 10.1101/pdb.prot095802

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, R. H., and Heckman, C. J. (2000). Adjustable amplification of synaptic input in the dendrites of spinal motoneurons in vivo. J. Neurosci. 20, 6734–6740. doi: 10.1523/JNEUROSCI.20-17-06734.2000

PubMed Abstract | CrossRef Full Text | Google Scholar

Levine, A. J., Lewallen, K. A., and Pfaff, S. L. (2012). Spatial organization of cortical and spinal neurons controlling motor behavior. Curr. Opin. Neurobiol. 22, 812–821. doi: 10.1016/j.conb.2012.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Horns, F., Wu, B., Xie, Q., Li, J., Li, T., et al. (2017). Classifying drosophila olfactory projection neuron subtypes by single-cell RNA sequencing. Cell 171, 1206–1220.e1222. doi: 10.1016/j.cell.2017.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., and Burke, R. E. (2001). Short-term synaptic depression in the neonatal mouse spinal cord: effects of calcium and temperature. J. Neurophysiol. 85, 2047–2062. doi: 10.1152/jn.2001.85.5.2047

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopez-Garcia, J. A., and King, A. E. (1994). Membrane properties of physiologically classified rat dorsal horn neurons in vitro: correlation with cutaneous sensory afferent input. Eur J Neurosci 6, 998–1007. doi: 10.1111/j.1460-9568.1994.tb00594.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Madisen, L., Zwingman, T. A., Sunkin, S. M., Oh, S. W., Zariwala, H. A., Gu, H., et al. (2010). A robust and high-throughput Cre reporting and characterization system for the whole mouse brain. Nat. Neurosci. 13, 133–140. doi: 10.1038/nn.2467

PubMed Abstract | CrossRef Full Text | Google Scholar

Madison, D. V., and Nicoll, R. A. (1984). Control of the repetitive discharge of rat CA 1 pyramidal neurones in vitro. J. Physiol. 354, 319–331. doi: 10.1113/jphysiol.1984.sp015378

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez, J. J., Rahsepar, B., and White, J. A. (2017). Anatomical and electrophysiological clustering of superficial medial entorhinal cortex interneurons. eNeuro. 4:ENEURO.0263-16.2017. doi: 10.1523/ENEURO.0263-16.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Maxwell, D. J., Belle, M. D., Cheunsuang, O., Stewart, A., and Morris, R. (2007). Morphology of inhibitory and excitatory interneurons in superficial laminae of the rat dorsal horn. J. Physiol. 584, 521–533. doi: 10.1113/jphysiol.2007.140996

PubMed Abstract | CrossRef Full Text | Google Scholar

McGarry, L. M., Packer, A. M., Fino, E., Nikolenko, V., Sippy, T., and Yuste, R. (2010). Quantitative classification of somatostatin-positive neocortical interneurons identifies three interneuron subtypes. Front. Neural. Circuits 4, 12. doi: 10.3389/fncir.2010.00012

PubMed Abstract | CrossRef Full Text | Google Scholar

Menon, V. (2018). Clustering single cells: a review of approaches on high-and low-depth single-cell RNA-seq data. Brief Funct. Genomics 17, 240–245. doi: 10.1093/bfgp/elx044

PubMed Abstract | CrossRef Full Text | Google Scholar

Mickelsen, L. E., Bolisetty, M., Chimileski, B. R., Fujita, A., Beltrami, E. J., Costanzo, J. T., et al. (2019). Single-cell transcriptomic analysis of the lateral hypothalamic area reveals molecularly distinct populations of inhibitory and excitatory neurons. Nat. Neurosci. 22, 642–656. doi: 10.1038/s41593-019-0349-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Miles, G. B., Dai, Y., and Brownstone, R. M. (2005). Mechanisms underlying the early phase of spike frequency adaptation in mouse spinal motoneurones. J. Physiol. 566, 519–532. doi: 10.1113/jphysiol.2005.086033

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitra, P., and Brownstone, R. M. (2012). An in vitro spinal cord slice preparation for recording from lumbar motoneurons of the adult mouse. J. Neurophysiol. 107, 728–741. doi: 10.1152/jn.00558.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Mrowczynski, W., Celichowski, J., Raikova, R., and Krutki, P. (2015). Physiological consequences of doublet discharges on motoneuronal firing and motor unit force. Front. Cell Neurosci. 9, 81. doi: 10.3389/fncel.2015.00081

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Leary, T., and Marder, E. (2016). Temperature-robust neural function from activity-dependent ion channel regulation. Curr. Biol. 26, 2935–2941. doi: 10.1016/j.cub.2016.08.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Paninski, L., and Cunningham, J. P. (2018). Neural data science: accelerating the experiment-analysis-theory cycle in large-scale neuroscience. Curr. Opin. Neurobiol. 50, 232–241. doi: 10.1016/j.conb.2018.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Patterson-Cross, R. B., Levine, A. J., and Menon, V. (2021). Selecting single cell clustering parameter values using subsampling-based robustness metrics. BMC Bioinform. 22, 39. doi: 10.1186/s12859-021-03957-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Perrenoud, Q., Geoffroy, H., Gauthier, B., Rancillac, A., Alfonsi, F., Kessaris, N., et al. (2012). Characterization of type I and type II nNOS-expressing interneurons in the barrel cortex of mouse. Front Neural Circuits 6, 36. doi: 10.3389/fncir.2012.00036

PubMed Abstract | CrossRef Full Text | Google Scholar

Petko, M., Veress, G., Vereb, G., Storm-Mathisen, J., and Antal, M. (2004). Commissural propriospinal connections between the lateral aspects of laminae III-IV in the lumbar spinal cord of rats. J. Comp. Neurol. 480, 364–377. doi: 10.1002/cne.20356

PubMed Abstract | CrossRef Full Text | Google Scholar

Pierani, A., Moran-Rivard, L., Sunshine, M. J., Littman, D. R., Goulding, M., and Jessell, T. M. (2001). Control of interneuron fate in the developing spinal cord by the progenitor homeodomain protein Dbx1. Neuron 29, 367–384. doi: 10.1016/S0896-6273(01)00212-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Pohlkamp, T., David, C., Cauli, B., Gallopin, T., Bouche, E., Karagiannis, A., et al. (2014). Characterization and distribution of Reelin-positive interneuron subtypes in the rat barrel cortex. Cereb. Cortex 24, 3046–3058. doi: 10.1093/cercor/bht161

PubMed Abstract | CrossRef Full Text | Google Scholar

Prescott, S. A., and De Koninck, Y. (2002). Four cell types with distinctive membrane properties and morphologies in lamina I of the spinal dorsal horn of the adult rat. J. Physiol. 539, 817–836. doi: 10.1113/jphysiol.2001.013437

PubMed Abstract | CrossRef Full Text | Google Scholar

Puskar, Z., and Antal, M. (1997). Localization of last-order premotor interneurons in the lumbar spinal cord of rats. J. Comp. Neurol. 389, 377–389. doi: 10.1002/(sici)1096-9861(19971222)389:3<377::aid-cne2>3.0.co;2-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Robertson, R. M., and Money, T. G. (2012). Temperature and neuronal circuit function: compensation, tuning and tolerance. Curr. Opin. Neurobiol. 22, 724–734. doi: 10.1016/j.conb.2012.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruscheweyh, R., and Sandkuhler, J. (2002). Lamina-specific membrane and discharge properties of rat spinal dorsal horn neurones in vitro. J. Physiol. 541, 231–244. doi: 10.1113/jphysiol.2002.017756

PubMed Abstract | CrossRef Full Text | Google Scholar

Russell, S. N., Publicover, N. G., Hart, P. J., Carl, A., Hume, J. R., Sanders, K. M., et al. (1994). Block by 4-aminopyridine of a Kv1.2 delayed rectifier K+ current expressed in Xenopus oocytes. J. Physiol. 481, 571–584. doi: 10.1113/jphysiol.1994.sp020464

PubMed Abstract | CrossRef Full Text | Google Scholar

Russo, R. E., and Hounsgaard, J. (1999). Dynamics of intrinsic electrophysiological properties in spinal cord neurones. Prog. Biophys. Mol. Biol. 72, 329–365. doi: 10.1016/S0079-6107(99)00011-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Sathyamurthy, A., Johnson, K. R., Matson, K. J. E., Dobrott, C. I., Li, L., Ryba, A. R., et al. (2018). massively parallel single nucleus transcriptional profiling defines spinal cord neurons and their activity during behavior. Cell. Rep. 22, 2216–2225. doi: 10.1016/j.celrep.2018.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Shevtsova, N. A., and Rybak, I. A. (2016). Organization of flexor-extensor interactions in the mammalian spinal cord: insights from computational modelling. J. Physiol. 594, 6117–6131. doi: 10.1113/JP272437

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, C. C., and Brownstone, R. M. (2020). Spinal motoneuron firing properties mature from rostral to caudal during postnatal development of the mouse. J. Physiol. 598, 5467–5485. doi: 10.1113/JP280274

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, M., and Perrier, J. F. (2006). Intrinsic properties shape the firing pattern of ventral horn interneurons from the spinal cord of the adult turtle. J. Neurophysiol. 96, 2670–2677. doi: 10.1152/jn.00609.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, M. R., Sun, Y., Bryson, A., Gill, G. N., Evans, S. M., and Pfaff, S. L. (2009). Islet-to-LMO stoichiometries control the function of transcription complexes that specify motor neuron and V2a interneuron identity. Development 136, 2923–2932. doi: 10.1242/dev.037986

PubMed Abstract | CrossRef Full Text | Google Scholar

Stachowski, N. J., and Dougherty, K. J. (2021). Spinal inhibitory interneurons: gatekeepers of sensorimotor pathways. Int. J. Mol. Sci. 22, 2667. doi: 10.3390/ijms22052667

PubMed Abstract | CrossRef Full Text | Google Scholar

Talpalar, A. E., Bouvier, J., Borgius, L., Fortin, G., Pierani, A., and Kiehn, O. (2013). Dual-mode operation of neuronal networks involved in left-right alternation. Nature 500, 85–88. doi: 10.1038/nature12286

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanay, A., and Regev, A. (2017). Scaling single-cell genomics from phenomenology to mechanism. Nature 541, 331–338. doi: 10.1038/nature21350

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, C., Zhang, G., Xiong, Y., and Zhou, Y. (2015). Functional dissection of synaptic circuits: in vivo patch-clamp recording in neuroscience. Front. Neural Circuits 9, 23. doi: 10.3389/fncir.2015.00023

PubMed Abstract | CrossRef Full Text | Google Scholar

Tazerart, S., Vinay, L., and Brocard, F. (2008). The persistent sodium current generates pacemaker activities in the central pattern generator for locomotion and regulates the locomotor rhythm. J. Neurosci. 28, 8577–8589. doi: 10.1523/JNEUROSCI.1437-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Todd, A. J. (2010). Neuronal circuitry for pain processing in the dorsal horn. Nat. Rev. Neurosci. 11, 823–836. doi: 10.1038/nrn2947

PubMed Abstract | CrossRef Full Text | Google Scholar

Vergara, V. M., Salman, M., Abrol, A., Espinoza, F. A., and Calhoun, V. D. (2020). Determining the number of states in dynamic functional connectivity using cluster validity indexes. J. Neurosci. Methods 337, 108651. doi: 10.1016/j.jneumeth.2020.108651

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, J. M., Cowan, A. I., and Brownstone, R. M. (2007). Heterogeneous electrotonic coupling and synchronization of rhythmic bursting activity in mouse Hb9 interneurons. J. Neurophysiol. 98, 2370–2381. doi: 10.1152/jn.00338.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, J. M., Hartley, R., Maxwell, D. J., Todd, A. J., Lieberam, I., Kaltschmidt, J. A., et al. (2005). Conditional rhythmicity of ventral spinal interneurons defined by expression of the Hb9 homeodomain protein. J. Neurosci. 25, 5710–5719. doi: 10.1523/JNEUROSCI.0274-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Yasaka, T., Tiong, S. Y. X., Hughes, D. I., Riddell, J. S., and Todd, A. J. (2010). Populations of inhibitory and excitatory interneurons in lamina II of the adult rat spinal dorsal horn revealed by a combined electrophysiological and anatomical approach. Pain 151, 475–488. doi: 10.1016/j.pain.2010.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, H., and Sanes, J. R. (2017). Neuronal cell-type classification: challenges, opportunities and the path forward. Nat. Rev. Neurosci. 18, 530–546. doi: 10.1038/nrn.2017.85

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Lanuza, G. M., Britz, O., Wang, Z., Siembab, V. C., Zhang, Y., et al. (2014). V1 and v2b interneurons secure the alternating flexor-extensor motor activity mice require for limbed locomotion. Neuron 82, 138–150. doi: 10.1016/j.neuron.2014.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Narayan, S., Geiman, E., Lanuza, G. M., Velasquez, T., Shanks, B., et al. (2008). V3 spinal neurons establish a robust and balanced locomotor rhythm during walking. Neuron 60, 84–96. doi: 10.1016/j.neuron.2008.09.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, G., Droho, S., Crone, S. A., Dietz, S., Kwan, A. C., Webb, W. W., et al. (2010). Electrophysiological characterization of V2a interneurons and their locomotor-related activity in the neonatal mouse spinal cord. J. Neurosci. 30, 170–182. doi: 10.1523/JNEUROSCI.4849-09.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, G., Sharma, K., and Harris-Warrick, R. M. (2011). Frequency-dependent recruitment of V2a interneurons during fictive locomotion in the mouse spinal cord. Nat. Commun. 2, 274. doi: 10.1038/ncomms1276

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, G., Shevtsova, N. A., Rybak, I. A., and Harris-Warrick, R. M. (2012). Neuronal activity in the isolated mouse spinal cord during spontaneous deletions in fictive locomotion: insights into locomotor central pattern generator organization. J. Physiol. 590, 4735–4759. doi: 10.1113/jphysiol.2012.240895

PubMed Abstract | CrossRef Full Text | Google Scholar

Ziskind-Conhaim, L., and Hochman, S. (2017). Diversity of molecularly defined spinal interneurons engaged in mammalian locomotor pattern generation. J. Neurophysiol. 118, 2956–2974. doi: 10.1152/jn.00322.2017

PubMed Abstract | CrossRef Full Text | Google Scholar