Identification of marker genes expressed in specific cell types is essential for the genetic dissection of neural circuits. Here we report a new strategy for classifying heterogeneous populations of neurons into functionally distinct types and for identifying associated marker genes. Quantitative single-cell expression profiling of genes related to neurotransmitters and ion channels enables functional classification of neurons; transcript profiles for marker gene candidates identify molecular handles for manipulating each cell type. We apply this strategy to the mouse medial vestibular nucleus (MVN), which comprises several types of neurons subserving cerebellar-dependent learning in the vestibulo-ocular reflex. Ion channel gene expression differed both qualitatively and quantitatively across cell types and could distinguish subtle differences in intrinsic electrophysiology. Single-cell transcript profiling of MVN neurons established six functionally distinct cell types and associated marker genes. This strategy is applicable throughout the nervous system and could facilitate the use of molecular genetic tools to examine the behavioral roles of distinct neuronal populations.
A major neuroscientific challenge is to determine how cellular signaling and plasticity mediate behavioral performance and learning. The ability to link neuronal firing directly with quantifiable behavioral performance and sensory-motor learning makes the vestibulo-ocular reflex (VOR) a particularly attractive model. The VOR stabilizes images on the retina during self-motion by transforming head motion signals from the inner ear into oculomotor commands. The primary circuit is simple: a “three neuron arc” comprising vestibular ganglion cells, medial vestibular nucleus (MVN) neurons, and ocular motoneurons mediates gaze stability in the horizontal plane via a pair of extraocular muscles. Excellent behavioral performance is maintained throughout life via cerebellar learning (du Lac et al., 1995; Blazquez et al., 2004; Broussard and Kassardjian, 2004; Schubert and Zee, 2010). Progress at the cellular mechanistic level has been limited, however, by the challenges of classifying and manipulating functionally distinct cell types which are spatially intermingled within central vestibular nuclei.
In vivo recordings have revealed several types of vestibular nucleus neurons with differential synaptic connections and/or firing responses during head and eye movements (Highstein and Ito, 1971; Keller and Kamath, 1975; Lisberger and Miles, 1980; McCrea et al., 1987; Scudder and Fuchs, 1992). Cerebellar learning was originally thought to be mediated by Purkinje cell synapses onto one of these cell types (Lisberger et al., 1994; Ramachandran and Lisberger, 2008). Although initial in vitro recordings classified MVN neurons into two types on the basis of action potential waveforms (Serafin et al., 1991; Johnston et al., 1994), subsequent studies demonstrated that intrinsic physiological properties are distributed continuously rather than bimodally (du Lac and Lisberger, 1995; Bagnall et al., 2007). Many studies have now demonstrated specific differences in intrinsic excitability across MVN neurons with different axonal projections (Sekirnjak and du Lac, 2006; Kolkman et al., 2011a,b) and/or transmitter expression (Takazawa et al., 2004; Bagnall et al., 2007; Shin et al., 2011). The observation that Purkinje cells differentially innervate the somata and dendrites of several types of MVN neurons (Sekirnjak et al., 2003; Shin et al., 2011) complicates analyses of cerebellar learning. Collectively, these findings indicate that the central VOR circuit comprises multiple functionally and molecularly distinct cell types. Given this heterogeneity, how can we target recordings and experimental manipulations to neurons playing distinct roles in behavior and learning?
In this study, we employ a new strategy for classifying heterogeneous neurons and for identifying “marker genes” which can be used in conjunction with modern molecular tools (Dymecki and Kim, 2007; Luo et al., 2008; Knöpfel et al., 2010; Deisseroth, 2011) to label, monitor, and manipulate distinct cell types. In our approach, the quantitative, single-cell expression profiling of genes related to neurotransmitters, ion channels, and candidate molecular markers enables mature neurons to be classified into molecularly distinct types. Anterograde and retrograde tracings provide insights into the functional role of each molecularly defined class of neurons. Application of this approach to the MVN resulted in the identification of six major neuronal classes and associated marker genes. This strategy for classifying neurons and identifying cell-type-specific marker genes is widely applicable to other brain regions and could facilitate the use of molecular genetic tools to examine behavioral roles of distinct neuronal populations.
Materials and Methods
Single-cell sampling and global cDNA amplification.
Experimental protocols for this study were approved by Salk Institute Animal Care and Use Committee. Coronal brainstem slices through the middle part of the rostrocaudal extent of the MVN (bregma from −5.80 to −6.24 mm) were prepared from male juvenile mice of YFP16, GIN, and GlyT2–13 lines in C57BL/6 background. Neurons were acutely dissociated from mouse brainstem slices. Following anesthesia with Nembutal and decapitation, slices (thickness, 250 μm) were prepared by vibratome in ice-cold modified artificial CSF (ACSF) containing the following (in mm): 125 NaCl, 1.25 KCl, 3 MgCl2, 25 NaHCO3, 1 CaCl2, 1.25 NaH2PO4, and 25 dextrose aerated with 95% CO2/5% O2. The central portion of the MVN (including both parvocellular and magnocellular regions) was punched out using a micropipette (660 μm, i.d.), enzymatically dissociated using papain (40 U/ml; Worthington) with cysteine (2 mm) in HEPES-buffered ACSF containing the following (in mm): 140 NaCl, 1.25 KCl, 1.25 NaH2PO4, 10 HEPES, 25 dextrose, 3 MgCl2, 1 CaCl2, 0.0001 tetrodotoxin, 0.02 6-cyano-7-mitroquinoxaline-2,3-dione, 0.05 d-APV at 37°C for 10 min, then triturated using a micropipette (350 μm, o.d.) in 1% BSA in HEPES-buffered ACSF. Dissociated cells were carefully washed one by one to avoid contamination of any tissue debris using glass capillaries (Hempel et al., 2007) and individually transferred to tubes containing cell lysis buffer (Kurimoto et al., 2007). In each sampling batch, we processed eight cells, one negative control (including HEPES-buffered ACSF in the dish used to wash cells), and one positive control (including 10 pg of mouse brain total RNA). Both controls were processed in the same way as the sample containing cells. The negative control was helpful in detecting potential contamination via cell suspension buffer. Each tube also contained spike-in RNAs (Lys, Dap, Phe, and Thr, amounting to 1000, 100, 20, and 5 copies, respectively). The following steps including cell lysis, reverse transcription (RT), exonuclease I treatment, poly(dA) addition, second-strand cDNA synthesis, and the 20 cycle PCR amplification were exactly as described previously (Kurimoto et al., 2007), except for exonuclease I treatment (1.5 U/μl exonuclease I, 37°C for 60 min, followed by 95°C for 30 min). Note that this method amplifies ∼700 bp of the 3′ end of cDNAs. To avoid change in transcriptional state and RNA degradation, tissues and cells were kept on ice whenever possible, and the procedure from decapitation to RT was finished in less than 2 h.
After the initial 20 cycle PCR, we ran an additional 10 cycle PCR to yield a sufficient amount of cDNA fragments for a large number of quantitative PCRs (qPCRs). One microliter of the 20 cycle PCR product was added to 200 μl of reaction mixture of 1× ExTaq buffer containing 0.25 mm each of dATP, dCTP, cGTP, and dTTP, 0.02 μg/μl primer V1(dT)24, 0.02 μg/μl primer V3(dT)24, and 0.05U/μl ExTaq Hot Start Version (Takara Bio), which was divided into four 50 μl scale reactions. The 10 cycle PCR was performed as follows: 95°C for 5 min, 10 cycles of 95°C for 30 s, 67°C for 1 min, and 72°C for 5 min (+ 6 s for each cycle), followed by 72°C for 10 min. The 10 cycle PCR products were mixed together after the reaction, purified using QIAquick PCR purification kit (Qiagen), and eluted into 400 μl of TE buffer. We used the 10 cycle PCR products for qPCR. No obvious distortion in transcript abundance representation was observed after the 10 cycle PCR (see Fig. 1D).
Quality control of single-cell cDNA.
The expression level of a housekeeping gene (Gapdh) was examined by qPCR with the 20 cycle PCR products. We rejected single-cell cDNAs showing an extraordinarily low Gapdh level, which were defined as those scoring qPCR threshold cycle (Ct) greater than mean + SD (4.3%; 9 of 208 cells). The rejected samples were assumed to have suffered severe RNA degradation or cDNA amplification failure. To exclude RNA contamination via the cell-suspending solution during the cell harvest, the sample batches in which Gapdh was detected in the negative control were all rejected (15.4%; 32 of 208 cells). Qualified single-cell cDNAs were further amplified by the 10 cycle PCR.
qPCR was performed using the 7900 Real-Time PCR System (Applied Biosystems) according to the manufacturer's instruction. Reaction mixture (10 μl per well) includes 1 μl of the 10 cycle PCR product, 200 nm of forward and reverse primers, and 100 nm of probe conjugated with 6-FAM at the 5′ end and Iowa Black FQ (Integrated DNA Technologies) or MGBNFQ (Applied Biosystems) at the 3′ end. For the quality control experiment (qPCR for Gapdh), 1 μl of 1/40 dilution of the 20 cycle PCR product was added as a template. Reaction preparation was assisted by the pipetting robot, epMotion (Eppendorf). The sequence of primers and probes are shown in Table 1; genes examined in this study are listed in Table 2. Each reaction plate included a duplicate or triplicate of the standard reaction which was qPCR for Gapdh with cDNA corresponding to 25 ng of mouse brain total RNA. The threshold to determine Ct was set so that Ct values of the standard reaction became identical across all the reaction plates. All reactions were done in duplicate or triplicate, and averaged Ct values were used for further analyses. If Ct values showed a large variation between duplicates (difference >1) or among triplicates (SD > 1), the same reaction was retried until consistent Ct values were obtained. Note that Ct values >30 generally showed larger variation because of the stochastic error associated with PCR and were not subjected to retest.
qPCR primer/probe design.
The target cDNA sequences for which primer/probe sets were designed were selected from 3′-end sequences flanking to the polyadenylation signal identified by polyadq (Tabaska and Zhang, 1999) or PolyA_DB 2 (Lee et al., 2007). In case multiple polyadenylation signals were found, the one flanked by the greatest number of 3′ expressed sequence tag (EST) annotations in UniGene (built 174, http://www.ncbi.nlm.nih.gov/unigene) was chosen. Notably, the qPCR primer/probe for Kcnma1 was designed for the transcript (A830039N20Rik) flanking to the designated 3′ end of Kcnma1 gene in the mouse genome, which was assigned as a part of the KCNMA1 gene in the human genome. This is because the primer/probe designed for the 3′ end of Kcnma1 designated in the mouse genome database (i.e., MGSCv37) detected no transcript in most cells, in striking contrast to the ubiquitously observed BK channel current in MVN neurons (Smith et al., 2002; Gittis and du Lac, 2007; van Welie and du Lac, 2011). We speculate that the recognition of the 3′ end of Kcnma1 in the mouse genome database is not correct, and that A830039N20Rik correctly corresponds to the actual 3′ end, as is the case with the human genome. Primers and probes were designed by Primer Express 3.0 (Applied Biosystems). The specificity of the primer/probe sequence was confirmed by BLAT homology search implemented in Ensembl transcript database (http://www.ensembl.org).
Slice preparation and electrophysiology.
The procedure for obtaining brainstem slices for electrophysiological recordings has been described previously (Bagnall et al., 2007; Shin et al., 2011). Mice were deeply anesthetized with Nembutal and decapitated. Coronal slices (thickness, 250 μm) were prepared by vibratome in ice-cold ACSF containing the following (in mm): 124 NaCl, 5 KCl, 1.3 MgSO4, 26 NaHCO3, 2.5 CaCl2, 1 NaH2PO4, and 11 dextrose aerated with 95% CO2/5% O2, then incubated at 34°C for 30 min, and held at RT for at least 30 min before recording. Kynurenic acid (2 mm), picrotoxin (100 μm), and strychnine (10 μm) were added to the ACSF to block glutamatergic, GABAergic, and glycinergic synaptic transmission during recording, respectively. Whole-cell patch-clamp recording from green fluorescent protein (GFP)-positive cells in the magnocellular region of the middle third of the MVN (in the rostrocaudal plane: −5.85 mm to −6.21 mm from bregma), was performed with glass electrodes (4–8 MΩ) filled with the internal recording solutions containing (in mm): 140 K-gluconate, 8 NaCl, 10 HEPES, 0.1 EGTA, 2 Mg-ATP, and 0.3 Na-GTP (pH 7.2, 280–290 mOsm). Data were analyzed as described by Bagnall et al. (2007).
In situ hybridization.
Double fluorescent in situ hybridization was performed using the protocol used in the Allen Mouse Brain Atlas (Lein et al., 2007), with some modification described below. Experiments were performed in 2- to 6-month-old mice. Two mRNAs were simultaneously detected by riboprobes labeled by digoxigenin (DIG), fluorescein (FLU), or biotin (BIO). The riboprobes were synthesized to cover almost the entire sequence of target mRNA, the specificity of which was assessed by comparing the results with those of the Allen Mouse Brain Atlas (Lein et al., 2007). Labeled riboprobes were sequentially processed with peroxidase-conjugated antibodies against DIG, FLU, or BIO, amplified by tyramide-BIO and tyramide-DIG (Hopman et al., 1998), then visualized by streptavidin conjugated with Alexa 488 (Invitrogen) and anti-DIG antibody conjugated with alkaline phosphatase followed by reaction with HNPP/Fast red TR (Roche Diagnostics), respectively. In situ hybridization with retrograde labeling by fluorogold was performed as described previously (Watakabe et al., 2010). Fluorogold 2.5% (Fluorochrome) in 0.9% NaCl was unilaterally injected into the flocculus by pressure through a craniotomy on the petrosal bone. After 7 d, animals were perfused by 4% paraformaldehyde (PFA) in PBS, and the brainstem was processed for in situ hybridization.
Retrograde tracer injection and immunostaining.
Ten percent biotinylated dextran amine (BDA) (3000 MW; Invitrogen) dissolved in 0.1 m PBS (pH 7.4) was stereotaxically injected into the abducens nucleus, the oculomotor nucleus (OMN), or the MVN of mice (2–6 months old) as described previously (Shin et al., 2011). Rabies virus expressing GFP instead of glycoprotein (Osakada et al., 2011) was unilaterally injected into the flocculus as was fluorogold. Five to seven days after injection, mice were transcardially perfused by 0.1 m PBS, followed by 4% PFA in 0.1 m PBS, under deep anesthesia by Nembutal. Coronal brainstem sections (30 μm) were prepared and subjected to immunostaining and BDA staining. Sections were first blocked and permeabilized by 3% normal donkey serum in 0.3% Triton X-100 in PBS for 1 h at room temperature, then incubated in primary antibody solution including anti-Spp1 (1:400, goat; R & D Systems #AF808) and/or anti-Pcp2 (1:200, rabbit; Abgent #AP6356a) for overnight at 4°C. Immunoreactivity was visualized by anti-goat-DyLight549, anti-goat-DyLight649, and/or anti-rabbit-Cy3 (for all the secondary antibodies, 1:400, raised in donkey; Jackson ImmunoResearch). Subsequently, if needed, BDA was visualized by streptavidin-Alexa Fluor 647 (1 μg/ml; Invitrogen). Distributions of immunoreactivity and retrogradely labeled cells were first mapped by fluorescent images (1297 × 1038 μm) taken with 10× objective lens, which was further confirmed by direct inspection by eye at higher magnification. GFP signal intensity in the GlyT2 line greatly varied across cells; we regarded any trace of GFP signal as GFP positive. Analysis was performed using the Cell Counter plug-in of ImageJ (http://rsbweb.nih.gov/ij/).
Data analysis and statistics.
Ct values are presented in the figures (Figs. 1, 2, 3) and subjected to the hierarchical clustering analysis without normalization. For dendrogram computation with the qPCR results (Figs. 2, 3), the distance (dissimilarity) between leafs (single-cell cDNAs) was defined as “1-Pearson's correlation coefficient of Ct values.” In the cases where the Ct value was greater than average Ct of spike-in RNA Thr (23.07) or not determined (ND), they were assigned Ct values of 40 so that gene expression fluctuation at low levels (<5 copies) did not influence the analyses. All of the dendrograms were calculated by Ward's hierarchical clustering method. In Figure 2A, mutual information (MI) was calculated as follows: where i is the index of the primary dichotomy of the dendrogram (YFP-16-rich cluster/GIN and GlyT2-rich cluster; i = 1, 2), and j corresponds to rounded Ct value (j = 1, …, 40). Multiple pairwise comparison was done by the Wilcoxon rank sum test followed by Benjamini–Hochberg false discovery rate correction. All analyses were done by R (http://www.r-project.org/). Variances were reported in the text and the figures are SDs.
Our strategy for identifying cell types and associated marker genes was to perform single-cell expression profiling of specific sets of genes and to classify neurons based on quantitative gene expression patterns. We chose three distinct classes of genes for analysis: transmitter-related genes, ion channel genes, and candidate marker genes identified in the Allen Mouse Brain Atlas (Lein et al., 2007). Motivated by observations that functionally distinct cell types are partially segregated in different regions of the MVN (Scudder and Fuchs, 1992; Bagnall et al., 2007), we selected genes as marker candidates if they exhibited spatially distinct patterns of hybridization signal (Fig. 1A). Transcripts were analyzed in neurons from the MVN, a brainstem nucleus containing several cell types which subserve the horizontal VOR.
Preparation of single-cell cDNAs from MVN neurons
MVN neurons were harvested from acute brainstem slices obtained from three strains of transgenic mice (YPF-16, Feng et al., 2000; GIN, Oliva et al., 2000; GlyT2, Zeilhofer et al., 2005) which express fluorescent protein in spatially distributed subsets of neurons (Fig. 1B). Fluorescently labeled neurons were isolated with an enzymatic and mechanical dissociation method which preserves neuronal physiology (Gittis and du Lac, 2007). Juvenile mice, age 24–32 d (27.2 ± 2.2), were used to attain a higher yield of healthy neurons; at this stage, performance and plasticity of the horizontal VOR are mature (Faulstich et al., 2004). Neurons sampled primarily from the central portion of the vestibular complex (Fig. 1C) were identified with a fluorescent stereoscope, washed, and individually harvested into PCR tubes by a manual sorting method (Hempel et al., 2007). After RT, single-cell cDNAs were globally amplified using a PCR-based method (Kurimoto et al., 2007), in which the 3′ rapid amplification of cDNA ends method (Tietjen et al., 2003) was modified to attain high amplification of cDNAs with minimal distortion of the transcript abundance distribution (see Materials and Methods). Using this method, a single neuron can yield as much as 50 μg of cDNA (fragments of ∼700 bp from 3′ end), enabling repetition of sufficient numbers of qPCR to ensure reliable, extensive transcript profiling. In total, we obtained cDNAs from 167 individual neurons (YPF-16: n = 61; GIN: n = 52; GlyT2: n = 54) which passed stringent quality control tests (see Materials and Methods).
To evaluate the detection power of the single-cell transcript profiling, the abundance of spike-in control RNAs added to each single-cell sample was examined by qPCR. Lys (1000 copies), Dap (100 copies), and Phe (20 copies) were detected in virtually all the samples (Lys, 167/167; Dap, 164/167; Phe, 164/167), whereas Thr (five copies) was detected in ∼60% of the samples (Thr, 102/167). This indicates that transcripts of >20 copies are nearly always represented in the single-cell transcript profiling. The qPCR results of the spike-in RNAs generally followed the initial copy numbers, supporting the notion that transcript abundance representation in a single cell was maintained after amplification (Fig. 1E).
Expression profiling of neurotransmitter-related genes
To confirm that the methods for dissociation and molecular amplification could reliably evaluate gene expression and to determine neurotransmitter profiles across neurons, single-cell cDNAs were probed with qPCR for five genes related to neurotransmitter synthesis or transport. Genes for vesicular glutamate transporters (VGluT) 1 and/or 2 (Slc17a7 and/or Slc17a6) and glycine transporter 2 (Slc6a5) were used to identify glutamatergic and glycinergic neurons, respectively. GABAergic neurons were identified as those expressing genes for glutamic acid decarboxylase (Gad1 and/or Gad2).
Neurotransmitter-related gene expression in 167 MVN neurons from three lines of transgenic mice is shown as a heat map in Figure 1F. Hierarchical clustering analyses divided the population into two distinct groups, here defined as “excitatory” and “inhibitory” on the basis of mutually exclusive expression of Slc17a6 and Gad1, respectively. All glutamatergic neurons expressed Slc17a6, and a small subset of glutamatergic neurons coexpressed Slc17a7. In contrast, all inhibitory neurons expressed both Gad1 and Gad2. Surprisingly, Slc6a5 was coexpressed in most inhibitory neurons, suggesting that GABA and glycine may be coreleased (Walberg et al., 1990; Wentzel et al., 1993; Tanaka and Ezure, 2004; Lu et al., 2008). Excitatory neurons were obtained exclusively from YFP-16 mice while inhibitory neurons were obtained from all three lines of mice, as is consistent with previous reports indicating that glutamatergic or glycinergic neurons are labeled in the YFP-16 line while GABAergic and glycinergic neurons are labeled in the GIN line (Bagnall et al., 2007).
Expression profiling of ion channel genes
Previous reports demonstrate that functionally distinct MVN neurons differ in their expression of ionic currents and intrinsic excitability. Relative to interneurons targeted in GIN mice, projection neurons targeted in YFP-16 mice have larger repolarizing currents and narrower action potentials which enable faster firing (Bagnall et al., 2007; Gittis and du Lac, 2007; Gittis et al., 2010). To determine whether ion channel genes reflect these differences, we performed quantitative analyses on 36 transcripts encoding α subunits of voltage-gated or calcium-dependent ion channels and on 4 sodium channel β subunits.
Ion channel expression patterns divided the population of 167 MVN neurons into two major groups, each comprising both excitatory and inhibitory neurons, with one group dominated by YFP-16 neurons and the other by GIN and GlyT2 neurons (Fig. 2A). Qualitative and quantitative differences in gene expression across cells were evident. The contribution of each gene to this dichotomy was quantified as MI (see Materials and Methods). Genes related to the initiation and repolarization of action potentials were primarily responsible for the classification into the two groups. Genes with the highest MI values included Na channel α and β subunits (Scn8a, Scn1a, Scn4b, and Scn1b) and repolarizing K channel subunits (Kcna1, Kcnc1). Surprisingly, a “silent” KV α subunit (Kcng4) also ranked high in MI, revealing a potentially important role for Kv6.4 channels in regulating excitability.
Recordings of ionic currents from dissociated MVN neurons demonstrated that Na and KV3, but not BK, currents are significantly larger in YFP-16 than GIN neurons (Gittis and du Lac, 2007). A comparison of whole-cell conductances from that study and transcript levels of the related genes is shown in Figure 2B. Larger Na channel conductance in YFP-16 than GIN neurons is paralleled by significantly higher expression of Na channel genes Scn1a and Scn8a but no difference in Scn2a1 or Scn3a transcript levels. Larger KV3 conductances in YFP-16 versus GIN neurons are matched by significantly higher expression of associated Kcnc1, Kcnc2, and Kcnc3 transcripts. BK channel conductances and corresponding Kcnma1 gene expression were equivalent in YFP-16 and GIN neurons. The parallels between ionic currents measured electrophysiologically and gene expression levels indicates that the relative strength of ionic conductances can be inferred from transcript levels.
Action potential speed and the ability to sustain high firing rates are determined primarily by the amplitude and kinetics of Na and KV3 channels in YFP-16 and GIN neurons (Gittis and du Lac, 2007; Gittis et al., 2010). Expression levels of genes associated with action potentials vary widely in GlyT2 neurons (Fig. 2A), predicting similar variations in action potentials and firing properties. To test this prediction, GlyT2 neurons were targeted for whole-cell recordings in the MVN in acute brainstem slices. As expected from ion channel expression profiling, GlyT2 neurons ranged widely in intrinsic physiology (Fig. 2C). A subset of GlyT2 neurons with narrow action potentials and relatively high maximum firing rates overlapped with the distribution of YFP-16 neurons, while the majority overlapped with GIN neurons (Fig. 2C). The congruence between results from transcript profiling and physiological experiments indicates that ion channel expression profiles can distinguish neurons which exhibit graded differences in intrinsic physiology.
Postinhibitory rebound firing varies widely within vestibular nucleus neurons (Sekirnjak and du Lac, 2002); YFP-16 neurons exhibit significantly more pronounced rebound firing than do GIN neurons (Bagnall et al., 2007). Rebound firing in vestibular and cerebellar nucleus neurons relies predominantly on cationic H channels and T-type Ca channels (Aizenman and Linden, 1999; Sekirnjak and du Lac, 2002; Molineux et al., 2008). Consistent with these findings, transcripts associated with these channels differed significantly between YFP-16 and GIN neurons (Fig. 2A). H channel genes (Hcn1, Hcn2, and Hcn4) were expressed solely or in combination with each other in 88% of YFP-16 neurons but only 57% of GIN neurons (p < 0.001 Fisher exact test; n = 52 and 61, respectively). In neurons expressing H channel genes, transcript levels were significantly higher in YFP-16 than GIN neurons (p < 0.0001; sum of threshold cycle (Ct) values, 90 ± 16 vs 109 ± 12; note that the smaller the Ct is, the more transcript there is). At least one T-type calcium channel gene (Cacna1 g, Cacna1 h, Cacna1i) was expressed in most neurons, but combined T-type channel transcript levels were significantly higher in YFP-16 than GIN neurons (p < 0.0001; sum of Ct values, 90 ± 17 vs 122 ± 11). Together, these results demonstrate that transcript profiling of ion channel genes can reliably differentiate between physiologically distinct cell types.
Neuronal classification incorporating marker gene candidates
Marker gene candidates were identified in the Allen Mouse Brain Atlas (Lein et al., 2007) via systematic screening of the Anatomical Gene Expression Atlas (Ng et al., 2009) and manual queries. Hierarchical clustering analyses incorporating quantitative expression patterns of 14 marker gene candidates together with the 5 neurotransmitter and 40 ion channel-related genes classified neurons into functionally and molecularly distinct groups (Fig. 3). Six major cell types emerged from this classification; three excitatory and three inhibitory (designated E1–3 and I1–3, respectively).
A subset of marker gene candidates identified in the Allen Mouse Brain Atlas (Lein et al., 2007) were expressed predominantly in one cell type, including Adcyap1 (E3), Crh (E2), Crhpb (I3), and Coch (I2). The gene for VGlutT1 (Slc17a7) was also expressed exclusively in E2. Nrn1 was expressed in virtually all glutamatergic neurons, while Spp1 was expressed predominantly in E1 as well as in subsets of E2, I1, and I2. As such, E1 could be defined by the intersection of Spp1 and Nrn1 or Slc17a6. Several marker gene candidates, including neuropeptide genes Sst and Cck, were expressed in multiple cell types.
Two excitatory cell types and one inhibitory cell type (E1, E2, and I1) expressed patterns of ion channel transcripts associated with narrow action potentials and faster firing neurons (e.g., high levels of Scn8a and Kcnc1). Calcium channel transcripts distinguished these cell types; Cacna1e was coexpressed with Cacna1i in E2 but not in E1 or I1 (Fig. 3, E1 vs E2). Inhibitory cell types expressing low levels of narrow action potential-related genes indicative of slower firing neurons could be distinguished by differential expression of Kcnc2 (I2), Cacna1 h (I2), and Kcnf1 (I3) (Fig. 3, I2 vs I3). These results support the principle that molecular classification based on quantitative expression of genes related to neurotransmitters and ion channels can define functionally distinct cell types.
Validation and anatomical distribution of marker genes via in situ hybridization
To confirm marker genes identified with single-cell transcript profiling, we examined their coexpression with neurotransmitter-related genes (excitatory: Slc17a7, Slc17a6; inhibitory: Gad1) using double fluorescent in situ hybridization (DFISH). As predicted, marker genes for excitatory groups E2 (Crh) and E3 (Adycap1) were coexpressed with Slc17a7 (Fig. 4A; 83.0% of Crh-positive neurons, n = 141) and Slc17a6 (Fig. 4B; 81.1% of Adcyap1-positive neurons, n = 74), respectively. Note that Slc17a7 is another marker for E2. Similarly, marker genes for inhibitory groups I2 and I3 (Coch and Crhbp) were expressed exclusively in Gad1-positive neurons (Fig. 4C,D, 98.7% of Coch-positive neurons, n = 75; 97.5% of Crhbp-positive neurons, n = 118).
The anatomical distribution of marker genes provides further evidence that they are expressed in distinct cell types. Adycap1-positive cells (E3) were distributed near the fourth ventricle (Fig. 4B). In contrast, Crh- and Slc17a7-positive neurons (E2) were sparsely distributed in the central portion of the MVN (Fig. 4A) and formed a dense cluster in the supragenual nucleus (data not shown). Of the inhibitory cell markers, Crhbp (I3) was densely expressed near the fourth ventricle (Fig. 4D), dorsal to Coch-positive neurons (I2, Fig. 4C). No cells coexpressed Crhbp and Coch (n = 167 expressing either but not both), confirming that these genes are expressed in distinct populations of neurons.
Functional identity of cell types classified by gene expression profiles
The MVN comprises neurons that differentially project to the cerebellum, ocular motoneurons or locally within the bilateral vestibular system (Highstein and Holstein, 2006). To determine the axonal projections of molecularly defined cell classes, we combined retrograde tracer labeling with in situ hybridization or immunohistochemistry. Injection of fluorogold or glycoprotein-deleted rabies virus expressing enhanced GFP (Osakada et al., 2011) into the cerebellar flocculus labeled neurons in the central portion of the MVN and in the supragenual nucleus which expressed Slc17a7 and/or Crh (Fig. 5A; 21/25 and18/28 fluorogold-positive cells, respectively, n = 2 animals), but not Spp1 (3/47, 0/51, 3/79, and 4/56 rabies virus infected neurons, n = 4 animals). This result indicates that the E2 population includes precerebellar mossy fiber neurons.
To identify premotor cell types, we injected BDA into motor nuclei subserving horizontal eye movements: the abducens and OMN. Experiments were performed in the GlyT2 mouse line, in which GFP-negative neurons are predominantly glutamatergic (Zeilhofer et al., 2005); DFISH confirmed that virtually all MVN neurons expressing Slc6a5 (gene for glycine transporter 2) were GFP positive in the GlyT2 line (data not shown). Neurons retrogradely labeled from the contralateral abducens or OMN were devoid of GFP expression in the GlyT2 line, as expected from previous studies demonstrating that MVN neurons projecting contralaterally to ocular motoneurons are glutamatergic (Wentzel et al., 1995; Shin et al., 2011). Neurons projecting contralaterally to the abducens and OMN were immunopositive for Spp1 (Fig. 5B; for each of two animals, abducens: 64/77 and 88/96 neurons; OMN: 50/52 and 54/69 neurons). These results demonstrate that the E1 population comprises premotor neurons that control movements of the contralateral eye.
Neurons in the MVN that project axons to the ipsilateral abducens nucleus are predominantly glycinergic (Spencer et al., 1989; Shin et al., 2011). A subset of these neurons receive dense somatic innervation from the cerebellar flocculus (Sekirnjak et al., 2003; Shin et al., 2011). We identified such flocculus target neurons in GlyT2 mice with antibodies against the Purkinje cell marker Pcp2 (aka L7) and used immunostaining to determine whether they coexpress Spp1. The majority of glycinergic neurons that were densely surrounded by Pcp2 immunoreactivity were also immunostained for Spp1 (Fig. 5C,D; 40/43 and 69/79 neurons in two animals). These results suggest that flocculus target neurons that mediate cerebellar influence of ipsilateral eye movements are enriched in Spp1-positive neurons included in I1 and/or I2 (Fig. 3).
A summary of major cell types, axonal projections, and marker genes identified in this study is shown in Figure 6. The E1 population comprises glutamatergic neurons that project to contralateral motor nuclei and are identified by coexpression of Spp1 and Slc17a6 (or Nrn1). The E2 population includes glutamatergic precerebellar mossy fiber neurons which can be identified by expression of Crh or Slc17a7. The projections of E3, glutamatergic neurons which express the marker gene Adycap1, have yet to be identified; ion channel transcripts suggest that these neurons have slower firing properties than those in E1 and E2. The I1 population of inhibitory neurons expresses genes related to fast firing; these neurons can be identified by coexpression of Scn4b and inhibitory marker such as Slc6a5, and include cerebellar target neurons projecting to the ipsilateral abducens motor nucleus. I2 and I3 comprise anatomically distinct populations of inhibitory neurons distinguished by marker genes Coch (I2) and Crhbp (I3). The I2 population largely comprises GlyT2 neurons, some of which are likely to be cerebellar target neurons and project ipsilaterally to the abducens (Shin et al., 2011). Previous reports demonstrated that GIN neurons project axons locally or to the contralateral vestibular nucleus (Bagnall et al., 2007; McElvain et al., 2010). The I3 population, which is dominated by GIN neurons and includes transcriptionally similar GlyT2 neurons, is likely to comprise local and commissurally projecting interneurons. Although additional subdivision of these cell types may be prompted by further analyses, each of the previously known classes of vestibular nucleus neurons are represented in our molecularly defined classification scheme. Thus, single-cell transcript profiling was able to successfully classify a heterogeneous population of intermingled brainstem neurons into six anatomically and functionally distinct cell types.
This study presents a new strategy for classifying mature neurons into functionally distinct types and for identifying differentially expressed marker genes via quantitative single-cell transcript analyses. By applying this strategy to a heterogeneous population of MVN neurons, we identified six major cell types together with associated marker genes. In situ hybridization and immunohistochemical analyses in conjunction with retrograde and anterograde tracing confirmed that cell types expressing different marker genes occupy distinct positions within the vestibulo-motor circuits. This approach provides an efficient means of identifying differentially expressed genes which can be used to label and manipulate distinct cell types for functional analyses of neural circuitry and behavior.
Strategic and technical considerations
Several criteria differ across functionally distinct classes of neurons, including dendritic morphology, axonal projections, innervation patterns, intrinsic excitability, and neurotransmitter types. In designing our strategy, we reasoned that genes directly related to excitability and transmitter types have the best chance of distinguishing mature cell types with clear functional differences. To identify marker gene candidates, we used the Allen Mouse Brain Atlas (Lein et al., 2007) and selected 14 genes that appeared to be expressed in subsets of MVN neurons. Although most of these candidates were expressed in several cell types, four were expressed almost exclusively in single populations and thus served as useful markers genes. Additional differentially expressed genes which may further subdivide the major classes defined here could be identified via genome-wide analyses of pooled cDNAs from neurons of each type.
Quantitative analyses of ion channel transcripts serve as a powerful means to classify neurons, even in populations with subtle differences in electrophysiological profiles (Fig. 2). Our findings demonstrate that small variations in action potential width between cell types are evident in expression levels of several ion channel transcripts. In addition, variations in postinhibitory rebound firing within and between cell types is matched by variations in genes encoding channels activated at hyperpolarized membrane potentials (H and T-type Ca channel). Most brain regions comprise cell types which exhibit larger differences in intrinsic excitability than do vestibular nucleus neurons; our strategy should be valuable for identifying subtle differences among cell types even in the absence of prior electrophysiological information.
Although we used juvenile animals (P24–32) to prepare single-cell cDNAs, several lines of evidence indicate that MVN neurons are mature by that stage. First, action potential waveforms and intrinsic excitability attains developmental maturity before P20 (Murphy and du Lac, 2001). Second, performance and plasticity of the horizontal VOR are mature by P21–26 (Faulstich et al., 2004). Finally, in situ hybridization data from 2- to 6-month-old mice matched the gene expression profiles obtained with qPCR in neurons harvested from juveniles.
Single-cell transcript analyses performed with qPCR have several potential sources of variation and error. First, the stochastic nature of gene expression regulation (Elowitz et al., 2002) can result in variability in transcript numbers. This may account for atypical patterns of sparsely expressed transcripts (e.g., Slc17a6 in inhibitory neurons, shown as blue-gray in Fig. 1F). Second, inappropriate design of primers and probes can yield false negatives. We selectively amplified 3′-end sequences of the cDNAs, and, for primer/probe design, used the distribution of 3′ ESTs to identify the most frequently transcribed 3′-end sequence; however, annotation of 3′ ESTs to genes in the database (i.e., UniGene) is not always perfect. In one case (Kcnma1; see Materials and Methods) we redesigned primers and probes that were initially discrepant with previous electrophysiological or anatomical data. Third, any method involving several time-critical steps and biological tissue can be compromised by accumulation of errors. To minimize such errors, we used a robotic pipetting system for qPCR preparation, and performed rigorous quality control experiments, rejecting all samples that might be compromised by contamination or transcript degradation. Mouse lines in which neurons express fluorescent proteins are not necessary for our strategy but can be useful both for distinguishing neurons from debris during cell sorting and for enriching samples with cell types of interest.
Classification of vestibular nucleus neurons into six distinct types
Conjoint expression of transmitter-related, ion channel, and marker candidate genes defined six molecularly distinct classes of vestibular nucleus neurons. Three cell types (E1–3) expressed the glutamate transporter VGluT2 and were defined as excitatory. The remaining three cell types (I1–3) coexpressed genes encoding synthetic enzymes for GABA (Gad65 and Gad67) and for the glycine transporter GlyT2 and were thus defined as inhibitory. Although distinct types of vestibular nucleus neurons are thought to mediate GABAergic versus glycinergic inhibition (Precht et al., 1973; Spencer et al., 1989; Furuya et al., 1991; Camp et al., 2006; Bagnall et al., 2007), the coexpression of these transmitters suggests the possibility of corelease, as has been shown to speed the kinetics of synaptic responses in the auditory brainstem neurons (Lu et al., 2008).
Ion channel transcript levels defined two predominant types of MVN neurons (Fig. 2A), but genes responsible for initiating and repolarizing the action potential (Na and Kv3 channel α subunits) exhibited continuously graded differences across neurons in each type. These results bridge previously discrepant findings (Serafin et al., 1991; Johnston et al., 1994; du Lac and Lisberger, 1995; Straka et al., 2005; Bagnall et al., 2007) by demonstrating both bimodal and continuous patterns of ion channel gene expression. Differences in transmitter expression and axonal projections across neurons within each of the two cell types defined with gene expression pattern underscores the importance of a classification scheme which incorporates multiple functional properties.
Several types of vestibular nucleus neurons have been distinguished on the basis of firing responses during behavior, axonal projections, and/or synaptic inputs (Highstein and Holstein, 2006; Cullen, 2012). Two cell types are thought to directly control the VOR: excitatory “position-vestibular-pause (PVP)” neurons provide strong drive to motoneurons (Scudder and Fuchs, 1992) and flocculus target neurons (FTNs; Lisberger and Pavelko, 1988), subsequently identified as glycinergic (Shin et al., 2011), mediate cerebellar learning. The firing properties of these cell types differ from each other (Scudder and Fuchs, 1992; Lisberger et al. 1994; Ramachandran and Lisberger, 2008), from vestibular nucleus neurons projecting mossy fiber axons to the cerebellar flocculus (Zhang et al., 1995; Cheron et al., 1996), and from other vestibular nucleus neurons (Scudder and Fuchs, 1992, Cullen, 2012). Additional physiological (Precht et al., 1973; Furuya et al., 1991; Malinvaud et al., 2010) and anatomical (Barmack et al., 1993; De Zeeuw et al., 1993; Balaban and Beryozkin, 1994; Holstein et al., 1999a,b; Bagnall et al., 2007; McElvain et al., 2010) studies indicate the presence of several inhibitory cell types with axons projecting locally, commissurally, to the inferior olive, or to ocular motoneurons.
Two of the excitatory cell types defined by our transcript analyses have a clear relationship with vestibular neuronal classes previously identified in vivo. The E1 population, comprising large neurons which project contralaterally to abducens or OMNs, corresponds with PVP interneurons which provide the predominant excitatory drive to motoneurons. (Scudder and Fuchs, 1992; Ramachandran and Lisberger, 2008). The E2 population, which includes precerebellar neurons, corresponds with mossy fiber neurons projecting to the cerebellar flocculus (Zhang et al., 1995; Cheron et al., 1996). The remaining excitatory neurons (E3) have no known function in the VOR but might correspond with vestibulo-sympathetic neurons projecting to autonomic circuits in the medulla (Holstein et al., 2011). The marker gene Adycap1 provides a molecular handle for circuit and behavioral analyses of the role of E3 neurons.
Inhibitory cell types identified in our molecular classification scheme include premotor neurons densely targeted by Purkinje cell terminals (I1) and spatially segregated populations of neurons (I2 and I3) expressing distinct marker genes (Coch and Crhbp, respectively). The I1 population includes FTNs which were identified in vivo as vestibular nucleus neurons receiving powerful, monosynaptic cerebellar inhibition (Lisberger and Pavelko, 1988; Ramachandran and Lisberger, 2008). Glycinergic neurons with sparse somatic innervation by Purkinje cells also project axons to the abducens nucleus (Shin et al., 2011) and are likely to be included in the I2 population which express the premotor marker Spp1. The I3 population predominantly comprises GIN neurons known to restrict axonal projections to the bilateral vestibular complex (McElvain et al., 2010) and are likely to be a predominant source of feedforward inhibition (Biesdorf et al., 2008).
The strategy for classifying neurons and identifying associated marker genes presented here can be applied throughout the nervous system. The extent of neuronal diversity in the cerebellar nuclei (Bagnall et al., 2009; Uusisaari and Knöpfel, 2012), which are anatomical and functionally similar to vestibular nuclei, could be readily addressed with single-cell transcript profiling. Heterogeneity in intrinsic excitability across vestibular ganglion cells (Kalluri et al., 2010; Eatock and Songer, 2011) and anatomical and functional differences between ocular motoneurons (Büttner-Ennever et al., 2001; Davis-López de Carrizosa et al., 2011) suggest that corresponding molecular differences could be identified in VOR sensory inputs and motor outputs. This strategy, which requires no pre-existing knowledge about neuronal diversity, should be particularly effective for fundamental brainstem circuits in which experimental traction was previously limited by an inability to unambiguously identify different cell types.
Cerebellar and vestibular circuits exhibit remarkable plasticity throughout life, and several studies indicate that mechanisms of plasticity differ across neurons (Straka et al., 2005; Dutia, 2010; McElvain et al., 2010; Pettorossi et al., 2011). Future studies exploiting the molecular markers identified here to generate viruses and/or mouse lines for cell-type-specific manipulations will make it possible to analyze cerebellar learning and vestibulo-motor plasticity at the level of gene expression in functionally identified classes of neurons.
This work was supported by the Howard Hughes Medical Institute, and National Institute of Health Grant EY-11027. We thank F. Osakada and E.M. Callaway for the recombinant rabies virus, K. Kurimoto for technical advice, S. Fenstermacher and M. Fuentes for technical assistance, and S. Sahara and P. Slesinger for valuable discussions.
- Correspondence should be addressed to Dr. Sascha du Lac, Howard Hughes Medical Institute, Salk Institute for Biological Studies, La Jolla, California, 92037.