Combinatorial Expression of Grp and Neurod6 Defines Dopamine Neuron Populations with Distinct Projection Patterns and Disease Vulnerability

Abstract Midbrain dopamine neurons project to numerous targets throughout the brain to modulate various behaviors and brain states. Within this small population of neurons exists significant heterogeneity based on physiology, circuitry, and disease susceptibility. Recent studies have shown that dopamine neurons can be subdivided based on gene expression; however, the extent to which genetic markers represent functionally relevant dopaminergic subpopulations has not been fully explored. Here we performed single-cell RNA-sequencing of mouse dopamine neurons and validated studies showing that Neurod6 and Grp are selective markers for dopaminergic subpopulations. Using a combination of multiplex fluorescent in situ hybridization, retrograde labeling, and electrophysiology in mice of both sexes, we defined the anatomy, projection targets, physiological properties, and disease vulnerability of dopamine neurons based on Grp and/or Neurod6 expression. We found that the combinatorial expression of Grp and Neurod6 defines dopaminergic subpopulations with unique features. Grp+/Neurod6 + dopamine neurons reside in the ventromedial VTA, send projections to the medial shell of the nucleus accumbens, and have noncanonical physiological properties. Grp+/Neurod6- dopamine neurons are found in the VTA as well as in the ventromedial portion of the SNc, where they project selectively to the dorsomedial striatum. Grp-/Neurod6 + dopamine neurons represent a smaller VTA subpopulation, which is preferentially spared in a 6-OHDA model of Parkinson’s disease. Together, our work provides detailed characterization of Neurod6 and Grp expression in the midbrain and generates new insights into how these markers define functionally relevant dopaminergic subpopulations.


Introduction
Midbrain dopamine (DA) neurons of the substantia nigra pars compacta (SNc) and ventral tegmental area (VTA) make widespread projections throughout the brain and modulate a host of behaviors from motor function to reward learning to cognition (Björklund and Dunnett, 2007). Although they represent only ϳ0.03% of neurons in the mouse brain, DA neurons are heterogeneous, as they vary significantly in their circuitry (Lammel et al., 2011;Watabe-Uchida et al., 2012;Beier et al., 2015;Menegas et al., 2015), physiology (Lammel et al., 2008;Lerner et al., 2015), gene expression (Poulin et al., 2014;La Manno et al., 2016), and response to disease (Damier et al., 1999;Di Salvio et al., 2010a;Blesa and Przedborski, 2014;Brichta and Greengard, 2014). VTA DA neurons are particularly diverse, comprising multiple subcircuits that project to different brain regions and have distinct electrophysiological properties according to their projection target (Lammel et al., 2008;Margolis et al., 2008;Morales and Margolis, 2017). Depending on their connectivity, VTA neurons can also mediate opposing behaviors, such as reward and aversion (Lammel et al., 2012), necessitating tools that can parse this functional heterogeneity to allow selective manipulation of specific VTA subpopulations (Poulin et al., 2016;Vogt Weisenhorn et al., 2016). Building on prior studies that identified genetic differences between SNc and VTA neurons (Grimm et al., 2004;Chung et al., 2005;Greene et al., 2005;Di Salvio et al., 2010b), recent single-cell gene profiling studies have uncovered further genetic heterogeneity in the DA system, including several subtypes of VTA neurons (Poulin et al., 2014;La Manno et al., 2016). However, it is currently unknown how these genetically defined classes of DA neurons map onto subtypes defined by their circuitry and physiology.
A notable feature of dopaminergic subpopulations is their differential vulnerability to disease. For example, in the neurodegenerative disorder Parkinson's disease (PD), SNc DA neurons degenerate earlier and to a greater degree than VTA DA neurons (Damier et al., 1999;Alberico et al., 2015). The reason for this selective vulnerability is not well understood, although current hypotheses point to differences in the expression of ion channels and metabolic proteins between SNc and VTA neurons (Chung et al., 2005;Greene et al., 2005;Liss et al., 2005;Chan et al., 2007;Brichta and Greengard, 2014;Liu et al., 2014). Despite the relative sparing of the VTA compared to the SNc, 40%-77% of VTA DA neurons still degenerate (Alberico et al., 2015) and the molecular features that define susceptible versus spared VTA neurons are unknown.
Here our goal was to define dopaminergic subpopulations based on gene expression and determine how these populations map onto DA subtypes defined by physiology and circuitry. To do this, we analyzed DA neuron populations marked by two genes, Grp and Neurod6, that we identified by single-cell RNA-sequencing (RNA-seq), and which have previously been reported to mark subpopulations of VTA DA neurons (Chung et al., 2005;Greene et al., 2005;Poulin et al., 2014;La Manno et al., 2016;Viereckel et al., 2016;Khan et al., 2017). With a combination of anatomy, retrograde tracing, and physiology, we show that these genes define overlapping yet distinct DA neuron populations. We further demonstrate that the combinatorial expression of these two genes influences susceptibility to degeneration in a 6-OHDA mouse model of PD. Together, our findings further our understanding of dopaminergic cell type diversity and validate genetic approaches for defining functional cell types in the brain.

Mice
Animal procedures were conducted in accordance with protocols approved by the University of California, Berkeley Institutional Animal Care and Use Committee (IACUC) and Office of Laboratory Animal Care (OLAC).
For single-cell RNA-seq experiments, DAT IRES Cre mice (Bäckman et al., 2006; Jackson Laboratory strain #006660, RGD_12905031) were crossed and maintained with the Ai9 tdTomato Cre-reporter line (Madisen et al., 2010; Jackson Laboratory strain #007909). For physiology experiments, NEX-Cre mice were obtained from Dr. Klaus-Armin Nave (Goebbels et al., 2006) and crossed with the Ai9 mouse line. C57BL/6J mice were used for retrograde bead injections. The ages, sexes, and numbers of mice used are indicated for each experiment in the results and figure legends.
This single-cell suspension was sorted on a BD Influx cell sorter in the Flow Cytometry Facility at UC Berkeley. Cells were gated for PI-/tdTomato ϩ and sorted into a PCR tube. Based on the number of neurons sorted and cell viability count, neurons were brought to ϳ200,000 neurons/ml. Neurons were then put into a large Fluidigm C1 chip, and each of the 96 wells was visually inspected to verify cell presence, cell health, and tdTomato expression. Wells containing cell doublets were excluded from further processing.
Cells went through single-cell mRNA extraction using the Fluidigm C1 system in the Functional Genomics Laboratory at UC Berkeley. Single-cell cDNA was removed and measured via Qubit. Any cell that gave Ͻ0.3 g/ml of cDNA was removed due to likely low quality. cDNA from single cells that passed the initial quality check was diluted to 0.3 g/ml. 379 single-neuron cDNA extracts were library prepped using the Nextera DNA library prep protocol (Illumina #FC-121-1012). The cDNA was then sequenced on a HiSeq 2500 in the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley.

Single-cell RNA-seq preprocessing
Reads were aligned to the GRCm38.3 (mm10, patch release 3) mouse genome assembly with Tophat2 (v. 2.1.1; Kim et al., 2013), and low-quality reads were removed with Trimmomatic (v. 0.3.2; Bolger et al., 2014). Gene expression was quantified using featureCounts (v. 1.5.0-p3; Liao et al., 2014) and RefSeq transcript annota-tion. Reads that aligned to more than one gene as well as chimeric fragments were excluded. We used a quality control (QC) pipeline that computes an extensive set of quality metrics, relying in part on FastQC (v. 0.3.2) and the Picard tools (v. 2.5.0 with samtools 1.3.1) as done previously (Fletcher et al., 2017). We used the Bioconductor package scone (https://bioconductor.org/packages/scone; v. 0.99.6) to perform data-adaptive cell and gene filtering. This yielded the following exclusion criteria: any cell with fewer than 500,000 aligned reads or a percentage of aligned reads below 85%. In addition, we filtered out cells with large dropout rates, as defined by the "false-negative curves" of scone. This procedure resulted in a total of 232 retained cells (out of 379). Finally, we retained only those genes having at least 10 reads in at least 10 cells (10,983 genes).

Single-cell RNA-seq statistical analysis
We performed and assessed several normalization schemes using scone  and selected fullquantile normalization (Bolstad et al., 2003;Bullard et al., 2010). We then applied principal component analysis on the normalized data and retained the first 50 principal components, which explained 41% of the variance. Cluster analysis was performed on the first 50 principal components using the RSEC method (Risso et al., 2018b) implemented in the Bioconductor package clusterExperiment (https://bioconductor.org/packages/clusterExperiment; v. 1.0.0), as previously described (Fletcher et al., 2017). We used RSEC with the following specific parameters: alphas ϭ 0.3, minSizes ϭ 5, combineProportions ϭ 0.5; all other parameters were left at their default values. RSEC found 9 stable clusters. We used limma (v. 3.30.13) with voom correction weights (Law et al., 2014), as implemented in the clusterExperiment function getBestFeature, to find marker genes for each cluster. To visualize the clustering results, we applied ZINB-WaVE (Risso et al., 2018a; with K ϭ 10) for dimensionality reduction, followed by t-distributed stochastic neighbor embedding (t-SNE; perplexity parameter set to 20; van der Maaten and Hinton, 2008).

Accession number and code accessibility
The accession number for the RNA-seq data reported in this paper is GEO: GSE115070. The computer code for the single-cell RNA-seq analysis is available at https:// github.com/drisso/striatum. This code was run on an Apple Mac computer with macOS Sierra 10.12.4 operating system.

Fluorescent in situ hybridization
To visualize mRNA we used the RNAScope fluorescent in situ hybridization (FISH) method (ACDBio). Fresh-frozen tissue was processed as per RNAScope instructions. Briefly, mouse brain tissue from male and female mice aged P21-P120 was fresh-frozen in OCT on dry ice and stored at least overnight at -80˚C. Tissue was then cut on a cryostat (Leica Microm HM550) into 12 m sections and mounted onto slides. Slides were fixed in 4% PFA in 1ϫ PBS for 15 min. Slides were dehydrated using 5 min incubations in 50%, 70%, and twice with 100% ethanol.

Microscopy and image analysis
Two confocal microscopes were used to take Z-stack images of FISH-labeled or immunostained sections: an Olympus FV1000 with a 20ϫ Nikon objective and a Zeiss LSM 710 AxioObserver with Zeiss 10ϫ, 20ϫ, and 63ϫ objectives housed in the Molecular Imaging Center at UC Berkeley. Images were analyzed using the Fiji image analysis toolbox. Cells were considered positive for Neurod6, Grp, or Lypd1 if they contained three fluorescent puncta within the boundary created by a cell marker: Slc6a3 (DAT), Th, or DAPI.
To allow for sufficient DA neuron labeling, mice were sacrificed at various time points following injection: 7 d for the DMS and DLS, 21 d for the prefrontal cortex, and 14 d for the NAc (medial shell, core, and lateral shell), lateral septum and amygdala. Brains were harvested and cut to separate the injection site and midbrain. The midbrain was frozen for cryostat sectioning as described above for FISH. The brain region containing the injection site was incubated in 4% PFA overnight at 4˚C, then cryoprotected in 30% sucrose in 1ϫ PB until the tissue sank. Injection site tissue was sectioned on a freezing microtome, mounted with Vectashield hardset mounting media with DAPI (Vector laboratories #H-1500), and analyzed for bead expression at the injection site. Brains with correctly targeted injection sites and minimal off-target bead expression were chosen for analysis.
Mice were monitored daily following the injection to ensure recovery. Kitten Milk Replacement (Santa Cruz #sc-362120) was fed to mice daily for up to 2 wk following the injection to aid recovery and meloxicam was injected subcutaneously to alleviate pain if necessary. Motor function was assessed using the cylinder test each week following the injection (see below). 4 wk following 6-OHDA injection, mice were quickly anesthetized using isoflurane and decapitated, and their brains were fresh-frozen as described above for FISH.

Cylinder test
To test the severity of Parkinsonian-like symptoms following unilateral 6-OHDA injection, we used the cylinder test to score limb use asymmetry. Mice were habituated to the behavior room for a minimum of 30 min during their dark cycle under red light illumination. Mice were placed into a clear plastic cylinder 12 cm in diameter and 20 cm in height. The cylinder was placed next to two mirrors to visualize paw use. The mouse was both video recorded and observed by the experimenter while it was allowed to move around freely in the cylinder for 10 min. Full 360°i psiversive and contraversive rotations (relative to the 6-OHDA injection side) were counted. Forelimb asymmetry was measured by counting the number of times the ipsilateral paw, contralateral paw, or both paws were used for support when the mouse reared against the wall of the cylinder. The percentage of ipsilateral or contralateral paw use was calculated based on total rears (e.g., ipsilateral paw touches/ipsilateral ϩ contralateral ϩ both paw touches). A greater number of ipsiversive turns and ipsilateral paw use indicates a successful 6-OHDA injection and unilateral Parkinsonian-like symptoms.

Experimental design and statistical analysis
A one-way ANOVA was used to compare the means of three or more groups. Post hoc pairwise comparisons were made using either Bonferroni's or Tukey's multiple comparisons tests. Unpaired or paired two-tailed t tests were used to compare the means of two groups. A paired, one-way ANOVA with Dunnett's multiple comparisons post hoc test was used to compare DA neuron subpopulations to the entire VTA DA population for the 6-OHDA experiments. Data are reported as the mean Ϯ SEM. Superscript letters listed with p-values correspond to the statistical tests shown in Table 1. p-values were corrected for multiple comparisons.

Single-cell RNA-seq of midbrain dopamine neurons
To define subclasses of DA neurons in an unbiased way, we performed single-cell RNA sequencing of DA neurons from P26 -P34 male and female mice in which dopamine transporter (DAT)-expressing neurons were labeled with a tdTomato reporter (DAT IRES Cre;Ai9, Fig.  1-1A). A bioinformatic and statistical workflow revealed nine clusters of DA neurons based on differential gene expression ( Fig. 1-1B-E and Table 1-1). We identified markers for each cluster and validated eight of the nine clusters. Two of the clusters corresponded to DA subpopulations in the hypothalamus, four defined subclasses of VTA neurons, and two corresponded to the SNc ( Fig.  1-1F). These populations are consistent with recent single-cell DA neuron profiling studies (Poulin et al., 2014;La Manno et al., 2016;Romanov et al., 2017).

Grp and Neurod6 define anatomically overlapping but distinct midbrain subpopulations
Consistent with prior studies (Poulin et al., 2014;La Manno et al., 2016), we identified a cell cluster that showed relatively high and selective expression of Grp, Neurod6, and Gpr83 (cluster #8, Fig. 1-1C). Given recent interest in Neurod6 and Grp as markers that label ventromedial VTA DA neurons (La Manno et al., 2016;Viereckel et al., 2016;Khan et al., 2017), we chose to quantitatively analyze their expression patterns in the midbrain and examine the extent of their overlap. Using multiplex FISH, we found that in adult mice, Grp-expressing neurons represented 29.9% Ϯ 0.5% of the total midbrain DA neuron population ( Neurod6 expression defined a more restricted DA population, accounting for 26.8% Ϯ 0.5% of VTA DA neurons (1172/4377 neurons from 8 mice), with the highest density of Neurod6 ϩ DA neurons in the ventromedial IF and PN/ PIF regions (p Ͻ 0.0001, b one-way ANOVA with Tukey's post hoc test, Fig. 1E). In contrast to Grp, Neurod6 mRNA was not expressed in the SNc (Fig. 1A-F). 93.1% Ϯ 0.6% of Neurod6 ϩ VTA neurons were dopaminergic (948/1016 cells from 4 mice). Consistent with our RNA-seq data, we found that most, but not all, Neurod6 ϩ DA neurons coexpressed Grp (77.5% Ϯ 0.8%, 909/1172 cells from 8 mice). Broken down by subregion, the majority of Neurod6 ϩ DA neurons in the IF and PN/PIF coexpressed Grp; however, a third of Neurod6 ϩ DA neurons in the PBP did not have   (Fig. 1F). Grp-expressing but Neurod6-lacking DA neurons were found throughout the VTA, and nearly all of the Grp ϩ DA neurons in the SNc were lacking Neurod6 (Fig. 1F). Together, these data show that Grp is expressed in a subpopulation of midbrain DA neurons spanning the VTA and ventromedial portion of the SNc. Neurod6 is expressed in a more restricted midbrain population, defining a subgroup of the Grp ϩ DA population, which is located exclusively in the VTA. In addition, there is a small subset of VTA DA neurons that express Neurod6 but not Grp.
To determine if these cell populations were present throughout development, we assessed Grp and Neurod6 expression at multiple ages. We found that while there were subtle increases in the percentages of Grp ϩ VTA DA neurons from 8 to 16 wk (p ϭ 0.036, c one-way ANOVA with Tukey's post hoc test, Fig. 1G) and Neurod6 ϩ VTA DA neurons from 2 to 16 wk (p ϭ 0.024, d one-way ANOVA with Tukey's post hoc test, Fig. 1H), the proportions of DA neurons expressing these markers were largely stable from 2 wk postnatal through adulthood (P14 -P112). We also investigated whether these cell populations were similar between sexes. We found that the Grp-and Neurod6-expressing DA subpopulations were present in male and female mice in similar proportions throughout VTA and SNc subregions (see statistics worksheet for unpaired t test p values, e-m Fig. 1I,J). Furthermore, in the RNA-seq data, we found no significant differences in expression levels of Grp or Neurod6 between males and females (log2-fold-change of 1.70 for Neurod6, p ϭ 0.51 n and 0.94 for Grp, p ϭ 0.68 o in males versus females, likelihood ratio test). These data demonstrate that the Grp ϩ and Neurod6 ϩ DA neuron populations are present in both sexes and stable over time.

Grp and Neurod6-expressing neurons project to the nucleus accumbens medial shell
VTA DA neurons comprise multiple subcircuits, which send projections to different brain regions with distinct functional consequences (Roeper, 2013). To determine if the Grp and Neurod6 DA populations have specific projection targets, we combined retrograde labeling from eight primary DA neuron projection sites with Grp and Neurod6 FISH ( Fig. 2A and Fig. 2-1A). We found that of the DA neurons projecting to the medial shell of the nucleus accumbens (NAc), 75.0% Ϯ 1.1% were Grppositive and 70.4% Ϯ 2.5% were Neurod6-positive, indicating that these markers were expressed in the majority of medial shell-projecting DA neurons (Fig. 2B,C,F,H). Grp ϩ and Neurod6 ϩ DA neurons also projected to the NAc core and lateral shell but represented a smaller fraction of the neurons projecting to these regions compared to the medial shell (Fig. 2F,H). When quantified as a percentage of the total marker-positive bead-labeled DA neurons across all injection sites, the NAc medial shell was the primary target region for both Grp ϩ and Neurod6 ϩ DA neurons (Fig. 2G,I). Together, these findings indicate a strong mesoaccumbens projection from Neurod6 ϩ and Grp ϩ VTA DA neurons.
We found that Neurod6 mRNA was largely absent from DA neurons projecting to regions outside of the NAc, indicating a selective output of the Neurod6 ϩ DA neuron population (Fig. 2H,I). By contrast, Grp ϩ DA neurons represented a substantial percentage (47.3% Ϯ 2.4%) of DA neurons projecting to the dorsomedial striatum (DMS; Fig.  2D-G). 87.7% Ϯ 2.5% of the DMS-projecting Grppositive neurons were located in the ventromedial portion of the SNc (Fig. 2-1B), consistent with prior reports mapping dopaminergic projections to the DMS (Lerner et al., 2015). These neurons lacked expression of Neurod6, as essentially no DMS-projecting neurons were Neurod6positive (Fig. 2H,I). Taken together, these results indicate that there are at least two populations of DA neurons that express Grp: those located in the VTA that project primarily to the medial NAc and those in the ventromedial SNc, which project selectively to the DMS (Fig. 2-1B,C).
A recent study showed that Neurod6-expressing DA neurons, labeled by a Cre reporter in NEX-Cre mice [Neu-  rod6 was previously referred to as NEX (Goebbels et al., 2006)], project to the lateral septum (Khan et al., 2017). We found relatively few lateral septum-projecting neurons in the midbrain, of which 37.4% Ϯ 3.1% were nondopaminergic (92/252 cells from 4 mice). Of the DA neurons projecting to the lateral septum, 17.0% Ϯ 2.3% were Grp ϩ and 24.0% Ϯ 3.0% were Neurod6 ϩ (Fig. 2F,H). Quantified as a percentage of total bead-labeled neurons, Ͻ2% of Grp ϩ or Neurod6 ϩ DA neurons projected to the lateral septum (Fig. 2G,I). These results indicate that compared to other brain regions, the lateral septum is not a major target for Grp-or Neurod6-expressing DA neurons.

Neurod6 ؉ DA neurons have unique physiological properties
DA neurons projecting to different target areas possess distinct electrophysiological profiles (Lammel et al., 2008;Margolis et al., 2008;Lerner et al., 2015). To investigate whether Neurod6 expression defines a physiologically distinct subclass of DA neurons, we used NEX-Cre knock-in mice (Goebbels et al., 2006). We injected virus expressing a Cre-dependent tdTomato reporter (AAV-Flex-tdTomato) into the midbrain and found NEX-Cre ϩ DA neurons along the ventromedial portion of the VTA (Fig.  3A), consistent with the expression pattern of Neurod6 mRNA. In agreement with our tracing data, we found that tdTomato-labeled NEX-Cre ϩ neurons projected to the NAc medial shell and core (Fig. 3B). To visualize NEX-Cre ϩ neurons for physiology, we crossed NEX-Cre mice to the Ai9 tdTomato Cre-reporter mouse line (Madisen et al., 2010). We performed FISH for Neurod6 mRNA in NEX-Cre;Ai9 mice ( Fig. 3C) and found that 93.6% Ϯ 1.2% of tdTomato-positive DA neurons in the VTA coexpressed Neurod6 (331/353 cells from 4 mice), making this a suitable model to use for targeted electrophysiology recordings. Consistent with the 77.5% of Neurod6 ϩ neurons that coexpressed Grp, 73.7% Ϯ 2.9% of the NEX-Cre; tdTomato-positive VTA DA neurons expressed Grp (175/ 232 cells from 4 mice). We did observe that not all Neu-rod6 ϩ VTA DA neurons expressed tdTomato in the NEX-Cre;Ai9 mice (30.4% Ϯ 1.7% of Neurod6 ϩ DA neurons were tdTomato-positive, 361/1092 cells from 4 mice). In addition, a third of the tdTomato-positive midbrain neurons were non-dopaminergic (36.9% Ϯ 3.4% Th negative, 193/546 cells from 4 mice). This discrepancy may be due to transient Cre expression in non-DA neurons during development. These data suggest that NEX-Cre mice may not be appropriate for studies requiring selective access to the entire Neurod6 ϩ VTA DA subpopulation, but can be used to target Neurod6 ϩ cells for whole-cell recordings in which dopaminergic identity can be confirmed post hoc.
To determine if Neurod6-expressing neurons represent a functionally distinct cell class, we recorded from tdTomato-labeled NEX-Cre ϩ neurons in the VTA and analyzed their intrinsic membrane properties, action potential waveform, and response to hyperpolarizing current injection (see Fig. 3-1 for a summary of the physiology data including sample sizes). We confirmed that the NEX-Cre ϩ neurons analyzed were dopaminergic by filling patched neurons with neurobiotin and costaining for TH (Fig. 3D,E). We found that NEX-Cre ϩ DA neurons had a distinct electrophysiological signature compared to "classic" SNc DA neurons. Specifically, NEX-Cre ϩ DA neurons exhibited a more depolarized membrane potential (Vm, p ϭ 0.0075, p unpaired t test), had higher membrane resistance (Rm, p Ͻ 0.0001, q unpaired t test), and reduced membrane capacitance (Cm, p Ͻ 0.0001, r unpaired t test) compared to SNc DA neurons (Fig. 3F-H). The action potential height of NEX-Cre ϩ DA neurons was also significantly shorter (p Ͻ 0.0001, s unpaired t test), and they had a less pronounced afterhyperpolarization (AHP, p ϭ 0.0002, t unpaired t test; Fig. 3I-K). Due to their high membrane resistance, NEX-Cre ϩ DA neurons required less negative current to hyperpolarize to -100 mV compared to SNc DA neurons (-25 to -50 pA for NEX-Cre ϩ neurons versus -150 pA for SNc neurons). NEX-Cre ϩ DA neurons had a significantly smaller sag component, which is indicative of smaller I h (p Ͻ 0.0001, u unpaired t test, Fig.  3L-N) and less rebound depolarization (p Ͻ 0.0001, v unpaired t test, Fig. 3L,M,O). The noncanonical electrophysiological characteristics of NEX-Cre ϩ neurons are consistent with those reported for NAc medial shellprojecting DA neurons (Lammel et al., 2008), suggesting that Neurod6 is a useful marker for this VTA subpopulation.

Neurod6-lacking VTA DA neurons show increased susceptibility to degeneration in a 6-OHDA mouse model
In addition to their anatomic location, projection targets, and physiology, vulnerability to degeneration is a key feature of DA neurons that differs by subtype. Previous in vitro studies have implicated both Neurod6 and Grp as being potentially neuroprotective (Chung et al., 2005;Uittenbogaard and Chiaramello, 2005;Uittenbogaard et al., 2010;Baxter et al., 2012). In mice, expression of Neurod6 and the related transcription factor Neurod1 are important for survival of DA neurons during development (Khan et al., 2017). Grp-expressing cells have been observed in postmortem tissue from PD patients, suggesting that Grp expression may be related to cell survival (Viereckel et al., 2016). We therefore investigated the sensitivity of Neu-rod6 and/or Grp-expressing VTA DA neurons to degeneration in a mouse model of PD. To do this, we injected the dopaminergic toxin 6-hydroxydopamine (6-OHDA) unilaterally into the medial forebrain bundle of adult (P120) continued Bars represent mean Ϯ SEM, each dot represents one mouse. I, Summary of the projection targets of Neurod6 ϩ DA neurons expressed as a percentage of total bead ϩ /Neurod6 ϩ /DAT ϩ neurons. For panels F-I: NAc MSh n ϭ 850 cells from 2 male and 2 female mice, NAc Core n ϭ 858 cells from 2 males and 2 females, NAc LSh n ϭ 1215 cells from 3 males and 1-2 females, DMS n ϭ 637 cells from 3 males and 2 females, DLS n ϭ 1038 cells from 2-4 males and 2-3 females, mPFC n ϭ 211 cells from 3-4 males and 1 female, BLA n ϭ 297 cells from 3 males, and LS n ϭ 345 cells from 2 males and 4 females. See also Fig. 2-1.  (Thiele et al., 2012). 6-OHDA resulted in a progressive, unilateral loss of DA neurons with a 96.0% Ϯ 0.6% reduction in SNc DA neurons and a 69.3 Ϯ 1.1% loss of VTA DA neurons after 4 wk (Fig. 4A). We confirmed DA axon denervation throughout the striatum in the 6-OHDA-injected hemisphere, with notable sparing of DA projections to the NAc medial shell and core (Fig. 4B). DA neuron loss led to unilateral motor impairment as measured by the cylinder test (p ϭ 0.0002, w one-way ANOVA with Tukey's post hoc test: saline paw vs. 6-OHDA paw p ϭ 0.0163, saline paw vs. both paws p ϭ 0.0157, 6-OHDA paw vs. both paws p ϭ 0.0002, n ϭ 4 mice).
To determine how Grp-and Neurod6-expressing DA neurons in the VTA responded to 6-OHDA ( Fig. 4C-F), we compared the reduction of Grp ϩ and/or Neurod6 ϩ VTA DA neurons between the saline and 6-OHDA hemisphere to all VTA DA neurons defined by Th expression (Fig. 4G). For this analysis, we included only VTA DA neurons, as essentially all SNc neurons (including Grp ϩ SNc neurons) degenerated in the 6-OHDA-injected hemisphere (Fig.  4A,C,D). We found that 4 wk following 6-OHDA injection, 26.1% Ϯ 2.7% of Grp ϩ /Neurod6 ϩ DA neurons survived, which was similar to the percentage of total Th ϩ VTA neurons surviving (30.7% Ϯ 1.1%; p Ͻ 0.0001, x paired one-way ANOVA; p ϭ 0.3169, Tukey's post hoc test; Fig.  4G). Notably, VTA DA neurons that expressed only Neu-rod6 and not Grp (Grp-/Neurod6 ϩ ) were significantly spared compared to the rest of VTA DA neurons, with 82.7% Ϯ 4.2% of neurons surviving (p ϭ 0.0022, Tukey's post hoc test, Fig. 4G). By contrast, VTA DA neurons that expressed Grp but not Neurod6 (Grp ϩ /Neurod6 -) showed slightly increased vulnerability compared to all VTA DA neurons, with 21.7% Ϯ 2.0% surviving (p ϭ 0.0246, Tukey's post hoc test, Fig. 4G). Therefore, the DA neuron subpopulations marked by expression of Grp and/or Neu-rod6 have different responses to 6-OHDA, and VTA neurons lacking Neurod6 are more susceptible to degeneration.
An increase in the proportion of Neurod6 ϩ DA neurons in the VTA following 6-OHDA could be due to selective sparing of the Neurod6 cell population or from Neurod6 expression turning on in surviving neurons of other VTA subpopulations. To attempt to distinguish these possibilities, we performed FISH for another VTA DA neuron subtype marker Lypd1 (La Manno et al., 2016).
Under normal conditions, Lypd1 was expressed almost exclusively in the PBP subregion of the VTA as well as throughout the SNc (Fig. 4D,F). Neurod6 was generally not coexpressed in Lypd1 ϩ DA neurons (7.0% Ϯ 1.6% of Lypd1 ϩ DA neurons coexpressed Neurod6, Fig. 4D,F,H), indicating that these markers define distinct cell popula-tions. In response to 6-OHDA, we found that 20.4% Ϯ 1.3% of Lypd1 ϩ /Neurod6-VTA DA neurons survived, which was significantly lower than the total VTA DA population (p ϭ 0.0403, Tukey's post hoc test, Fig. 4G). The surviving Lypd1 ϩ VTA DA neurons did not turn on expression of Neurod6, as the proportion of Neurod6 ϩ /Lypd1 ϩ VTA DA neurons remained low and was the same between the saline-injected and 6-OHDA-injected hemispheres (p ϭ 0.8091, y paired t test, Fig. 4H). These results suggest that the VTA DA subpopulation defined by Neu-rod6 expression, which lacks Grp or Lypd1, is preferentially spared in response to a neurotoxic challenge.

Discussion
Midbrain DA neurons are small in number but vast in their behavioral influence (Björklund and Dunnett, 2007;Roeper, 2013). As such, dopaminergic dysfunction is associated with numerous psychiatric and neurologic disorders ranging from drug addiction to PD (Damier et al., 1999;Alberico et al., 2015;Nutt et al., 2015). Recent studies have revealed that the dopaminergic system is heterogeneous at multiple levels from gene expression to circuitry to physiology to behavior (Roeper, 2013;Morales and Margolis, 2017). To tackle this heterogeneity, genetic markers that define functional DA neuron subtypes are needed. This would enable the generation of tools that allow selective manipulation of dopaminergic subpopulations. Here, we investigated two markers that we and others have identified as labeling dopaminergic subpopulations, Grp and Neurod6. We show that the combinatorial expression of these genes defines the anatomic location, projection target, physiology, and disease susceptibility of DA neurons.
We found that Grp, which encodes the neuropeptide gastrin-releasing peptide (Roesler and Schwartsmann, 2012), was expressed in a third of midbrain DA neurons, of which more than half coexpressed Neurod6. These Grp ϩ /Neurod6 ϩ neurons resided in the VTA and projected to the medial portions of the NAc. This is consistent with prior reports showing that Grp is expressed in a subpopulation of VTA DA neurons that shows overlap with Neurod6-expressing neurons (Chung et al., 2005;Greene et al., 2005;Poulin et al., 2014;La Manno et al., 2016;Viereckel et al., 2016). The fact that these neurons project to the NAc corroborates a projection-specific translational profiling study reporting that ribosome-bound Grp mRNA was enriched in the population of VTA DA neurons projecting to the NAc (Ekstrand et al., 2014). Notably, we also identified a previously undiscovered population of Grp ϩ DA neurons that lack Neurod6, which were located in the continued neurons positive for TH. F-I, Graphs display the mean Ϯ SEM membrane potential (F, Vm, ‫,ءء‬ p ϭ 0.0075), membrane resistance (G, Rm, ‫,ءءء‬ p Ͻ 0.0001), capacitance (H, Cm, ‫,ءءء‬ p Ͻ 0.0001), and afterhyperpolarization (I, AHP, ‫,ءءء‬ p ϭ 0.0002) for NEX-Cre ϩ and SNc neurons. J, K, Representative action potential traces and phase plots (mVms Ϫ1 /mV) from a NEX-Cre ϩ (J) and SNc (K) neuron. L, M, Representative responses to a 500 ms negative current step from a NEX-Cre ϩ (-25 pA, L) and SNc (-150 pA, M) neuron. N, Quantification of the mean Ϯ SEM sag component in NEX-Cre ϩ and SNc neurons ‫,ءءء(‬ p Ͻ 0.0001). O, Quantification of the mean Ϯ SEM rebound depolarization from neurons hyperpolarized to -100 Ϯ 7 mV ‫,ءءء(‬ p Ͻ 0.0001). An unpaired t test was used for all comparisons. For all panels, dots represent values from individual cells, NEX-Cre ϩ n ϭ 7 male and 7 female mice, NEX-Cre-/SNc n ϭ 4 male mice and 1 female mouse. See Fig. 3-1 for a complete summary of the electrophysiology results. , Neurod6 ϩ (red), and Lypd1 ϩ /Neurod6 ϩ (yellow) DA neurons in the midbrain following saline and 6-OHDA injection. E, F, Confocal images of FISH for the indicated markers. E, Images show Grp ϩ (blue circle), Neurod6 ϩ (red circles), and Grp ϩ /Neurod6 ϩ (purple circles) DA neurons in the VTA following saline and 6-OHDA injection. F, Images show Lypd1 ϩ (green circles), Neurod6 ϩ (red circles), and Lypd1 ϩ /Neurod6 ϩ (yellow circle) DA neurons in the VTA following saline and 6-OHDA injection. G, Box-and-whisker plots (min to max) show the percentage of VTA DAT neurons expressing each set of markers in the 6-OHDA-injected hemisphere compared to the saline-injected hemisphere. Dotted line represents the percentage of all VTA DA neurons (Th ϩ ) surviving in the 6-OHDA hemisphere. Dots represent data from individual mice, n ϭ 4 female mice. A paired one-way ANOVA with Dunnett's multiple comparisons test was used to compare each subpopulation to all Th ϩ VTA DA neurons (n ϭ 4376 saline and 1352 6-OHDA cells): Grp -/Neurod6 ϩ (n ϭ 184 saline and 151 6-OHDA cells), ‫ءء‬ p ϭ 0.0022; Grp ϩ /Neurod6 ϩ (n ϭ 763 saline and 207 6-OHDA cells), p ϭ 0.3169; Grp ϩ /Neurod6 -(n ϭ 749 saline and 165 6-OHDA cells), ‫ء‬ p ϭ 0.0246; Lypd1 ϩ /Neurod6 -(n ϭ 1056 saline and 171 6-OHDA cells), ‫ء‬ p ϭ 0.0313; Nd6 ϭ Neurod6. H, Bar graphs display the percentage of Lypd1 ϩ VTA DA neurons that coexpress Neurod6 mRNA in the saline-injected (n ϭ 102/3750 cells) and 6-OHDA-injected (n ϭ 16/904 cells) hemispheres. Bars represent mean Ϯ SEM, dots represent the values from individual mice, n ϭ 4 female mice; paired t test, n.s., not significant, p ϭ 0.8091. ventromedial portion of the SNc. These neurons sent projections to the DMS with very little innervation of the DLS. The anatomic location of these cells in the ventral SNc is consistent with DA neurons that project to the DMS, which have unique physiological and behavioral properties compared to DLS-projecting DA neurons located in the dorsal SNc (Lerner et al., 2015).
Neurod6 expression defined a smaller population of DA neurons that were located in the ventromedial VTA, projected selectively to the NAc, and exhibited noncanonical physiologic properties. While Neurod6 has previously been identified as a VTA marker (Chung et al., 2005;La Manno et al., 2016;Viereckel et al., 2016;Khan et al., 2017), our study is the first to systematically map the projection sites of these cells, revealing a strong medial shell NAc projection with very little output to other DA target regions. This medial shell NAc projection is consistent with the physiology of Neurod6 ϩ DA neurons, which showed unique characteristics similar to those previously reported for medial NAc-projecting DA neurons defined by retrograde labeling (Lammel et al., 2008). Mesoaccumbens projections are important for mediating reward-seeking behaviors (Yun et al., 2004;Ikemoto, 2007;Lammel et al., 2011Lammel et al., , 2014. Therefore, our identification of Neurod6 as a marker for this cell population enables future mechanistic investigations into how this DA subcircuit controls motivated behaviors, the dysfunction of which may be important for the pathophysiology of psychiatric disorders such as drug addiction. It was previously reported that Neurod6 ϩ neurons project to the lateral septum based on the axon projections of NEX-Cre mice and fluorogold retrograde labeling (Khan et al., 2017). However, projections from the VTA to the lateral septum are relatively sparse compared to the striatum and NAc (Beckstead et al., 1979;Swanson, 1982;Beier et al., 2015). To investigate this further, we performed retrobead injections into the dorsal and intermediate regions of the lateral septum and found that relatively few VTA DA neurons projected to this area. The lateral septum therefore represented only 2% of total bead-labeled Neurod6 ϩ DA cells across all injection sites. This indicates that under our conditions, the lateral septum was not a primary projection site of Neurod6 ϩ DA neurons. One possible reason for the discrepancy is that the Khan et al. (2017) study relied exclusively on the NEX-Cre mouse line, as opposed to endogenous Neu-rod6 expression as done here. We found that a substantial proportion (37%) of the VTA neurons labeled in NEX-Cre mice are non-dopaminergic and that many of the VTA neurons projecting to the lateral septum are also nondopaminergic (37%), which may have contributed to the differing results.
A defining feature of DA neurons is their susceptibility or resilience to degeneration in the context of PD (Damier et al., 1999;Alberico et al., 2015). As such, significant effort has been made to identify the molecules that determine vulnerability, as this has clear clinical importance (Brichta and Greengard, 2014). Interestingly, in addition to being genes that define specific DA neuron subtypes, both Neurod6 and Grp have been shown to confer neu-roprotection in cell culture models (Chung et al., 2005;Uittenbogaard and Chiaramello, 2005;Uittenbogaard et al., 2010;Baxter et al., 2012). To determine if the neuronal populations defined by Neurod6 and/or Grp are preferentially spared in the context of PD, we performed unilateral 6-OHDA injections and compared the relative abundance of these markers in the 6-OHDA versus control hemisphere. We found that expression of either of these genes alone was not sufficient to confer neuroprotection, as Grp ϩ neurons in the SNc degenerated completely and Grp ϩ neurons in the VTA showed either similar (Grp ϩ / Neurod6 ϩ ) or greater (Grp ϩ /Neurod6-) susceptibility relative to all VTA DA neurons. Notably, we did find that the small population of Neurod6 ϩ VTA DA neurons that lack Grp was significantly spared compared to neighboring DA neurons. This finding is consistent with a potential neuroprotective effect of Neurod6 but indicates that other factors are likely involved, since Neurod6 ϩ VTA DA neurons that coexpressed Grp were not preferentially spared. Together, these results refine our understanding of the genetic factors contributing to vulnerable and spared DA cell types and suggest that the combinatorial expression of genes in a given cell population is important for defining vulnerability.
In summary, our work provides in-depth characterization of Neurod6 and Grp expression in the midbrain and reveals previously unappreciated complexity in how these markers define specific DA subpopulations. Our results provide new insights into the genetic and functional heterogeneity of the DA system, which is just beginning to be unraveled. Future studies can use this information to design intersectional genetic tools based on the expression of two or more genes that will allow access to specific dopaminergic subpopulations. These types of tools represent powerful approaches to dissecting the complex ways in which the DA system contributes to behavior and disease.