Multi-Scale Molecular Deconstruction of the Serotonin Neuron System
Introduction
Synthesis of the neurotransmitter serotonin (5-hydroxytryptamine, 5HT) broadly defines a class of long-range projection neurons with cell bodies distributed across nine brainstem nuclei—referred to as the raphe system—and axonal projections that collectively innervate most regions of the brain and spinal cord (Azmitia and Gannon, 1986). These neurons regulate a wide array of homeostatic, sensorimotor, cognitive, and behavioral processes, and their dysfunction is likewise associated with a broad range of human disorders. Though often viewed as a single neuron type based on neurochemical (5HT) identity, phenotypic diversity within the 5HT neuron system has been observed at various levels, suggesting the existence of specialized subtypes of 5HT neurons that regulate specific biological functions (Allers and Sharp, 2003; Brust et al., 2014; Calizo et al., 2011; Gaspar and Lillesaar, 2012; Hale and Lowry, 2011; Kiyasova et al., 2011; O’Hearn and Molliver, 1984). However the full extent of molecular, cellular, and functional variation present across the entirety of the 5HT neuronal system, as well as organizing principles, is currently unknown. Here we present progress in these areas, aiming not only to improve understanding of this critical and highly conserved neuronal system, but also ultimately to inform strategies for therapies tailored to particular serotonergic dysfunction.
Historically, subcategorization of 5HT neurons has been based primarily on anatomical location within the raphe system as well as by gross axonal trajectories, as exemplified by common subdivision of the raphe nuclei into rostral versus caudal subgroups that project primarily to forebrain (“ascending pathway”) versus spinal cord (“descending pathway”), respectively. These rostral and caudal anatomical domains can each be further broken down into subdomains (raphe nuclei and subfields), and variation in 5HT neuron morphology, innervation patterns, electrophysiological properties, neurochemistry, and gene expression has been observed both within and between these anatomical domains (Bang et al., 2012; Brust et al., 2014; Calizo et al., 2011; Crawford et al., 2010; Hale and Lowry, 2011; O’Hearn and Molliver, 1984; Spaethling et al., 2014; Wylie et al., 2010). In part, these phenotypic differences likely reflect differences in developmental history, as our earlier rodent fate-mapping studies have shown that rostral and caudal 5HT neurons derive from distinct rhombomeric sublineages in the developing embryo (Jensen et al., 2008). Further, we found that these major rostral and caudal subdivisions are themselves comprised of multiple genetic sublineages of 5HT neurons, which in some cases intermingle within a single anatomical domain. How these differences in 5HT neuron developmental history relate to differences in mature molecular identity and other cellular phenotypes is presently unknown. Furthermore, which parameters best define a 5HT neuron subtype functionally, behaviorally, and clinically is unclear—whether it be mature anatomy, developmental neuron sublineage, molecular expression, or some combination of features. Lacking is a framework for characterizing 5HT neuronal diversity across the entirety of the 5HT system, for identifying underlying principles shaping 5HT neuron subtype organization, and further for linking these neuron subtypes to specific cognitive, behavioral, and/or physiological functions.
Toward this end, we employed a novel approach that iteratively combined intersectional (dual-recombinase) genetic fate mapping (Awatramani et al., 2003; Jensen et al., 2008), manual cell sorting (Hempel et al., 2007; Sugino et al., 2006), and genome-wide RNA-Seq to comprehensively deconstruct the mature mouse 5HT system at multiple levels of granularity. We separately sampled 5HT neurons from each anatomically defined raphe nucleus, and from each embryonically (genetically) defined 5HT neuron sublineage co-populating and intermingling within each nucleus, down to individual 5HT neurons. Using both population and single-neuron scale transcriptome profiles as measures of cell identity, we identify salient principles of 5HT neuron organization and reveal thousands of genes significantly differentially expressed across identified 5HT neuron subtypes. Using in vitro slice electrophysiology to measure intrinsic properties and drug responses, as well as in vivo subtype-specific synaptic suppression and conditional gene knockout approaches, we further show that these molecularly defined 5HT subtypes are also functionally distinct. This resource thus redefines the 5HT neuronal system and offers robust predictions for subtype-specific functions with potential relevance to human disease.
Results
RNA-Seq profiles of 5HT neuron subpopulations sampled separately across anatomical raphe nuclei and cellular sublineages
To query molecular variation across the 5HT neuronal system and how that variation relates to specific 5HT neurons and cellular traits, we iteratively combined intersectional (dual-recombinase) genetic fate mapping (Awatramani et al., 2003; Dymecki et al., 2010; Jensen et al., 2008), manual fluorescent cell sorting (Hempel et al., 2007), and genome-wide RNA sequencing. The first step of fate mapping allowed us to resolve adult 5HT neuron subsets along two informative dimensions: adult anatomical location of cell bodies and originating progenitor cell subset in the embryonic hindbrain, i.e. genetic sublineage (Figure 1A,B, Figure S1). Serotonergic sublineages were identified in the adult brainstem by green fluorescent protein (GFP) expression from the Cre- and Flp-dependent reporter, RC::FrePe (Brust et al., 2014; Dymecki et al., 2010) (see Methods; Figure S1A), which reflects earlier embryonic expression of two marker genes that, in overlap, define the identity and location of the respective parent progenitor cells. One marker gene of the pair—Pet1—is expressed in the CNS exclusively in serotonergic precursors through mature 5HT neurons in adulthood (Hendricks et al., 1999) and is used to drive Flpe recombinase via the transgene Pet1::Flpe (Jensen et al., 2008)(Figure S1A); the other marker gene is expressed in a specific rhombomere in the embryonic hindbrain and used to drive Cre recombinase (Figure S1A). Serotonergic neurons arising from a given rhombomere (R)—R1, R2, R3, R5, or R6 and posterior—are then distinguished in the adult brain as a GFP-labeled 5HT neuron sublineage (reflecting the “intersectional” expression of both marker genes). Serotonergic Pet1::Flpe-expressing neurons falling outside of a given rhombomeric lineage, which we refer to as the “subtractive” population, are marked by mCherry (Figure S1A).
(A) Schematic of mouse embryonic hindbrain (left) with serotonergic primordium marked in dark blue; (center) rhombomeric domains from which 5HT neurons arise (R1, R2, R3, R5, R6P) shaded in different colors (note these same colors are used throughout all figures); (right) their mature anatomical distribution presented in a sagittal brainstem schematic with labeled B1–B9 nuclei. (B) Coronal representations of brainstem sections (adapted from Franklin and Paxinos, 2007) indicating anatomical distribution of lineage-mapped 5HT neurons. Note anatomical mixing of R1, R2, and R3 5HT neuron lineages in the MR (corresponding to B5 and B8 nuclei in (A)). (C) Schematic of the population-scale experimental workflow. (D) Hierarchical clustering of 5HT neuron transcriptomes (using all 13,505 detected genes). (E) Multidimensional scaling using the top 1,000 genes based on log fold change across samples. See also Figure S1 and Supplemental Information for details of clustering analyses.
Distribution of the five major serotonergic developmental sublineages (Figure 1A,B and S1B–H) are as follows: R1-derived 5HT neurons constitute the bulk of the anatomically defined dorsal raphe (DR) serotonergic population and also distribute into the median raphe (MR). We refer to these two 5HT sublineages as R1DR and R1MR, respectively. R2-derived 5HT neurons populate the MR where they are intermixed with R1MR neurons, as well as with a third sublineage, derived from R3. R5-derived 5HT neurons populate parts of the raphe magnus (RMg). Finally, serotonergic neurons derived from R6 and posterior (denoted R6P) are intermixed with R5-derived 5HT neurons in the RMg and also extend caudally, encompassing the raphe obscurus (ROb) and raphe pallidus (RPa). The anatomical distribution of R6P-derived 5HT neurons was determined from the mCherry-expressing subtractive population of 5HT neurons (Figure S1), as, prior to this study, we lacked a Cre driver permitting R6P-specific intersectional (GFP+) labeling. To effectively identify and sort R6P 5HT neurons using GFP, thereby ensuring that all sorted populations expressed the identical fluorescent transgene, we crossed Pet1::Flpe with the Flp-only responsive reporter, RC::Fe (a derivative of RC::FrePe) (Dymecki et al., 2010) and restricted our selection of R6P 5HT neurons to those residing in the ROb and RPa (Figure S1H).
To selectively harvest labeled 5HT neuron subpopulations with high precision, we employed a manual cell sorting approach (Hempel et al., 2007; Sugino et al., 2006) (Figure 1C; see Methods). This method has demonstrated high cell type-specific purity and no overt evidence of expression artifacts due to dissociation (Okaty et al., 2011a, b). Additionally, several of these 5HT neuron subpopulations are small in number, making them less tractable to automated sorting methods such as FACS but ideal for manual cell sorting (Okaty et al., 2011a). Thus, through a combination of genetic sublineage, anatomy, and manual cell sorting, we distinguished and selectively purified six 5HT neuron subsets: R1DR, R1MR, R2, R3, R5, and R6P. For our population-scale transcriptome analyses, pools of on average 55 cells were collected from each subset, with 3 biological replicate pools collected per condition, for a total of 18 sample pools of manually sorted cells spanning the rostral-caudal extent of 5HT neuronal anatomy and developmental sublineage. RNA-Seq libraries were prepared from each pool and sequenced on an Illumina HiSeq2000 or HiSeq2500 platform at an average depth of 30 million reads per sample (see Methods).
Unbiased classification of population-scale 5HT neuron transcriptomes reveals discrete molecular subtypes shaped by both developmental lineage and adult anatomical location
To determine if samples could be grouped into coherent cell types based solely on global gene expression profiles, we applied two different methods of unsupervised clustering: (1) hierarchical clustering (Figure 1D) and (2) multidimensional scaling (MDS; 1E), using TMM-normalized log2(CPM) values as measures of transcript abundance (see Methods; (Anders et al., 2013; Robinson et al., 2010)). As described below, both methods revealed similar relationships among samples, suggesting the existence of robustly defined 5HT neuron molecular subtypes.
Hierarchical clustering (Figure 1D) shows that samples group both by anatomy and lineage. The most generic level of separation is between the medullary (caudal) and midbrain/pontine (rostral) raphe nuclei, as might be expected from gross anatomy; however, a more refined organization is apparent in the dendrogram. Samples within a shared anatomical domain segregate by lineage, indicating that lineage is a more informative classifier than anatomy per se. For example, within the medullary raphe, R5 and R6P 5HT neurons cluster separately, and are nearly as distant from each other as they are from the rostral raphe taken as a whole. Likewise R1MR, R2, and R3 5HT neurons all form distinct clusters despite being anatomically intermixed in the MR, with R2 and R1MR 5HT neurons appearing more similar to each other than to R3 5HT neurons. Notably, R1MR neurons cluster more closely with R2 neurons than to R1DR neurons, suggesting an interplay between lineage and anatomy; i.e. not only does developmental lineage distinguish molecular 5HT neuron subtypes within a given anatomical domain, but in some cases anatomy appears to further shape and distinguish the molecular identities of 5HT neurons derived from a common developmental domain. MDS analysis (Figure 1E) reveals the same overall organizational structure as hierarchical clustering. Projecting the data onto two scaled dimensions shows tight clustering between biological replicates and reveals clear divisions among sublineages and anatomical domains. Thus, at the population scale, 5HT neurons can be classified into six distinct molecular subtypes in an unbiased manner, revealing an organization that bears the imprint of both developmental lineage and mature anatomy.
Unique repertoires of differentially expressed genes distinguish 5HT neuron subtypes and show enrichment for gene categories essential for neuron function
Next we set out to identify genes that are significantly differentially expressed across these identified 5HT neuron subtypes. Using the count-based analysis platform edgeR (Anders et al., 2013; Robinson et al., 2010), we performed “ANOVA-like” assessments, testing for differential expression (DE) between any sample groups, as well as pairwise tests for all pairwise combinations of groups (see Methods). Both analyses revealed thousands of differentially expressed genes at <5% false discovery rate (FDR) (i.e. Benjamini-Hochberg corrected p<0.05), and a minimum fold change of 2 between at least 2 subtypes (Figure 2A–C). Two-way clustering of DE genes identified by the ANOVA-like analysis revealed 5HT neuron subtype-specific expression patterns, as well as complex patterns of gene expression shared by multiple subtypes (Figure 2B). Consistent with the subtype relationships identified by unsupervised clustering methods, pairwise tests showed that the fewest number of DE genes is between R2 and R1MR 5HT neurons, whereas the largest number is between R3 and R6P which both consistently showed the highest number of DE genes in pairwise comparisons with other sample groups (Figure 2C). Bi-clustering of the top 100 most differentially expressed genes is shown in Figure 2D.
(A) Maximum log fold-change (FC) plotted against Benjamini-Hochberg corrected p-values for differential expression between any subtypes (lineage- and anatomy-defined groups) for all detected genes using edgeR (“ANOVA-like” test). The un-shaded upper left quadrant formed by the intersection of the red lines demarcates significantly differentially expressed (DE) genes. Vertical red line corresponds to a corrected p=0.05, thus genes represented by points falling to the left of the line are significantly DE at an FDR<5%, and the horizontal red line indicates logFC=1, corresponding to a two-fold difference in expression between at least two subtypes. (B) Z-normalized clustergram of genes (rows) significantly differentially expressed across sample groups (columns). Leaf node colors in the dendrogram on top of the heatmap correspond to the rhombomere of origin for each sample, as outlined in Figure 1A. (C) Heatmap/table (upper and lower diagonals, respectively) indicating the number of significantly DE genes determined by pairwise contrasts using edgeR. (D) Clustergram of top 100 subtype enriched marker genes as determined by a combination of corrected p-value (<0.01), maximum fold change (>100-fold), mean CPM<1 in at least one sample group (i.e. not detected in at least one other subtype), and low biological replicate variance.
Toward more closely evaluating the functional importance of differentially expressed genes, we applied gene ontology over-expression analysis using the Bioconducter package GOseq (Young et al., 2010). We found that several gene categories critical for neuron function are significantly over-represented among differentially expressed genes; these categories include synaptic transmission, G protein-coupled receptor (GPCR) activity, ion channel activity, cell adhesion, axon guidance, and neuron differentiation. Some overrepresented categories were indicative of highly specific functions, such as sensory perception of pain, multicellular organismal response to stress, and regulation of circadian sleep/wake cycle. Additionally, genes with a direct role in gene transcription were also overrepresented among DE genes and are particularly interesting insofar as they may be involved in 5HT neuron subtype-specific gene regulatory networks. Clustergrams of each of these gene categories reveal distinct patterns of subtype enrichments (Figure 3A–D; an extended heatmap of DE genes encoding co-transmitters and neuropeptide receptors is given in Figure S2). We further validated the expression profile of a subset of these genes at protein or transcript levels, using immunohistology, intersectional genetics, and/or two-color single-molecule fluorescence mRNA in situ hybridization (smFISH) (Figure 4, Figure S3).
(A) G protein-coupled receptors, (B) ion channels and synaptic genes, (C) transcriptional regulators, and (D) cell adhesion and axon guidance molecules. As in Figure 2, displayed expression levels are Z-normalized, and dendrogram colors correspond to rhombomere of origin (lineage) for each sample, as outlined in Figure 1A. Only the top 40 genes are shown for each category. See also Figure S2.
(A) Sagittal brain schematic showing rostral-caudal position (dashed lines) of numbered (1–3) coronal brain schematics shown in B–F. (B) Solid navy circles represent Tac1 5HT neurons (abbreviated 5HT N) as found populating the DR modestly and the ROb, RPa, and lateral paragigantocellular nuclei (LPGC) more substantially. (B′) Confocal photomicrographs showing the density of GFP-marked (green) Tac1 5HT neurons in the DR and ROb in triple transgenic Tac1-IRES2-cre, Pet1::Flpe, RC::FrePe tissue; non-Tac1-cre-expressing 5HT neurons are marked by immunodetection of mCherry (red). (C–F) Confocal photomicrographs showing (C) Pax5 immunodetection (red) in R1DR 5HT neurons (green GFP-marked cells in triple transgenic En1::cre, Pet1::Flpe, RC::FrePe tissue), (D) Nr2f2 immunodetection (red) in R2 5HT neurons (green GFP-marked cells in triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::FrePe tissue), (E) Meis immunodetection (red) in R5 5HT neurons (green GFP-marked cells in triple transgenic Egr2::cre, Pet1::Flpe, RC::FrePe tissue), and (F) Mu opioid receptor (MOR) immunodetection in R5 5HT neurons (green GFP-marked cells in triple transgenic Egr2::cre, Pet1::Flpe, RC::FrePe tissue). (G) Workflow of two-color single-molecule (mRNA) fluorescence in situ hybridization (smFISH), enabling quantitative estimate of relative transcript level in Tph2 mRNA-labeled 5HT neurons across DR, MR, and RMg. (H) Oxtr mRNA per cell area is enriched in 5HT neurons in DR (p<0.0001) and MR (p<0.0001) 5HT neurons versus RMg. (I) Gad1 mRNA is enriched in RMg-residing 5HT neurons compared to 5HT neurons in DR (p<0.0001). (J) Fluorescent photomicrographs showing Oxtr mRNA expression in Tph2+ 5HT neurons in DR (top) and MR (middle), versus RMg (bottom). (K) Fluorescent photomicrographs showing elevated Gad1 mRNA expression in 5HT neurons in RMg (bottom) as compared to DR- and MR-residing 5HT neurons (top and middle, respectively). Individual puncta (red or green) denote single mRNA molecules, while larger yellow spots represent nonspecific autofluorescent lipofuscin particles observed in both channels; see (Lyubimova et al., 2013) for explanation. Outlined (dashed lines) cells indicate Tph2+ 5HT neurons with clearly associated nuclei in the image plane. Black and white panels (middle and right columns) are threshold-transformed masks of the green and red channels, respectively, in which lipofuscin particles, which produce broad-spectrum fluorescent signal, have been removed to permit quantitative analysis of mRNA puncta. Scale bar = 50 μm (B–F), 5 μm (J, K). n.s. = p>0.05; *** = p<0.0001. Bar plots show mean ± SEM. See also Figure S3.
One gene associated with both sensory perception of pain and respiratory control—Tac1, encoding the preprotachykinin protein, a precursor for neuropeptides including neurokinin A and substance P—showed striking enrichment in R6P 5HT neurons (>1000-fold; Figure 2D), with only low transcript levels detected in other 5HT neurons. This suggested that Tac1 might be a suitable marker gene for capturing R6P 5HT neurons intersectionally (a previously unmet need noted earlier). To test this possibility, we generated triple transgenic mice in which the Cre knock-in driver Tac1-IRES2-cre (Zheng, 2013) was partnered with Pet1::Flpe and the RC::FrePe reporter and found that within the 5HT neuron population, Tac1-IRES2-cre expression delineates a subset within the caudal medullary raphe, predominately in RPa and ROb (Figure 4B). A small number of GFP+ Tac1 neurons were also detected in the DR, however the proportion of DR-labeled 5HT neurons was much less than observed in the caudal medullary raphe, consistent with our RNA-Seq findings. Thus our anatomy- and sublineage-specific transcriptomics allowed us to identify a driver line enabling intersectional genetic access to the bulk of R6P 5HT neurons residing in the caudal-most portions of the medullary raphe.
Molecularly distinct 5HT neuron subtypes have distinct cellular properties and regulate distinct functions
Having identified 5HT neuron molecular subtypes through intersectional transcriptomics, we next wanted to test whether these identified subtypes also exhibit unique functional properties, as suggested by overrepresented gene categories among differentially expressed genes. As proof of concept, we focused on the R2 5HT neuron subtype because it enabled testing the function of a particular lineage-defined subtype within a given anatomical domain—the median raphe—as compared to other lineages populating the same region. First, we assayed the intrinsic electrophysiological properties of R2 5HT neurons (intersectional population labeled by GFP, Figure 5A) in comparison with R1MR and R3 5HT neurons (labeled by mCherry) in acute brain slices using whole-cell patch clamp recordings. We found that R2 5HT neurons exhibit a more depolarized resting membrane potential and increased excitability, as demonstrated by their steeper FI curve (Figure 5B–B′).
(A) Confocal photomicrograph (left) showing GFP-labeled R2 5HT neuron soma (green) in MR, with other 5HT neurons expressing mCherry (red) in triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::FrePe mice and threshold-transformed mask (right) showing GFP+ cell soma only. (B) Whole-cell current clamp recordings of MR neurons reveal R2 5HT neurons (GFP+; n=37; −64.3±2.7 mV) exhibit a more depolarized resting membrane potential (p=0.0006) than subtractive 5HT neurons (mCherry+; n=29, −71.7±2.9 mV). (B′) Comparison of current-induced spiking indicates greater excitability of R2 versus other MR 5HT neurons(Phenotype X Current interaction, F(8, 504) = 2.4, p=0.015). Posthoc tests (Fisher’s LSD) indicate R2 5HT neurons spike more than subtractive (mCherry-expressing) neurons at current injections greater than 100 pA (p<0.05). (C) Schematic of intersectional RC::PFtox allele. (C′) Bright-field photomicrograph showing GFP-tox expression in MR of triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::PFtox tissue (right) versus sibling control (left). (C″) Confocal photomicrograph showing enlarged 5HT-labeled axon varicosities in innervation target region (medial septal nucleus) of R2 5HT neurons in triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::PFtox tissue (right) relative to a wild-type (WT) sibling (left); enlarged varicosities are hallmarks of tox-mediated blockade of exocytic release of neurotransmitter. Black and white panels are threshold-transformed masks of 5HT signal. (D) Previously collected data (Kim et al., 2009) reanalyzed and confirmed here show that tox-mediated suppression of neurotransmission from double transgenic Pet1::Flpe, RC::Ftox male mice (n=15; tox expressed in all Pet1::Flpe-expressing 5HT neurons) exhibit enhanced PPI compared to control littermates (n=17) (for a prepulse of 85 dB, p=0.016). (E) Data reanalyzed and confirmed from the same prior study show no discernable effect on PPI upon tox-mediated perturbation of R1 5HT neurons (triple transgenic En1::cre, Pet1::Flpe, RC::PFtox male mice, n=10) as compared to control littermates (n=10) (for a prepulse of 85 dB, p=0.234). (F) By contrast, triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::PFtox male mice (n=10) with tox perturbation of R2 5HT neurons manifest a similar enhancement of PPI as observed upon perturbation of all Pet1::Flpe-expressing 5HT neurons as compared to littermate controls (n=11) (for a prepulse of 85 dB, p=0.045). (D). Scale bar = 50 μm (A,C). n.s. = p>0.05; * = p<0.05; ** = p<0.001. Plots show mean ± SEM. See also Figure S4.
Next, in order to test whether R2 5HT neurons serve specialized behavioral functions relative to other 5HT neuron subtypes, we suppressed neurotransmission from these neurons selectively in vivo and subjected the mice to a battery of behavioral tests, as we had done previously for all 5HT neurons and for R1 5HT neurons specifically (Kim et al., 2009). Here we drove expression of the light chain from tetanus toxin (tox) selectively in R2 5HT neurons to block their exocytic release of neurotransmitters through cleavage of Vamp2—an approach employing the intersectional allele RC::PFtox (Kim et al., 2009), which we previously established as effective (Kim et al., 2009) here partnered with drivers Rse2(Hoxa2)::cre and Pet1::Flpe (Figure 5C–C″). Behavior assays included tests for depressive- and anxiety-like behaviors (grooming test and open field, respectively), social behaviors (three-chamber sociability), and sensorimotor gating (prepulse inhibition (PPI) of the acoustic startle reflex, reflecting an organism’s ability to inhibit neural impulses related to extraneous sensory cues). We found that tox perturbation of R2 5HT neurons produced greater inhibition of the acoustic startle reflex when a startle stimulus (120 dB) is preceded by a quieter prepulse stimulus (85 dB) compared to control littermates (normalized PPI, with a greater value indicating greater inhibition of the startle reflex: 1.48±0.23 versus 1.00±0.36; n= 10 triple transgenics, n= 11 WT, respectively; p=0.045; Figure 5F), suggesting enhanced sensory gating in triple transgenics. There was no effect of genotype (tox-mediated perturbation) on other assayed behaviors (Figure S4A–M; see Supplemental Information, Table S1 summarizing behavior test results), including, notably, baseline startle reflex (Figure S4A). Hence, suppression of exocytic neurotransmission by R2 5HT neurons in mice specifically affects sensorimotor gating, but not other behaviors linked to serotonergic function. Moreover, within the rostral raphe, this function appears to map specifically to R2 5HT neurons, as we have previously shown that silencing all 5HT neurons results in a similar PPI phenotype (Figure 5D, re-analyzed data from (Kim et al., 2009) and re-plotted here as a normalized measure) whereas silencing all R1 5HT neurons has no effect (Figure 5E, re-analyzed data from Kim, et al. 2009 and re-plotted here as a normalized measure). While we did not test the remaining 5HT neuron subtypes (R3, R5, and R6P), the finding that the magnitude of the PPI alteration is roughly equivalent between silencing all 5HT neurons (Figure 5D) versus silencing R2 5HT neurons specifically (Figure 5F) supports the hypothesis that this function may map uniquely to R2 5HT neurons above other 5HT neuron subtypes, however this will need to be proven by further studies.
RNA-Seq profiles at single 5HT neuron resolution recapitulate population-scale sublineage groupings
Having characterized population-scale transcriptomes of 5HT neurons subdivided by lineage and anatomy, we next sorted and performed RNA-Seq on intersectionally marked single cells (Figure 6A) in order to examine patterns of co-expression at the single-neuron level and to investigate the extent of cell-to-cell variability, particularly within a given lineage-defined subtype. We asked whether 5HT neuron subtypes defined by unsupervised clustering of pooled-cell transcriptomes would hold at the single-neuron level, or whether a different organization would emerge. Moreover, we wanted to know whether our defined 5HT neuron subtypes could be further subdivided based on single-neuron gene expression.
(A) Schematic of the single-neuron scale experimental workflow. (B) Pair-wise correlation heatmap of pooled and single-cell transcriptomes. (C) Hierarchical clustering of single 5HT neuron transcriptomes (using all 12,908 detected genes). See Supplemental Information for additional details of clustering analysis. (D) Gene CPM averaged across single 5HT neurons plotted against gene CPM averaged across pools of 5HT neurons for each subtype. Data points are colored based on the percentage of single cells in which the corresponding gene was detected (CPM>1).
In our single-cell libraries, detectable expressed genes (CPM >1) fell from ~12,000 (observed in our pooled samples) to 8,000, as would be expected from a lower complexity library. Notably, pair-wise correlation coefficients between single-cell transcriptomes were lower and more variable than for pooled cells (Figure 6B), consistent with the growing observation of cell-to-cell variability among cells regarded as a somewhat homogenous class. However, despite greater variability, unsupervised clustering of single-cell transcriptomes still robustly resolved subtype identity (Figure 6C), as defined by sublineage and anatomy, and largely recapitulated the subtype groupings found in the analysis of pooled cells (Figure 1D,E). Importantly, averaged single-cell profiles were highly correlated with averaged pooled samples, suggesting that our sampled single cells faithfully represented the larger population (Figure 6D) and serving as further validation of our RNA-Seq data.
Identification of low variance “true” subtype marker genes versus high variance genes indicating further subdivisions
We next explored cell-to-cell variability in transcript expression, which as recently noted (Shalek et al., 2013; Shalek et al., 2014) can be thought of in terms of “analog” and “digital” variation, i.e. differences in relative expression level of a gene across single cells versus on/off expression. We quantified the analog variation as the square of the coefficient of variation of transcript CPM across single neurons, and quantified the digital variation as the percent of neurons in which a given gene was detected at a CPM>1. Both measures of variation were strongly dependent on the population mean. Genes that were highly abundant in pooled-cell RNA-Seq libraries were detected in 80–100% of single neurons, whereas lower abundance genes were generally detected in fewer single-neuron libraries (Figure 6D). This may reflect a combination of technical noise and true biological variation. As testament to the sensitivity of our single-neuron transcriptome profiling, however, in some cases we were able to detect genes that were very low abundance in pooled 5HT neuron libraries (CPM~=2) in 80–100% of single neurons. A recently proposed bootstrap approach to highlighting biological variation over technical variation is to identify genes that vary more than would be expected based solely on their mean expression level (Brennecke et al., 2013), which we implemented for our single-neuron RNA-Seq libraries (Figure 7A). We were particularly interested in identifying cell-to-cell variation of transcripts we found to be subtype enriched in our pooled analysis to assess whether these truly mark all neurons of a given subtype or are expressed by only a subset within.
(A) Smoothed scatterplot (kernel density estimate) of the squared coefficient of variation (CV) plotted against the mean R2 5HT single-neuron expression level for all genes. Solid black line shows regressed fit of the mean-variance relationship; dotted lines indicate the 95% confidence interval; red dots indicate mRNAs varying more than expected based on their mean expression level (implementation of methods outlined in Brennecke, 2013) (B) Heatmap of top 40 R2 5HT neuron-enriched genes as identified by population scale analyses, showing (right to left) percent detection (CPM>1), CV, and expression level, across single R2 5HT neurons, as well as the single and pooled mean expression levels. Red font denotes genes that significantly vary, based on the analysis in (A). (C) Heatmap of R2 5HT single-neuron expression profiles that are highly positively and negatively correlated with Vglut3 expression. The displayed ordering of single neurons was based on clustering of single R2 5HT neuron transcriptomes using all genes, shown in the dendrogram above the heatmap. As in (B), red font denotes significantly variable genes. (D) Comparison of single 5HT neuron RNA-Seq data and smFISH for probes against Tph2 and Vglut3 corroborates existence of Vglut3 high/Tph2 low and Vglut3 low/Tph2 high neuron populations in MR. (E) Scatterplot of Benjamini-Hochberg corrected p-values for differential expression versus fold difference between R2 subgroup 1 and R2 subgroup 2 for all genes using edgeR to identify differentially expressed genes. Red lines indicate cutoffs for significance. (F) Representative current traces of voltage-clamp recordings from R2 5HT neurons in acute slice during bath application of senktide (top current trace) and oxytocin (bottom current trace), agonists for tachykinin 3 and oxytocin receptors, respectively. Bar plot to the right of the current traces shows the percentage of R2 5HT neurons responding to either oxytocin or senktide. Consistent with the RNA-Seq data, these recordings show that R2 5HT neurons respond to either oxytocin or senktide, but not both, and that a greater proportion responds to senktide. See also Figure S5.
For each 5HT neuron subtype we identified a number of “true” marker genes stably expressed across the majority of single neurons assayed (Figure 7B, Figure S5A–E), as well as subtype-enriched genes that varied across single neurons more than would be expected based on their mean expression level (genes in red font Figure 7B, Figure S5A–E). The set of R2 5HT neuron marker genes contained the greatest number of significantly variably expressed genes compared to other subtypes (Figure 7B). Specifically, Vglut3, C1ql3, Pter, Zeb2, and Stxbp6 were all found to be more variable than expected by their mean expression level across single neurons (which was highly correlated with the population-scale mean). If significantly variable genes truly represent markers of more refined subdivisions within a given 5HT neuron subtype, we reasoned that a constellation of genes, or subtype-specific gene regulatory network (GRN), should co-vary with these genes. Thus we computed the pairwise correlation coefficient between all genes, and ranked them by the number of genes for which the absolute value of the pairwise correlation was greater than 0.9 (i.e. highly positively or negatively correlated genes). Vglut3, in addition to being significantly variable, had the highest number of correlated genes by this criterion (Figure 7C), suggesting that Vglut3 is indeed a marker of a finer subdivision within R2 5HT neurons. Notably, Tph2 and Sert (Slc6a4), which were also significantly variable within R2 5HT neurons, were strongly negatively correlated with Vglut3 expression levels, suggesting a subdivision among R2 5HT neurons between a Vglut3 low/Tph2 high subgroup (R2 subtype 1, R2-1) and a Vglut3 high/Tph2 low subgroup (R2 subtype 2, R2-2), as well as a number of other positively and negatively correlated genes (Figure 7C). Of note, both R2 subtypes expressed comparable levels of the serotonergic gene Pet1. Furthermore, the R2 5HT neurons corresponding to R2-1 and R2-2 were found to cluster separately based on their global transcriptomes (Figure 7C), further supporting the notion that these constitute distinct subtypes. The inverse relationship between Vglut3 and Tph2 expression levels in these groups appears to be unique to R2 5HT neurons, a novel finding which was further corroborated by smFISH (Figure 7D, Figure S5F,G). In R1DR neurons, for example, Tph2 and Vglut3 levels are more positively correlated (Figure S5F), whereas in R5 (and possibly R6P) RMg neurons, Vglut3 and Tph2 co-expression are not strongly correlated, and few RMg 5HT neurons express high levels of Vglut3 (Figure S5G).
To identify the full set of differentially expressed genes between R2-1 and R2-2, as with population-scale transcriptomes, we used edgeR to perform gene-by-gene pairwise tests for differential expression. We found that ~300 genes are differentially expressed between these subtypes at an FDR<5% (Figure 7E). Given that single neurons are more variable than pooled-neuron samples, we also applied an alternative test for differential expression designed specifically for single-cell RNA-seq data (Kharchenko et al., 2014), which identified fewer DE genes, however all of which were also identified by edgeR, indicating strong evidence for true differential gene expression between R2-1 and R2-2. To test whether these molecularly distinct subtypes are likewise functionally distinct, we performed whole-cell patch clamp recordings in acute slices from mice in which R2 5HT neurons were intersectionally labeled by GFP and applied ligands for two receptors we found to be differentially expressed between the two subtypes—the oxytocin receptor and the tachykinin 3 receptor. Consistent with our RNA-Seq data, we found that the R2 5HT neurons that respond to oxytocin do not respond to senktide, an agonist for the Tacr3 receptor, and vice versa (Figure 7F). Moreover, the proportion of oxytocin- versus senktide-responding R2 5HT neurons was roughly commensurate with the proportion of Oxtr- and Tacr3-expressing R2 5HT neurons. Thus single-neuron profiling allowed us to identify two molecular subdivisions of R2 5HT neurons that are functionally distinct at the level of neuropeptide response.
5HT neuron subtype-enriched genes are important for subtype-specific behavioral functions
Beyond identifying and characterizing 5HT neuron subtypes based on co-varying gene expression profiles, our transcriptome dataset allows identification of the 5HT neuron expression profile for any gene of interest. For example, we queried whether a clinically relevant gene displaying a distinct expression profile within a subset of 5HT neurons also conferred unique behavioral functions. We chose the Met gene based on the strength of evidence linking hypomorphic gene variants of MET to autism in humans (Campbell et al., 2006; Lambert et al., 2014; Rudie et al., 2012). By RNA-Seq, we detected the highest levels of Met transcript in R2 5HT neurons, with lower levels of expression in R1DR and R1MR and absence in other subtypes. Antibody co-labeling for Met and Tph2 confirmed this profile, as well showing a striking caudal bias to expression of Met protein in the adult DR (Figure 8A–D). Notably, the highest number of Met-expressing cells was found in the MR (Figure 8B). We generated male conditional knockout (CKO) mice that expressed the 5HT neuron-specific ePet1-cre transgene (Scott et al., 2005) and were homozygous for the Metflox allele (Huh et al., 2004), which when recombined by Cre behaves as a null. Deletion of Met in CKO mice was efficient, given that no Met protein was detectable in 5HT neurons as compared to Cre-negative littermate controls (Figure 8C–F″). In a three-chamber sociability test (Figure 8G,K), CKO mice displayed reduced social interest (~35% decreased) compared to sibling controls as measured by time spent in close interaction with an unfamiliar partner mouse (44.7±17.2 s versus 68.7±9.2 s; n=10 CKO, n=15 wild-type (WT), respectively; p=0.016; Figure 8H). Further, CKO mice showed no preference for a chamber containing an unfamiliar mouse compared to an empty chamber (time spent in chamber—89.2±23.1 s (social chamber) versus 60.7±26.6 s (empty chamber); p=0.287; Figure 8I), while WT littermates exhibited a strong preference for the social chamber (110.6±12.6 s (social chamber) versus 40.3±11.9 s (empty chamber); p<0.0001; Figure 8I). Moreover, social approach behaviors were impaired in CKO mice, with latency to approach an unfamiliar mouse being 68.5±34.0 s for CKO mice versus 33.5±11.9 s for sibling controls (p=0.038; Figure 8J). The ability to discern between a familiar and unfamiliar mouse was intact in CKO mice (time spent in close interaction with Stranger 2 (38.6±14.6 s) versus Stranger 1 (19.0±5.9 s), p=0.049; Figure 8L–N), suggesting intact olfaction. However, cumulative time spent in close interaction with any partner mouse (Stranger 1 or Stranger 2) across all test sessions remained significantly decreased in CKO mice compared to controls (102.4±21.6 s versus 140.8±19.0 s, respectively; p=0.017; Figure 8O). By contrast, we observed no effect of genotype on anxiety-related behaviors (Figure 8P), baseline activity (Figure 8Q), or grooming duration (Figure 8R). Thus, ePet1-cre, Metflox/flox CKO mice exhibited deficits in social interest and approach but not in other general behaviors.
(A) Sagittal (left) and coronal (middle, right) brainstem schematics showing location of Met+ 5HT neurons (solid orange circles). Dashed lines indicate rostral-caudal position of numbered (1,2) coronal sections. (B) The majority of Met+ 5HT neurons reside in the MR, with a smaller subset populating the caudal DR (cDR). (C–F″) Confocal photomicrographs showing absence of Met protein specific to CKO tissue (E,F) verus control sibling (C,D) in Tph2+ 5HT neurons in caudal DR (C–C″, E–E″) and MR (D–D″, F–F″). Occasional Met-expressing non-5HT neurons are found in the raphe (see arrows and inset, E). During the social interaction session of the three-chamber sociability test (schematic shown in (G)), (H) CKO mice (n=10) spent less time in close interaction with an unfamiliar mouse relative to littermate controls (n=15) (p=0.016); (I) CKO mice showed no preference for a social chamber relative to an empty chamber (p=0.287), while WT siblings showed a strong preference for the social chamber (p<0.0001); (J) CKO mice versus controls were slower to approach an unfamiliar mouse (p=0.038). During the social recognition session of the three-chamber sociability test (schematic shown in (K)), there was no effect of genotype on (L) time spent in close interaction with either Stranger 1 (p=0.100) or Stranger 2 (p=0.680), (M) time spent in a chamber containing either Stranger 1 (p=0.465) or Stranger 2 (p=0.758), or (N) latency to approach either Stranger 1 (p=0.190) or Stranger 2 (p=0.682). (O) Across all test sessions, cumulative time spent in close interaction with any partner mouse (Stranger 1 or Stranger 2) was significantly reduced in CKO mice compared to control siblings (p=0.017). Relative to WT siblings, CKO mice showed no differences in (P) percent distance traveled in the center of an open field (p=0.823), (Q) total distance traveled in an open field (p=0.326), or (R) grooming duration (p=0.968). Scale bar = 50 μm (C–F). n.s. = p>0.05; * = p<0.05; *** = p<0.0001. Plots show mean ± SEM.
RNA-Seq profiles of 5HT neuron subpopulations sampled separately across anatomical raphe nuclei and cellular sublineages
To query molecular variation across the 5HT neuronal system and how that variation relates to specific 5HT neurons and cellular traits, we iteratively combined intersectional (dual-recombinase) genetic fate mapping (Awatramani et al., 2003; Dymecki et al., 2010; Jensen et al., 2008), manual fluorescent cell sorting (Hempel et al., 2007), and genome-wide RNA sequencing. The first step of fate mapping allowed us to resolve adult 5HT neuron subsets along two informative dimensions: adult anatomical location of cell bodies and originating progenitor cell subset in the embryonic hindbrain, i.e. genetic sublineage (Figure 1A,B, Figure S1). Serotonergic sublineages were identified in the adult brainstem by green fluorescent protein (GFP) expression from the Cre- and Flp-dependent reporter, RC::FrePe (Brust et al., 2014; Dymecki et al., 2010) (see Methods; Figure S1A), which reflects earlier embryonic expression of two marker genes that, in overlap, define the identity and location of the respective parent progenitor cells. One marker gene of the pair—Pet1—is expressed in the CNS exclusively in serotonergic precursors through mature 5HT neurons in adulthood (Hendricks et al., 1999) and is used to drive Flpe recombinase via the transgene Pet1::Flpe (Jensen et al., 2008)(Figure S1A); the other marker gene is expressed in a specific rhombomere in the embryonic hindbrain and used to drive Cre recombinase (Figure S1A). Serotonergic neurons arising from a given rhombomere (R)—R1, R2, R3, R5, or R6 and posterior—are then distinguished in the adult brain as a GFP-labeled 5HT neuron sublineage (reflecting the “intersectional” expression of both marker genes). Serotonergic Pet1::Flpe-expressing neurons falling outside of a given rhombomeric lineage, which we refer to as the “subtractive” population, are marked by mCherry (Figure S1A).
(A) Schematic of mouse embryonic hindbrain (left) with serotonergic primordium marked in dark blue; (center) rhombomeric domains from which 5HT neurons arise (R1, R2, R3, R5, R6P) shaded in different colors (note these same colors are used throughout all figures); (right) their mature anatomical distribution presented in a sagittal brainstem schematic with labeled B1–B9 nuclei. (B) Coronal representations of brainstem sections (adapted from Franklin and Paxinos, 2007) indicating anatomical distribution of lineage-mapped 5HT neurons. Note anatomical mixing of R1, R2, and R3 5HT neuron lineages in the MR (corresponding to B5 and B8 nuclei in (A)). (C) Schematic of the population-scale experimental workflow. (D) Hierarchical clustering of 5HT neuron transcriptomes (using all 13,505 detected genes). (E) Multidimensional scaling using the top 1,000 genes based on log fold change across samples. See also Figure S1 and Supplemental Information for details of clustering analyses.
Distribution of the five major serotonergic developmental sublineages (Figure 1A,B and S1B–H) are as follows: R1-derived 5HT neurons constitute the bulk of the anatomically defined dorsal raphe (DR) serotonergic population and also distribute into the median raphe (MR). We refer to these two 5HT sublineages as R1DR and R1MR, respectively. R2-derived 5HT neurons populate the MR where they are intermixed with R1MR neurons, as well as with a third sublineage, derived from R3. R5-derived 5HT neurons populate parts of the raphe magnus (RMg). Finally, serotonergic neurons derived from R6 and posterior (denoted R6P) are intermixed with R5-derived 5HT neurons in the RMg and also extend caudally, encompassing the raphe obscurus (ROb) and raphe pallidus (RPa). The anatomical distribution of R6P-derived 5HT neurons was determined from the mCherry-expressing subtractive population of 5HT neurons (Figure S1), as, prior to this study, we lacked a Cre driver permitting R6P-specific intersectional (GFP+) labeling. To effectively identify and sort R6P 5HT neurons using GFP, thereby ensuring that all sorted populations expressed the identical fluorescent transgene, we crossed Pet1::Flpe with the Flp-only responsive reporter, RC::Fe (a derivative of RC::FrePe) (Dymecki et al., 2010) and restricted our selection of R6P 5HT neurons to those residing in the ROb and RPa (Figure S1H).
To selectively harvest labeled 5HT neuron subpopulations with high precision, we employed a manual cell sorting approach (Hempel et al., 2007; Sugino et al., 2006) (Figure 1C; see Methods). This method has demonstrated high cell type-specific purity and no overt evidence of expression artifacts due to dissociation (Okaty et al., 2011a, b). Additionally, several of these 5HT neuron subpopulations are small in number, making them less tractable to automated sorting methods such as FACS but ideal for manual cell sorting (Okaty et al., 2011a). Thus, through a combination of genetic sublineage, anatomy, and manual cell sorting, we distinguished and selectively purified six 5HT neuron subsets: R1DR, R1MR, R2, R3, R5, and R6P. For our population-scale transcriptome analyses, pools of on average 55 cells were collected from each subset, with 3 biological replicate pools collected per condition, for a total of 18 sample pools of manually sorted cells spanning the rostral-caudal extent of 5HT neuronal anatomy and developmental sublineage. RNA-Seq libraries were prepared from each pool and sequenced on an Illumina HiSeq2000 or HiSeq2500 platform at an average depth of 30 million reads per sample (see Methods).
Unbiased classification of population-scale 5HT neuron transcriptomes reveals discrete molecular subtypes shaped by both developmental lineage and adult anatomical location
To determine if samples could be grouped into coherent cell types based solely on global gene expression profiles, we applied two different methods of unsupervised clustering: (1) hierarchical clustering (Figure 1D) and (2) multidimensional scaling (MDS; 1E), using TMM-normalized log2(CPM) values as measures of transcript abundance (see Methods; (Anders et al., 2013; Robinson et al., 2010)). As described below, both methods revealed similar relationships among samples, suggesting the existence of robustly defined 5HT neuron molecular subtypes.
Hierarchical clustering (Figure 1D) shows that samples group both by anatomy and lineage. The most generic level of separation is between the medullary (caudal) and midbrain/pontine (rostral) raphe nuclei, as might be expected from gross anatomy; however, a more refined organization is apparent in the dendrogram. Samples within a shared anatomical domain segregate by lineage, indicating that lineage is a more informative classifier than anatomy per se. For example, within the medullary raphe, R5 and R6P 5HT neurons cluster separately, and are nearly as distant from each other as they are from the rostral raphe taken as a whole. Likewise R1MR, R2, and R3 5HT neurons all form distinct clusters despite being anatomically intermixed in the MR, with R2 and R1MR 5HT neurons appearing more similar to each other than to R3 5HT neurons. Notably, R1MR neurons cluster more closely with R2 neurons than to R1DR neurons, suggesting an interplay between lineage and anatomy; i.e. not only does developmental lineage distinguish molecular 5HT neuron subtypes within a given anatomical domain, but in some cases anatomy appears to further shape and distinguish the molecular identities of 5HT neurons derived from a common developmental domain. MDS analysis (Figure 1E) reveals the same overall organizational structure as hierarchical clustering. Projecting the data onto two scaled dimensions shows tight clustering between biological replicates and reveals clear divisions among sublineages and anatomical domains. Thus, at the population scale, 5HT neurons can be classified into six distinct molecular subtypes in an unbiased manner, revealing an organization that bears the imprint of both developmental lineage and mature anatomy.
Unique repertoires of differentially expressed genes distinguish 5HT neuron subtypes and show enrichment for gene categories essential for neuron function
Next we set out to identify genes that are significantly differentially expressed across these identified 5HT neuron subtypes. Using the count-based analysis platform edgeR (Anders et al., 2013; Robinson et al., 2010), we performed “ANOVA-like” assessments, testing for differential expression (DE) between any sample groups, as well as pairwise tests for all pairwise combinations of groups (see Methods). Both analyses revealed thousands of differentially expressed genes at <5% false discovery rate (FDR) (i.e. Benjamini-Hochberg corrected p<0.05), and a minimum fold change of 2 between at least 2 subtypes (Figure 2A–C). Two-way clustering of DE genes identified by the ANOVA-like analysis revealed 5HT neuron subtype-specific expression patterns, as well as complex patterns of gene expression shared by multiple subtypes (Figure 2B). Consistent with the subtype relationships identified by unsupervised clustering methods, pairwise tests showed that the fewest number of DE genes is between R2 and R1MR 5HT neurons, whereas the largest number is between R3 and R6P which both consistently showed the highest number of DE genes in pairwise comparisons with other sample groups (Figure 2C). Bi-clustering of the top 100 most differentially expressed genes is shown in Figure 2D.
(A) Maximum log fold-change (FC) plotted against Benjamini-Hochberg corrected p-values for differential expression between any subtypes (lineage- and anatomy-defined groups) for all detected genes using edgeR (“ANOVA-like” test). The un-shaded upper left quadrant formed by the intersection of the red lines demarcates significantly differentially expressed (DE) genes. Vertical red line corresponds to a corrected p=0.05, thus genes represented by points falling to the left of the line are significantly DE at an FDR<5%, and the horizontal red line indicates logFC=1, corresponding to a two-fold difference in expression between at least two subtypes. (B) Z-normalized clustergram of genes (rows) significantly differentially expressed across sample groups (columns). Leaf node colors in the dendrogram on top of the heatmap correspond to the rhombomere of origin for each sample, as outlined in Figure 1A. (C) Heatmap/table (upper and lower diagonals, respectively) indicating the number of significantly DE genes determined by pairwise contrasts using edgeR. (D) Clustergram of top 100 subtype enriched marker genes as determined by a combination of corrected p-value (<0.01), maximum fold change (>100-fold), mean CPM<1 in at least one sample group (i.e. not detected in at least one other subtype), and low biological replicate variance.
Toward more closely evaluating the functional importance of differentially expressed genes, we applied gene ontology over-expression analysis using the Bioconducter package GOseq (Young et al., 2010). We found that several gene categories critical for neuron function are significantly over-represented among differentially expressed genes; these categories include synaptic transmission, G protein-coupled receptor (GPCR) activity, ion channel activity, cell adhesion, axon guidance, and neuron differentiation. Some overrepresented categories were indicative of highly specific functions, such as sensory perception of pain, multicellular organismal response to stress, and regulation of circadian sleep/wake cycle. Additionally, genes with a direct role in gene transcription were also overrepresented among DE genes and are particularly interesting insofar as they may be involved in 5HT neuron subtype-specific gene regulatory networks. Clustergrams of each of these gene categories reveal distinct patterns of subtype enrichments (Figure 3A–D; an extended heatmap of DE genes encoding co-transmitters and neuropeptide receptors is given in Figure S2). We further validated the expression profile of a subset of these genes at protein or transcript levels, using immunohistology, intersectional genetics, and/or two-color single-molecule fluorescence mRNA in situ hybridization (smFISH) (Figure 4, Figure S3).
(A) G protein-coupled receptors, (B) ion channels and synaptic genes, (C) transcriptional regulators, and (D) cell adhesion and axon guidance molecules. As in Figure 2, displayed expression levels are Z-normalized, and dendrogram colors correspond to rhombomere of origin (lineage) for each sample, as outlined in Figure 1A. Only the top 40 genes are shown for each category. See also Figure S2.
(A) Sagittal brain schematic showing rostral-caudal position (dashed lines) of numbered (1–3) coronal brain schematics shown in B–F. (B) Solid navy circles represent Tac1 5HT neurons (abbreviated 5HT N) as found populating the DR modestly and the ROb, RPa, and lateral paragigantocellular nuclei (LPGC) more substantially. (B′) Confocal photomicrographs showing the density of GFP-marked (green) Tac1 5HT neurons in the DR and ROb in triple transgenic Tac1-IRES2-cre, Pet1::Flpe, RC::FrePe tissue; non-Tac1-cre-expressing 5HT neurons are marked by immunodetection of mCherry (red). (C–F) Confocal photomicrographs showing (C) Pax5 immunodetection (red) in R1DR 5HT neurons (green GFP-marked cells in triple transgenic En1::cre, Pet1::Flpe, RC::FrePe tissue), (D) Nr2f2 immunodetection (red) in R2 5HT neurons (green GFP-marked cells in triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::FrePe tissue), (E) Meis immunodetection (red) in R5 5HT neurons (green GFP-marked cells in triple transgenic Egr2::cre, Pet1::Flpe, RC::FrePe tissue), and (F) Mu opioid receptor (MOR) immunodetection in R5 5HT neurons (green GFP-marked cells in triple transgenic Egr2::cre, Pet1::Flpe, RC::FrePe tissue). (G) Workflow of two-color single-molecule (mRNA) fluorescence in situ hybridization (smFISH), enabling quantitative estimate of relative transcript level in Tph2 mRNA-labeled 5HT neurons across DR, MR, and RMg. (H) Oxtr mRNA per cell area is enriched in 5HT neurons in DR (p<0.0001) and MR (p<0.0001) 5HT neurons versus RMg. (I) Gad1 mRNA is enriched in RMg-residing 5HT neurons compared to 5HT neurons in DR (p<0.0001). (J) Fluorescent photomicrographs showing Oxtr mRNA expression in Tph2+ 5HT neurons in DR (top) and MR (middle), versus RMg (bottom). (K) Fluorescent photomicrographs showing elevated Gad1 mRNA expression in 5HT neurons in RMg (bottom) as compared to DR- and MR-residing 5HT neurons (top and middle, respectively). Individual puncta (red or green) denote single mRNA molecules, while larger yellow spots represent nonspecific autofluorescent lipofuscin particles observed in both channels; see (Lyubimova et al., 2013) for explanation. Outlined (dashed lines) cells indicate Tph2+ 5HT neurons with clearly associated nuclei in the image plane. Black and white panels (middle and right columns) are threshold-transformed masks of the green and red channels, respectively, in which lipofuscin particles, which produce broad-spectrum fluorescent signal, have been removed to permit quantitative analysis of mRNA puncta. Scale bar = 50 μm (B–F), 5 μm (J, K). n.s. = p>0.05; *** = p<0.0001. Bar plots show mean ± SEM. See also Figure S3.
One gene associated with both sensory perception of pain and respiratory control—Tac1, encoding the preprotachykinin protein, a precursor for neuropeptides including neurokinin A and substance P—showed striking enrichment in R6P 5HT neurons (>1000-fold; Figure 2D), with only low transcript levels detected in other 5HT neurons. This suggested that Tac1 might be a suitable marker gene for capturing R6P 5HT neurons intersectionally (a previously unmet need noted earlier). To test this possibility, we generated triple transgenic mice in which the Cre knock-in driver Tac1-IRES2-cre (Zheng, 2013) was partnered with Pet1::Flpe and the RC::FrePe reporter and found that within the 5HT neuron population, Tac1-IRES2-cre expression delineates a subset within the caudal medullary raphe, predominately in RPa and ROb (Figure 4B). A small number of GFP+ Tac1 neurons were also detected in the DR, however the proportion of DR-labeled 5HT neurons was much less than observed in the caudal medullary raphe, consistent with our RNA-Seq findings. Thus our anatomy- and sublineage-specific transcriptomics allowed us to identify a driver line enabling intersectional genetic access to the bulk of R6P 5HT neurons residing in the caudal-most portions of the medullary raphe.
Molecularly distinct 5HT neuron subtypes have distinct cellular properties and regulate distinct functions
Having identified 5HT neuron molecular subtypes through intersectional transcriptomics, we next wanted to test whether these identified subtypes also exhibit unique functional properties, as suggested by overrepresented gene categories among differentially expressed genes. As proof of concept, we focused on the R2 5HT neuron subtype because it enabled testing the function of a particular lineage-defined subtype within a given anatomical domain—the median raphe—as compared to other lineages populating the same region. First, we assayed the intrinsic electrophysiological properties of R2 5HT neurons (intersectional population labeled by GFP, Figure 5A) in comparison with R1MR and R3 5HT neurons (labeled by mCherry) in acute brain slices using whole-cell patch clamp recordings. We found that R2 5HT neurons exhibit a more depolarized resting membrane potential and increased excitability, as demonstrated by their steeper FI curve (Figure 5B–B′).
(A) Confocal photomicrograph (left) showing GFP-labeled R2 5HT neuron soma (green) in MR, with other 5HT neurons expressing mCherry (red) in triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::FrePe mice and threshold-transformed mask (right) showing GFP+ cell soma only. (B) Whole-cell current clamp recordings of MR neurons reveal R2 5HT neurons (GFP+; n=37; −64.3±2.7 mV) exhibit a more depolarized resting membrane potential (p=0.0006) than subtractive 5HT neurons (mCherry+; n=29, −71.7±2.9 mV). (B′) Comparison of current-induced spiking indicates greater excitability of R2 versus other MR 5HT neurons(Phenotype X Current interaction, F(8, 504) = 2.4, p=0.015). Posthoc tests (Fisher’s LSD) indicate R2 5HT neurons spike more than subtractive (mCherry-expressing) neurons at current injections greater than 100 pA (p<0.05). (C) Schematic of intersectional RC::PFtox allele. (C′) Bright-field photomicrograph showing GFP-tox expression in MR of triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::PFtox tissue (right) versus sibling control (left). (C″) Confocal photomicrograph showing enlarged 5HT-labeled axon varicosities in innervation target region (medial septal nucleus) of R2 5HT neurons in triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::PFtox tissue (right) relative to a wild-type (WT) sibling (left); enlarged varicosities are hallmarks of tox-mediated blockade of exocytic release of neurotransmitter. Black and white panels are threshold-transformed masks of 5HT signal. (D) Previously collected data (Kim et al., 2009) reanalyzed and confirmed here show that tox-mediated suppression of neurotransmission from double transgenic Pet1::Flpe, RC::Ftox male mice (n=15; tox expressed in all Pet1::Flpe-expressing 5HT neurons) exhibit enhanced PPI compared to control littermates (n=17) (for a prepulse of 85 dB, p=0.016). (E) Data reanalyzed and confirmed from the same prior study show no discernable effect on PPI upon tox-mediated perturbation of R1 5HT neurons (triple transgenic En1::cre, Pet1::Flpe, RC::PFtox male mice, n=10) as compared to control littermates (n=10) (for a prepulse of 85 dB, p=0.234). (F) By contrast, triple transgenic Rse2(Hoxa2)::cre, Pet1::Flpe, RC::PFtox male mice (n=10) with tox perturbation of R2 5HT neurons manifest a similar enhancement of PPI as observed upon perturbation of all Pet1::Flpe-expressing 5HT neurons as compared to littermate controls (n=11) (for a prepulse of 85 dB, p=0.045). (D). Scale bar = 50 μm (A,C). n.s. = p>0.05; * = p<0.05; ** = p<0.001. Plots show mean ± SEM. See also Figure S4.
Next, in order to test whether R2 5HT neurons serve specialized behavioral functions relative to other 5HT neuron subtypes, we suppressed neurotransmission from these neurons selectively in vivo and subjected the mice to a battery of behavioral tests, as we had done previously for all 5HT neurons and for R1 5HT neurons specifically (Kim et al., 2009). Here we drove expression of the light chain from tetanus toxin (tox) selectively in R2 5HT neurons to block their exocytic release of neurotransmitters through cleavage of Vamp2—an approach employing the intersectional allele RC::PFtox (Kim et al., 2009), which we previously established as effective (Kim et al., 2009) here partnered with drivers Rse2(Hoxa2)::cre and Pet1::Flpe (Figure 5C–C″). Behavior assays included tests for depressive- and anxiety-like behaviors (grooming test and open field, respectively), social behaviors (three-chamber sociability), and sensorimotor gating (prepulse inhibition (PPI) of the acoustic startle reflex, reflecting an organism’s ability to inhibit neural impulses related to extraneous sensory cues). We found that tox perturbation of R2 5HT neurons produced greater inhibition of the acoustic startle reflex when a startle stimulus (120 dB) is preceded by a quieter prepulse stimulus (85 dB) compared to control littermates (normalized PPI, with a greater value indicating greater inhibition of the startle reflex: 1.48±0.23 versus 1.00±0.36; n= 10 triple transgenics, n= 11 WT, respectively; p=0.045; Figure 5F), suggesting enhanced sensory gating in triple transgenics. There was no effect of genotype (tox-mediated perturbation) on other assayed behaviors (Figure S4A–M; see Supplemental Information, Table S1 summarizing behavior test results), including, notably, baseline startle reflex (Figure S4A). Hence, suppression of exocytic neurotransmission by R2 5HT neurons in mice specifically affects sensorimotor gating, but not other behaviors linked to serotonergic function. Moreover, within the rostral raphe, this function appears to map specifically to R2 5HT neurons, as we have previously shown that silencing all 5HT neurons results in a similar PPI phenotype (Figure 5D, re-analyzed data from (Kim et al., 2009) and re-plotted here as a normalized measure) whereas silencing all R1 5HT neurons has no effect (Figure 5E, re-analyzed data from Kim, et al. 2009 and re-plotted here as a normalized measure). While we did not test the remaining 5HT neuron subtypes (R3, R5, and R6P), the finding that the magnitude of the PPI alteration is roughly equivalent between silencing all 5HT neurons (Figure 5D) versus silencing R2 5HT neurons specifically (Figure 5F) supports the hypothesis that this function may map uniquely to R2 5HT neurons above other 5HT neuron subtypes, however this will need to be proven by further studies.
RNA-Seq profiles at single 5HT neuron resolution recapitulate population-scale sublineage groupings
Having characterized population-scale transcriptomes of 5HT neurons subdivided by lineage and anatomy, we next sorted and performed RNA-Seq on intersectionally marked single cells (Figure 6A) in order to examine patterns of co-expression at the single-neuron level and to investigate the extent of cell-to-cell variability, particularly within a given lineage-defined subtype. We asked whether 5HT neuron subtypes defined by unsupervised clustering of pooled-cell transcriptomes would hold at the single-neuron level, or whether a different organization would emerge. Moreover, we wanted to know whether our defined 5HT neuron subtypes could be further subdivided based on single-neuron gene expression.
(A) Schematic of the single-neuron scale experimental workflow. (B) Pair-wise correlation heatmap of pooled and single-cell transcriptomes. (C) Hierarchical clustering of single 5HT neuron transcriptomes (using all 12,908 detected genes). See Supplemental Information for additional details of clustering analysis. (D) Gene CPM averaged across single 5HT neurons plotted against gene CPM averaged across pools of 5HT neurons for each subtype. Data points are colored based on the percentage of single cells in which the corresponding gene was detected (CPM>1).
In our single-cell libraries, detectable expressed genes (CPM >1) fell from ~12,000 (observed in our pooled samples) to 8,000, as would be expected from a lower complexity library. Notably, pair-wise correlation coefficients between single-cell transcriptomes were lower and more variable than for pooled cells (Figure 6B), consistent with the growing observation of cell-to-cell variability among cells regarded as a somewhat homogenous class. However, despite greater variability, unsupervised clustering of single-cell transcriptomes still robustly resolved subtype identity (Figure 6C), as defined by sublineage and anatomy, and largely recapitulated the subtype groupings found in the analysis of pooled cells (Figure 1D,E). Importantly, averaged single-cell profiles were highly correlated with averaged pooled samples, suggesting that our sampled single cells faithfully represented the larger population (Figure 6D) and serving as further validation of our RNA-Seq data.
Identification of low variance “true” subtype marker genes versus high variance genes indicating further subdivisions
We next explored cell-to-cell variability in transcript expression, which as recently noted (Shalek et al., 2013; Shalek et al., 2014) can be thought of in terms of “analog” and “digital” variation, i.e. differences in relative expression level of a gene across single cells versus on/off expression. We quantified the analog variation as the square of the coefficient of variation of transcript CPM across single neurons, and quantified the digital variation as the percent of neurons in which a given gene was detected at a CPM>1. Both measures of variation were strongly dependent on the population mean. Genes that were highly abundant in pooled-cell RNA-Seq libraries were detected in 80–100% of single neurons, whereas lower abundance genes were generally detected in fewer single-neuron libraries (Figure 6D). This may reflect a combination of technical noise and true biological variation. As testament to the sensitivity of our single-neuron transcriptome profiling, however, in some cases we were able to detect genes that were very low abundance in pooled 5HT neuron libraries (CPM~=2) in 80–100% of single neurons. A recently proposed bootstrap approach to highlighting biological variation over technical variation is to identify genes that vary more than would be expected based solely on their mean expression level (Brennecke et al., 2013), which we implemented for our single-neuron RNA-Seq libraries (Figure 7A). We were particularly interested in identifying cell-to-cell variation of transcripts we found to be subtype enriched in our pooled analysis to assess whether these truly mark all neurons of a given subtype or are expressed by only a subset within.
(A) Smoothed scatterplot (kernel density estimate) of the squared coefficient of variation (CV) plotted against the mean R2 5HT single-neuron expression level for all genes. Solid black line shows regressed fit of the mean-variance relationship; dotted lines indicate the 95% confidence interval; red dots indicate mRNAs varying more than expected based on their mean expression level (implementation of methods outlined in Brennecke, 2013) (B) Heatmap of top 40 R2 5HT neuron-enriched genes as identified by population scale analyses, showing (right to left) percent detection (CPM>1), CV, and expression level, across single R2 5HT neurons, as well as the single and pooled mean expression levels. Red font denotes genes that significantly vary, based on the analysis in (A). (C) Heatmap of R2 5HT single-neuron expression profiles that are highly positively and negatively correlated with Vglut3 expression. The displayed ordering of single neurons was based on clustering of single R2 5HT neuron transcriptomes using all genes, shown in the dendrogram above the heatmap. As in (B), red font denotes significantly variable genes. (D) Comparison of single 5HT neuron RNA-Seq data and smFISH for probes against Tph2 and Vglut3 corroborates existence of Vglut3 high/Tph2 low and Vglut3 low/Tph2 high neuron populations in MR. (E) Scatterplot of Benjamini-Hochberg corrected p-values for differential expression versus fold difference between R2 subgroup 1 and R2 subgroup 2 for all genes using edgeR to identify differentially expressed genes. Red lines indicate cutoffs for significance. (F) Representative current traces of voltage-clamp recordings from R2 5HT neurons in acute slice during bath application of senktide (top current trace) and oxytocin (bottom current trace), agonists for tachykinin 3 and oxytocin receptors, respectively. Bar plot to the right of the current traces shows the percentage of R2 5HT neurons responding to either oxytocin or senktide. Consistent with the RNA-Seq data, these recordings show that R2 5HT neurons respond to either oxytocin or senktide, but not both, and that a greater proportion responds to senktide. See also Figure S5.
For each 5HT neuron subtype we identified a number of “true” marker genes stably expressed across the majority of single neurons assayed (Figure 7B, Figure S5A–E), as well as subtype-enriched genes that varied across single neurons more than would be expected based on their mean expression level (genes in red font Figure 7B, Figure S5A–E). The set of R2 5HT neuron marker genes contained the greatest number of significantly variably expressed genes compared to other subtypes (Figure 7B). Specifically, Vglut3, C1ql3, Pter, Zeb2, and Stxbp6 were all found to be more variable than expected by their mean expression level across single neurons (which was highly correlated with the population-scale mean). If significantly variable genes truly represent markers of more refined subdivisions within a given 5HT neuron subtype, we reasoned that a constellation of genes, or subtype-specific gene regulatory network (GRN), should co-vary with these genes. Thus we computed the pairwise correlation coefficient between all genes, and ranked them by the number of genes for which the absolute value of the pairwise correlation was greater than 0.9 (i.e. highly positively or negatively correlated genes). Vglut3, in addition to being significantly variable, had the highest number of correlated genes by this criterion (Figure 7C), suggesting that Vglut3 is indeed a marker of a finer subdivision within R2 5HT neurons. Notably, Tph2 and Sert (Slc6a4), which were also significantly variable within R2 5HT neurons, were strongly negatively correlated with Vglut3 expression levels, suggesting a subdivision among R2 5HT neurons between a Vglut3 low/Tph2 high subgroup (R2 subtype 1, R2-1) and a Vglut3 high/Tph2 low subgroup (R2 subtype 2, R2-2), as well as a number of other positively and negatively correlated genes (Figure 7C). Of note, both R2 subtypes expressed comparable levels of the serotonergic gene Pet1. Furthermore, the R2 5HT neurons corresponding to R2-1 and R2-2 were found to cluster separately based on their global transcriptomes (Figure 7C), further supporting the notion that these constitute distinct subtypes. The inverse relationship between Vglut3 and Tph2 expression levels in these groups appears to be unique to R2 5HT neurons, a novel finding which was further corroborated by smFISH (Figure 7D, Figure S5F,G). In R1DR neurons, for example, Tph2 and Vglut3 levels are more positively correlated (Figure S5F), whereas in R5 (and possibly R6P) RMg neurons, Vglut3 and Tph2 co-expression are not strongly correlated, and few RMg 5HT neurons express high levels of Vglut3 (Figure S5G).
To identify the full set of differentially expressed genes between R2-1 and R2-2, as with population-scale transcriptomes, we used edgeR to perform gene-by-gene pairwise tests for differential expression. We found that ~300 genes are differentially expressed between these subtypes at an FDR<5% (Figure 7E). Given that single neurons are more variable than pooled-neuron samples, we also applied an alternative test for differential expression designed specifically for single-cell RNA-seq data (Kharchenko et al., 2014), which identified fewer DE genes, however all of which were also identified by edgeR, indicating strong evidence for true differential gene expression between R2-1 and R2-2. To test whether these molecularly distinct subtypes are likewise functionally distinct, we performed whole-cell patch clamp recordings in acute slices from mice in which R2 5HT neurons were intersectionally labeled by GFP and applied ligands for two receptors we found to be differentially expressed between the two subtypes—the oxytocin receptor and the tachykinin 3 receptor. Consistent with our RNA-Seq data, we found that the R2 5HT neurons that respond to oxytocin do not respond to senktide, an agonist for the Tacr3 receptor, and vice versa (Figure 7F). Moreover, the proportion of oxytocin- versus senktide-responding R2 5HT neurons was roughly commensurate with the proportion of Oxtr- and Tacr3-expressing R2 5HT neurons. Thus single-neuron profiling allowed us to identify two molecular subdivisions of R2 5HT neurons that are functionally distinct at the level of neuropeptide response.
5HT neuron subtype-enriched genes are important for subtype-specific behavioral functions
Beyond identifying and characterizing 5HT neuron subtypes based on co-varying gene expression profiles, our transcriptome dataset allows identification of the 5HT neuron expression profile for any gene of interest. For example, we queried whether a clinically relevant gene displaying a distinct expression profile within a subset of 5HT neurons also conferred unique behavioral functions. We chose the Met gene based on the strength of evidence linking hypomorphic gene variants of MET to autism in humans (Campbell et al., 2006; Lambert et al., 2014; Rudie et al., 2012). By RNA-Seq, we detected the highest levels of Met transcript in R2 5HT neurons, with lower levels of expression in R1DR and R1MR and absence in other subtypes. Antibody co-labeling for Met and Tph2 confirmed this profile, as well showing a striking caudal bias to expression of Met protein in the adult DR (Figure 8A–D). Notably, the highest number of Met-expressing cells was found in the MR (Figure 8B). We generated male conditional knockout (CKO) mice that expressed the 5HT neuron-specific ePet1-cre transgene (Scott et al., 2005) and were homozygous for the Metflox allele (Huh et al., 2004), which when recombined by Cre behaves as a null. Deletion of Met in CKO mice was efficient, given that no Met protein was detectable in 5HT neurons as compared to Cre-negative littermate controls (Figure 8C–F″). In a three-chamber sociability test (Figure 8G,K), CKO mice displayed reduced social interest (~35% decreased) compared to sibling controls as measured by time spent in close interaction with an unfamiliar partner mouse (44.7±17.2 s versus 68.7±9.2 s; n=10 CKO, n=15 wild-type (WT), respectively; p=0.016; Figure 8H). Further, CKO mice showed no preference for a chamber containing an unfamiliar mouse compared to an empty chamber (time spent in chamber—89.2±23.1 s (social chamber) versus 60.7±26.6 s (empty chamber); p=0.287; Figure 8I), while WT littermates exhibited a strong preference for the social chamber (110.6±12.6 s (social chamber) versus 40.3±11.9 s (empty chamber); p<0.0001; Figure 8I). Moreover, social approach behaviors were impaired in CKO mice, with latency to approach an unfamiliar mouse being 68.5±34.0 s for CKO mice versus 33.5±11.9 s for sibling controls (p=0.038; Figure 8J). The ability to discern between a familiar and unfamiliar mouse was intact in CKO mice (time spent in close interaction with Stranger 2 (38.6±14.6 s) versus Stranger 1 (19.0±5.9 s), p=0.049; Figure 8L–N), suggesting intact olfaction. However, cumulative time spent in close interaction with any partner mouse (Stranger 1 or Stranger 2) across all test sessions remained significantly decreased in CKO mice compared to controls (102.4±21.6 s versus 140.8±19.0 s, respectively; p=0.017; Figure 8O). By contrast, we observed no effect of genotype on anxiety-related behaviors (Figure 8P), baseline activity (Figure 8Q), or grooming duration (Figure 8R). Thus, ePet1-cre, Metflox/flox CKO mice exhibited deficits in social interest and approach but not in other general behaviors.
(A) Sagittal (left) and coronal (middle, right) brainstem schematics showing location of Met+ 5HT neurons (solid orange circles). Dashed lines indicate rostral-caudal position of numbered (1,2) coronal sections. (B) The majority of Met+ 5HT neurons reside in the MR, with a smaller subset populating the caudal DR (cDR). (C–F″) Confocal photomicrographs showing absence of Met protein specific to CKO tissue (E,F) verus control sibling (C,D) in Tph2+ 5HT neurons in caudal DR (C–C″, E–E″) and MR (D–D″, F–F″). Occasional Met-expressing non-5HT neurons are found in the raphe (see arrows and inset, E). During the social interaction session of the three-chamber sociability test (schematic shown in (G)), (H) CKO mice (n=10) spent less time in close interaction with an unfamiliar mouse relative to littermate controls (n=15) (p=0.016); (I) CKO mice showed no preference for a social chamber relative to an empty chamber (p=0.287), while WT siblings showed a strong preference for the social chamber (p<0.0001); (J) CKO mice versus controls were slower to approach an unfamiliar mouse (p=0.038). During the social recognition session of the three-chamber sociability test (schematic shown in (K)), there was no effect of genotype on (L) time spent in close interaction with either Stranger 1 (p=0.100) or Stranger 2 (p=0.680), (M) time spent in a chamber containing either Stranger 1 (p=0.465) or Stranger 2 (p=0.758), or (N) latency to approach either Stranger 1 (p=0.190) or Stranger 2 (p=0.682). (O) Across all test sessions, cumulative time spent in close interaction with any partner mouse (Stranger 1 or Stranger 2) was significantly reduced in CKO mice compared to control siblings (p=0.017). Relative to WT siblings, CKO mice showed no differences in (P) percent distance traveled in the center of an open field (p=0.823), (Q) total distance traveled in an open field (p=0.326), or (R) grooming duration (p=0.968). Scale bar = 50 μm (C–F). n.s. = p>0.05; * = p<0.05; *** = p<0.0001. Plots show mean ± SEM.
Discussion
The mammalian brain, by virtue of its complexity, poses an extreme case of the cartographer’s dilemma: at what level of abstraction does a map of the brain render its structure and function intelligible without sacrificing essential detail? In part, this is a question of “granularity”—what are the most informative functional units of which the brain is comprised and at what scales do they emerge? One approach to this question has been to reduce the brain to a composite of discrete cell types; i.e. groupings of cells that share a common set of properties that determine their function (Fishell and Heintz, 2013; Masland, 2004; Nelson et al., 2006; Sugino et al., 2006). Establishing robust criteria for defining neuronal cell types, however, poses yet another conundrum, particularly given the growing observation of cell-to-cell variation within what have previously been regarded as largely homogeneous cell types—if cell types are defined too broadly, functionally meaningful diversity may be obscured, however if defined too narrowly, the map may tend toward inscrutable complexity, rendering its organizing principles less apparent. Here we address aspects of these challenges for a particular group of brainstem neurons defined by their shared capacity to synthesize 5HT, but for which growing awareness of their diversity has hinted at the existence of a more elaborate and functionally relevant subtype organization.
Employing a novel approach that iteratively combined intersectional genetic fate mapping, anatomical microdissection, cell sorting, and RNA-Seq, we measured transcriptomes for an extensive sampling of serotonergic neurons across the 5HT system, at both population and single-cell levels. This allowed for a multi-scale molecular deconstruction of the 5HT system—i.e. gene expression profiles at system-wide, subanatomical, sublineal, and single-neuron levels of granularity—and thereby serves as a platform for directly quantifying how much underlying diversity is masked or uncovered when viewing the system at different levels. Further, it allows for identifying which variables best explain 5HT neuron variation and therefore provide a coherent framework for 5HT system organization.
Through analysis of this resource, we demonstrate a previously underappreciated degree of molecular variation across 5HT neurons present at all scales; i.e. profound differences in gene expression both between and within anatomical domains and 5HT neuron sublineages. Unsupervised clustering of 5HT neuron transcriptomes revealed that this molecular variation displays an organization suggestive of distinct 5HT neuron subtypes defined by a combination of lineage and anatomy. Within all anatomical domains examined, 5HT neuron transcriptomes cluster by rhombomeric lineage, however in some cases 5HT neuron transcriptomes sampled from a common lineage segregate by anatomy, suggesting a role for both developmental lineage and anatomy in shaping 5HT cell identity. This classification scheme is largely robust even at the scale of single 5HT neuron transcriptomes, despite cell-to-cell variation that nonetheless appears bounded by lineal subtype rather than varying without constraint. This might not have been the case; for example, in theory different rhombomeres could give rise to different proportions of common core 5HT neuron subtypes rather than to different subtypes, per se. In this case, neurons arising from each domain would appear distinct in the pooled-cell analysis, insofar as they were comprised of different proportions of underlying subtypes, however at the single-cell level these progenitor domain- and anatomy-defined groups would break down, and some other organization would become apparent. Alternatively, 5HT neurons could vary randomly, in which case single neuron profiles would poorly cluster. However, our data show that neither of these alternatives is correct; rather, a combination of lineage and anatomy is a meaningful way to define 5HT neuron subtypes in a manner that is consistent across population and single-neuron scales. In some cases, co-variation of gene expression observed across single neurons suggests finer-scale groupings within lineage and anatomy that is not apparent at the population scale, as demonstrated by R2 5HT neurons, which appear to resolve into two subtypes. Thus, lineage and anatomy classification serves as a first tier of organization, within which further subdivisions may be justified based on observed sets of co-varying genes. While present findings define overarching principles of 5HT neuron subtype organization, elaboration of finer scale principles will require additional work. For example, evidence supports regional variation in gene expression and functional properties within the DR, between the lateral wing region versus ventromedial portions, suggestive of additional distinct subtypes ((Spaethling et al., 2014); Niederkofler, Asher, and Dymecki, unpublished data).
Notably, the 5HT neuron subtypes defined here also appear to be functionally distinct. First, we found that the set of transcripts differentially expressed across 5HT neuron subtypes is enriched for gene categories critical for neuron function—ion channels, membrane receptors, neuropeptides, neurotransmitter synthetic enzymes, and cell adhesion molecules. Then, focusing on one particular 5HT neuron subtype as an example—R2 5HT neurons residing in the MR—we found that, in addition to molecular differences, these R2 5HT neurons displayed electrophysiological properties distinct from R1MR and R3 5HT neurons despite sharing anatomical location. Further, through synaptic silencing behavioral experiments, we found that regulation of sensorimotor gating (as assessed by PPI) maps to R2 5HT neurons, and not for example to R1DR and R1MR 5HT neurons. Deficits in sensorimotor gating as assessed by PPI are a hallmark of various neuropsychiatric disorders, such as schizophrenia (Braff et al., 2001); as such, the R2 5HT neuron subtype represents a previously undescribed functional node in sensory processing with potential disease relevance.
Further functional analysis of R2 5HT neurons showed that constituent subsets differentially respond to the ligands senktide and oxytocin, consistent with our single-neuron RNA-Seq data suggesting that the R2 5HT neuron subset comprises two distinct subtypes—one expressing the oxytocin receptor (defined as the R2-1 subtype) and another expressing the tachykinin 3 receptor (the R2-2 subtype). Thus, our multipronged functional probing of the R2 5HT neuron subtype supports the claim that these molecularly defined subsets of 5HT neurons are bona fide functionally distinct, and that many of the distinguishing gene expression differences are relevant to these functional differences. Our recent functional characterization of the R5 5HT neuron subtype further supports this claim, revealing subtype-specific electrophysiological properties, innervation patterns, and a unique role in the dynamic control of breathing (Brust et al., 2014)
In addition to identifying and characterizing 5HT neuron subtypes through unsupervised clustering and differential expression analyses, this resource permits the expression profile of any gene to be identified across the 5HT neuron system. Genes associated with human disorders are of particular clinical interest. Here we show that the autism-linked gene Met is expressed by a subset of 5HT neurons in the MR and caudal DR, in line with and extending prior findings (Wu and Levitt, 2013). Further, we show that loss of Met in 5HT neurons results in reduced social interest—a phenotype consistent with autism in humans—while having no effect on anxiety- or depression-like behaviors. Interestingly, in R2 5HT neurons, which show the highest overall levels of Met, Met expression is highly correlated with Oxtr, another gene associated with autism in multiple human populations (Campbell et al., 2011; Jacob et al., 2007; Lerer et al., 2008; Liu et al., 2010; Wu et al., 2005). These findings add to a growing body of evidence linking 5HT neuron-enriched genes with autism (Dougherty et al., 2013; Maloney et al., 2013), and further suggest that serotonergic regulation of social behaviors maps to specific subtypes of 5HT neurons, rather than to the system as a whole. Thus this resource is likely to facilitate discovery of not only function-specific 5HT neuron subpopulations, but also distinct disorder-related subpopulations—for example, differentially involved in anxiety disorders, sensory processing disorders, or disorders of respiratory control and physiological homeostasis, such as SIDS (Duncan et al., 2010; Kinney et al., 2009)—and therefore offers the potential to inform strategies for designing improved small molecule, antibody, and gene-based therapeutics for 5HT-associated disorders. While no single gene is likely to be expressed exclusively by a single 5HT neuron subtype and nowhere else in the body, therapies targeting a unique combination of gene products identified here may be sufficient for achieving 5HT subtype specificity of action and thus behavioral or physiological specificity either through multi-drug synergy or by designing therapeutic agents, such as allosteric modulators, that target subtype-specific heteromeric complexes. Additionally, the molecular profiles presented here may be used as transcriptional “goal states” for deriving specific 5HT neuron subtypes from stem cells or through transdifferentiation of suitable cell populations, as well as suggest combinations of subtype-specific transcription factors that can be used to drive the derivation process. Such 5HT neuron subtype-specific in vitro models will not only be valuable tools for high-throughput testing of candidate compounds, but offer the potential for deriving patient-specific 5HT neuron subtypes to study pathology and test patient-specific drug efficacy as well.
Experimental Procedures
Intersectional Genetic Fate Mapping
Triple transgenic mice were generated by breeding double transgenic Pet1::Flpe, RC::FrePe mice with En1-cre (En1Cki) (Kimmel et al., 2000), Rse2(Hoxa2)::cre (Awatramani et al., 2003), Egr2::cre (Voiculescu et al., 2000), or Tac1-IRES2-cre (Zheng, 2013). For further details, including primers for genotyping, see Supplemental Information.
Cell Sorting
Triple transgenic mice aged P25–P32 were euthanized in accordance with Harvard’s Institutional Animal Care and Use Committee Protocols. Fluorescently labeled cells were manually sorted as outlined in (Hempel et al., 2007). Sorted cells were deposited directly in cell lysis buffer (PicoPure RNA Isolation Kit; Applied Biosystems), heated for 30 minutes at 42° C, and stored at −80° C until continuing with the PicoPure RNA isolation protocol, which included an on-column DNAse digestion step, to eliminate possible contamination by genomic DNA.
RNA-Seq experiments
Total RNA harvested from pools (30–100 cells) or single cells was converted to cDNA and amplified using the Ovation RNA-Seq System v2 kit (Nugen), which utilizes a proprietary, non-PCR based, single primer isothermal amplification (SPIA) method that affords linear amplification of full length transcripts from low amounts of total RNA. Following SPIA amplification, ds cDNA was fragmented (~200 bp) using a Covaris S2 sonicator and ~50 ng of fragmented cDNA was input into the Ovation Ultralow DR Multiplex System (Nugen) to generate bar-coded libraries (note, only 12 PCR cycles were used in the final amplification step, as per manufacturer’s recommendation given input concentration). Quantification and quality control were assessed using TapeStation and Q-PCR. 50 base pair, single-end reads were generated on an Illumina HiSeq2000 or HiSeq2500 platform. Between 4 and 8 sample libraries were multiplexed in a single lane, resulting in an average of 29 million reads per library for pooled cell samples, and an average of 18 million reads for single-cell libraries. On average 91% of reads were successfully mapped for pooled cell libraries (75% uniquely), and 75% for single-cell libraries (62% uniquely). Sequencing data are available at GEO DataSets (accession number pending).
For details of RNA-Seq data analysis, electrophysiology, histology, smFISH, and behavioral methods, see Supplemental Information.
Intersectional Genetic Fate Mapping
Triple transgenic mice were generated by breeding double transgenic Pet1::Flpe, RC::FrePe mice with En1-cre (En1Cki) (Kimmel et al., 2000), Rse2(Hoxa2)::cre (Awatramani et al., 2003), Egr2::cre (Voiculescu et al., 2000), or Tac1-IRES2-cre (Zheng, 2013). For further details, including primers for genotyping, see Supplemental Information.
Cell Sorting
Triple transgenic mice aged P25–P32 were euthanized in accordance with Harvard’s Institutional Animal Care and Use Committee Protocols. Fluorescently labeled cells were manually sorted as outlined in (Hempel et al., 2007). Sorted cells were deposited directly in cell lysis buffer (PicoPure RNA Isolation Kit; Applied Biosystems), heated for 30 minutes at 42° C, and stored at −80° C until continuing with the PicoPure RNA isolation protocol, which included an on-column DNAse digestion step, to eliminate possible contamination by genomic DNA.
RNA-Seq experiments
Total RNA harvested from pools (30–100 cells) or single cells was converted to cDNA and amplified using the Ovation RNA-Seq System v2 kit (Nugen), which utilizes a proprietary, non-PCR based, single primer isothermal amplification (SPIA) method that affords linear amplification of full length transcripts from low amounts of total RNA. Following SPIA amplification, ds cDNA was fragmented (~200 bp) using a Covaris S2 sonicator and ~50 ng of fragmented cDNA was input into the Ovation Ultralow DR Multiplex System (Nugen) to generate bar-coded libraries (note, only 12 PCR cycles were used in the final amplification step, as per manufacturer’s recommendation given input concentration). Quantification and quality control were assessed using TapeStation and Q-PCR. 50 base pair, single-end reads were generated on an Illumina HiSeq2000 or HiSeq2500 platform. Between 4 and 8 sample libraries were multiplexed in a single lane, resulting in an average of 29 million reads per library for pooled cell samples, and an average of 18 million reads for single-cell libraries. On average 91% of reads were successfully mapped for pooled cell libraries (75% uniquely), and 75% for single-cell libraries (62% uniquely). Sequencing data are available at GEO DataSets (accession number pending).
For details of RNA-Seq data analysis, electrophysiology, histology, smFISH, and behavioral methods, see Supplemental Information.
Supplementary Material
1
2
3
4
1
2
3
4
Acknowledgments
The authors thank the Biopolymers Facility at HMS for assistance in next generation sequencing, especially K Waraska. Thank you to the HMS NBL and CHB IDDRC (P30 HD18655), including B Caldarone, NA Andrews, and G Gunner, for help with behavior assays. Thank you to MG Heiman for microscope use, T Huycke for smFISH support, P Tschopp for discussions on this manuscript, P Kharchenko for advice on single-cell RNA-Seq analysis, K Sugino for advice on RNA-seq workflows, R Kirchner for advice on RNA-seq analysis, A Karger for computing support, and JJ Mai for reagents and animal husbandry. Grants supporting this work include NIH P01 HD036379 (SD, BO), Harvard Stem Cell Institute Seed Grant (SD, BO), NARSAD Distinguished Investigator Grant from the Brain and Behavior Foundation (SD, BO), R21 MH083613 (SD, JCK), R21 DA023643 (SD, MF), American SIDS Institute research grant (BO), NIGMS T32GM007753 (MH), NIH R01 DA034022 (S.M.D.), and The Blavatnik Biomedical Accelerator at Harvard University Pilot Award (S.M.D., B.W.O., and B.D.R.).
Summary
Serotonergic (5HT) neurons modulate diverse behaviors and physiology and are implicated in distinct clinical disorders. Corresponding diversity in 5HT neuronal phenotypes is becoming apparent and is likely rooted in molecular differences, yet a comprehensive approach characterizing molecular variation across the 5HT system is lacking, as is concomitant linkage to cellular phenotypes. Here we combine intersectional fate mapping, neuron sorting, and genome-wide RNA-Seq to deconstruct the mouse 5HT system at multiple levels of granularity—from anatomy, to genetic sublineages, to single neurons. Our unbiased analyses reveal: principles underlying system organization, novel 5HT neuron subtypes, constellations of differentially expressed genes distinguishing subtypes, and predictions of subtype-specific functions. Using electrophysiology, subtype-specific neuron silencing, and conditional gene knockout, we show that these molecularly defined 5HT neuron subtypes are functionally distinct. Collectively, this resource classifies molecular diversity across the 5HT system and discovers new subtypes, markers, organizing principles, and subtype-specific functions with potential disease relevance.
Footnotes
Author Contributions: BO, MF, and SD designed experiments and wrote the manuscript; BO performed cell sorting, RNA-Seq, and RNA-Seq data analysis; MF performed immunohistology, smFISH, and behavior experiments and analysis; BR performed electrophysiology experiments and analysis. RB contributed to immunohistology and cell sorting; MH contributed to histology; DdB, JK, and MC contributed to behavior experiments.
Authors declare no conflicts of interest.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.