Introduction

Midbrain DA neurons are traditionally organized in three main anatomical areas – the ventral tegmental area (VTA), substantia nigra pars compacta (SNc) and retrorubral area (RR)¹, with further sub-divisions therein². To explain the disparate roles of DA in normal behavior and disease, additional heterogeneity of midbrain neurons within these areas has been postulated³. Recent studies have indeed revealed diversity within clusters in electrophysiological properties as well as responses during various behavioral paradigms^4–15. Complementing these studies, recent evidence using single cell classification has opened the possibility that DA neurons can be clustered based on their molecular signatures^16–25. Early studies have begun to suggest that molecularly distinct DA populations may have distinctive anatomical projection patterns, as well as functionally distinct activity patterns^{4, 25–27}. For example, in the SNc, a ventral Anxa1+ population is enriched in acceleration correlated activity, whereas dorsal Calb1+ or lateral Vglut2+ populations are enriched in deceleration correlated activity²⁵. These results suggest that molecularly defined DA neurons need to be considered as a key variable for electrophysiological and behavioral studies.

Based on these exciting results, and the fact that the functional interrogation of DA neurons is happening at rapid speed often agnostic to DA subpopulations, there is a clear need to define DA neuron subtypes with more granularity. Previous studies using single cell sequencing in mice have been limited by the number of cells analyzed, the quality of cDNA library, as well as bias during isolation, thus providing an incomplete picture of DA neuron subtypes^{17, 19–23}. Furthermore, the anatomical distribution of these subtypes in the midbrain remains only partially elucidated. Additionally, given the closely related nature of these subtypes, the stability of these subtypes across pathological conditions remains unclear. In other words, are some of these “subtypes” a representation of cell state rather than cell type?

The unbiased study of the transcriptomic landscape of DA neuron subtypes in the midbrain can inform about the molecular drivers of selective dysfunction of DA neurons in PD. Pathogenic mutations which increase leucine-rich repeat kinase 2 (LRRK2) activity are one of the most common causes of autosomal dominant PD and clinically similar to sporadic cases^{28, 29,30}. Additionally, even in patients with idiopathic PD, LRRK2 kinase activity is increased in DA neurons³¹, providing justification for studying LRRK2-driven mechanisms in DA neurons. Structural and functional synaptic changes are a recurring theme with LRRK2 mutations^32–34. Accordingly, impaired DA transmission is observed with several LRRK2 mouse models including knock-in (KI) mouse lines expressing the most common LRRK2 mutation (G2019S) at physiological levels^35–37, suggestive of functional synaptic deficits in dopaminergic circuits in the absence of overt DA neuron loss. Thus, LRRK2 KI mouse models are valuable for investigating early PD mechanisms in vulnerable DA neurons. However, despite the clinical relevance, the cell-autonomous functions of LRRK2 in nigral DA neurons are largely unknown.

Single-nuclear RNA sequencing (snRNA-seq) offers an avenue to begin to understand downstream transcriptomic alterations within DA neuron subtypes. We have used snRNA-seq to surveil DA neurons, which allows greater acquisition of DA numbers, and minimizes isolation bias²⁵. Here, we extend previous studies by optimizing our isolation protocol, enabling sequencing of large numbers of DA nuclei. We provide a high-resolution view of DA subtypes and an online portal to interrogate this dataset. We systematically map DA subtype distribution using MERFISH, providing unprecedented spatial resolution. Finally, we determine the stability of these subtypes in a preclinical mouse model of PD, LRRK2^G2019S, and reveal molecular pathways that are altered in locomotion-relevant DA neuron subtypes.

Results

A large snRNA-seq dataset of midbrain dopaminergic neurons reveals novel molecular subtypes

We performed snRNA-seq from dopaminergic nuclei (defined by Dat-Cre expression) in LRRK2^G2019S heterozygous mutants and control littermates using a similar technique to our previously published dataset²⁵ but modified to utilize a chip-based microfluidic sorting method which minimizes nuclear stress. This improved protocol enabled increased yield and sample quality (Fig 1A, Fig S1A). Doing so enabled us to generate a dataset with over 28,500 profiled cells, after quality control filtering, with a median UMI count of 7750.5 and median of 3056 genes. Downstream clustering resulted in 22 unique clusters (Fig 1B), of which two clusters (16 and 21) were determined to represent possible dopaminergic/glial doublets based on co-expression of glial markers in these populations and were excluded from downstream analyses (shown in grey in Fig 1B, Fig S1F). The 20 remaining clusters all showed expression of key pan-dopaminergic markers including Slc18a2 (Vmat2), Ddc, Th, and Slc6a3 (Dat) (Fig 1C). These clusters were observed in both males and females, and genotypes, in similar proportions (Fig S1B-E). These subtypes showed distinct expression signatures of many previously described markers of DA neuron subtypes (Fig 1D)¹⁷. However, our increased yield and sample quality also allowed explicit detection of known populations that have eluded most prior single-cell studies, such as DA neurons expressing Vip^{17, 38}, thereby demonstrating the strength of this pipeline in detecting rare DA neuron subtypes. Importantly, all the described subtypes showed high expression of midbrain floorplate markers (Fig 1E), with exception of clusters 0 and 17, which possibly represent Slc6a3+ neurons in the supramammillary and/or premamillary regions^{39, 40}. Thus, our increased number of subtypes over prior studies represents increased granularity among classic midbrain DA neuron subtypes, rather than detection of non-floorplate derived populations that express DA markers.

Figure 1.

Analysis of differentially expressed genes reveals families of subtypes and relationships to prior studies

Given the large number of clusters discovered in our dataset, we first sought to organize these subtypes into groups in an unbiased manner. By creating a cluster dendrogram (Fig 2A), we were able to map an approximation of the relative relatedness of different subtypes, which segregated into three main families based on major branch points. By running differential expression on each node of the cluster dendrogram (i.e., exploring differentially expressed genes (DEGs) between the two arms emerging from a given node), we were then able to create a stepwise expression code that can uniquely define any individual cluster (Fig 2A, G). The earliest division point among our clusters was defined by genes Prkcd and Lef1, which are highly expressed only in cluster 17 (Log2 fold changes = 6.279 and 7.246 respectively, BH-corrected p values ≈ 0). Next, a family of clusters defined by Gad2/Kctd16 (deemed Gad2 family, Log2 fold changes = 5.318 and 4.014 respectively, BH-corrected p values ≈ 0) separated from all other clusters, which expressed significantly higher levels of Ddc and Slc6a3 (Log2 fold changes = 1.419 and 1.643 respectively, BH-corrected p values ≈ 0). The first major branch point between Ddc-high/Slc6a3-high clusters was then defined by cells expressing Sox6/Col25a1 (Log2 fold changes = 2.720 and 2.305 respectively, BH-corrected p values ≈ 0) versus those expressing Calb1/Ndst3 (Log2 fold changes = 2.025 and 3.368 respectively, BH-corrected p values ≈ 0), thus creating the Sox6 family and Calb1 families respectively. This division between Sox6 and Calb1 is consistent with our previous work³⁸, as well as a prior large single cell profiling study that found this division to be the central branch point among DA neurons across several species including humans ⁴¹. Cluster division can be seen by visualizing coexpression of the DEGs, taken from the branch points labeled in red (Fig 2C-F), on the UMAP plot. Of note, while these branch points markers describe the overall trends of gene expression, these genes are not perfectly distinct or universal among the members of these cluster families. For example, clusters 1, 11 and 13 show expression of both Sox6 and Calb1 (Fig 1D). The co-expression of these markers has previously been described in a fraction of mouse DA neurons^{38, 42, 43}. Furthermore, because the calculations of differential expression at each node take into account only the descendants of that branch point, markers defining a branch point may also be strongly expressed in clusters outside this comparison. Thus, while the final marker(s) listed for each subtype in Fig 2A are not unique to that population, following the sequential branches leading to this marker in a stepwise fashion will specifically define this population (Fig 2G). By doing so, we provide a strategy for potential genetic access to any individual DA subtype using intersectional logic gates.

Figure 2.

To enable comparisons of DA populations across existing and future literature, we next sought to better define genetic signatures for each of the individual putative subtypes and to provide a common language for discussing these populations. To do so, we utilized two complementary approaches. First, we explored the top differentially expressed genes with positive relative expression in each cluster (Fig 2B). While many of these markers are expressed in more than one population, some (e.g. Vip) are almost entirely unique to a single cluster. As a complementary approach, we performed differential expression among clusters within each of our three subtype families, thereby creating a shorthand identity based on the cluster family and top DEGs (Fig 2G). Within the Sox6 family, these clusters are Sox6 ^Tmem132d (listed as cluster 1 in Fig 1 and 2), Sox6^Arhgap28 (cluster 2), Sox6^March3 (cluster 3), Sox6^Tafa1 (cluster 4), Sox6^Kcnmb2 (cluster 8), and Sox6^Vcan (cluster 10); in the Calb1 family, these are Calb1^Pde11a (cluster 5), Calb1^Kctd8 (cluster 6), Calb1^Ptprt (cluster 7), Calb1^Sulf1 (cluster 9), Calb1^Stac (cluster 11), Calb1^Chrm2 (cluster 12), Calb1^Sox6 (cluster 13), Calb1^Lpar1 (cluster 14), Calb1^Ccdc192 (cluster 18), Calb1^Gipr (cluster 19); in the Gad2 family, these are Gad2^Syndig1 (cluster 0), Gad2^Egfr (cluster 15), and Gad2^Ebf2 (cluster 20). Cluster 17 (defined by Lef1) did not have an associated cluster family (Fig 2A). These shorthand identities provide a simple nomenclature for our clusters, which we have utilized throughout the remainder of this paper.

After establishing genetic signatures for each putative subtype, we were next able to correlate our cluster identities to those described in another recent study²⁵, as well as the consensus subtypes proposed in Poulin et al, 2020 (Fig 2G)¹⁷. Doing so provides context for discussing our putative DA neuron subtypes with regards to prior literature. For example, cluster Sox6^Tafa1 was found to be equivalent to a population we previously defined by high expression of Anxa1, which hold notably distinct functional and anatomical properties²⁵. Several of the subtypes that emerged in our new dataset also appear to be novel divisions within clusters that previously showed evidence of internal heterogeneity based on their cluster stability metrics²⁵. For instance, a Vip+ cluster (Calb1^Gipr) emerged from within a highly variable cluster (defined as #10 in ²⁵), as well as additional subdivisions emerged from within Sox6+ clusters (most notably #3 described in ²⁵) (Fig 2G). This emergence of unique populations, such as Calb1^Gipr, is likely due to the increased yield and sample quality of our new dataset which has enabled better detection of small populations and gene markers with lower expression level. This is further supported by the high cluster stability of the new Calb1^Gipr (Fig S2A,B), resolving the instability of the equivalent cluster in Azcorra & Gaertner et al²⁵. We reclustered our snRNA-seq data into the clusters defined in Azcorra & Gaertner et al. and found our clusters aligned with or further subdivided many of the previously defined clusters (Fig S3A). Lastly, to further corroborate our cluster identities, we used an independent bioinformatic platform, Scanpy, to re-cluster our nuclei (Fig S3B). Comparing the Seurat and Scanpy clusters using a Sankey diagram (Fig S3C) shows a high correlation between the Seurat and Scanpy clusters.

Mapping the spatial location of midbrain DA neurons using MERFISH

Having identified multiple clusters of midbrain DA neurons, we next resolved to determine the spatial distribution of these populations using MERFISH⁴⁴, a form of spatial transcriptomics that offers subcellular resolution imaging-based RNA quantification. We assembled a 500 gene panel which included 120 of the top DEGs identified within the snRNA-seq dataset, as well as 380 general markers of neuronal and non-neuronal cell populations (Supp Table 1). Seven coronal slices were collected from adult mouse midbrains and processed for MERFISH (Fig 3A). The coronal slices were selected containing ventral midbrain regions between -2.9mm to -3.9mm from bregma. Cells were segmented using the CellPose algorithm based on DAPI and Poly-T staining and passed through quality control filters to remove cells with too large of a volume (potential doublets) or too few transcripts detected (false cells; Fig S4A-C). This resulted in 429,713 identified cells (Fig 3B). Cells from the seven sections displayed significantly correlated gene expression and minimal batch effects (Fig S4D). Cells were integrated in Seurat using reciprocal principal component analysis (RPCA) algorithm and then clustered, resulting in 40 final clusters comprising distinct classes of neurons, oligodendrocytes, oligodendrocyte precursor cells (OPCs), astrocytes, microglia and meningeal cells (Fig 3B). These clusters could be attributed to 25 functional groups based on the expression of cell type markers from Supp. Table 1 and their spatial distribution. Glial cell associated genes were mostly contained within their expected cell type: oligodendrocytes (Mog, Cnp), OPC (Pgfra, Olig2), astrocytes (Aqp4, Gfap), and microglia (Aif1, Tmem119). However, we noted the unexpected distribution of the oligodendrocyte-specific transcript Mbp in non-glial cells which we believe might indicate the presence of this mRNA in glial-processes that overlap non-glial cells⁴⁵. Within the neuronal class we identified 16 excitatory clusters (presence of Slc17a6 or Slc17a7), 7 inhibitory clusters (presence of Slc32a1, Gad1 or Gad2), one serotonergic cluster and one DAergic cluster (Fig 3B-D, S4G). Though our gene panel was designed to identify neuronal heterogeneity of the midbrain, it offered the resolution to delineate 7 excitatory neuron classes present in distinct layers of the cortex and 6 cortical interneuron classes (Fig S4E,F). Cortical excitatory neurons were absent from subcortical structures such as thalamus, hypothalamus, midbrain and hindbrain. This fundamental dichotomy between cortical and subcortical neurons in the adult mouse brain is supported by recent spatial transcriptomic studies^46–48.

Figure 3.

A population of DA neurons was identified among our 40 clusters based on the expression Th, Slc6a3 and Slc18a2 (Fig 3C). The location of the dopaminergic cluster (blue in Fig 3D) considerably overlaps with Th expression (Fig 3E). This cluster comprised 4,532 cells, was entirely located in the midbrain, and was also enriched for Ddc, Dlk1, and Drd2 (not shown). Within this population, we detected expression of genes associated with glutamatergic and GABAergic neurotransmission, in line with previous reports of multi-lingual characteristics of DA neurons^49–54, for instance the glutamatergic marker Slc17a6 (Vglut2) and GABAergic markers Gad2 and Slc32a1 (Vgat) (Fig S4G). Markers for cholinergic (Chat) and serotonergic (Tph2) neurotransmission were not detected in putative DA neurons and the glutamatergic transporters Slc17a7 (Vglut1) and Slc17a8 (Vglut3) were also absent. We next wanted to characterize the neuroanatomical location of this DAergic cluster. We manually delineated neuroanatomical boundaries on each of the slices for seven midbrain regions known to contain most dopamine neurons including SNc, VTA, retrorubral area (RR), caudal linear nucleus (CLi), interfascicular nucleus (IF), rostral linear nucleus (RLi), and periaqueductal grey (PAG). These boundaries were delineated based on cell density (dapi), white matter landmarks, and Th transcripts distribution using the Allen Brain Atlas as a reference (Fig S5A). We found that cells in our DA neuron cluster were located in the following locations: 35.2% in the SNc, 19.8% in the VTA, 7.8% in the CLi, 15.3% in the RR, 2.0% in the IF, 0.6% in the RLi, 4% in the PAG and 15.3% located elsewhere (Fig S5B).

Finally, to explore the diversity of DAergic neurons identified by MERFISH we subclustered the 4,532 cells to resolve 12 distinct subgroups (Fig S5C-F). For clarity, we are designating the clusters identified solely with MERFISH transcript levels as MER-0 to 11. Many of the clusters were characterized by increased expression of known subtype markers such as Calb1 (MER-6, MER-8 and MER-10), Aldh1a1 (MER-5), Ndnf (MER-0), Gad2 (MER-7 and MER-9) and Slc17a6 (MER-2).

Cell label transfer allows the mapping of snRNA-seq clusters

To provide the spatial distribution of the clusters identified using snRNA-seq, we sought to project the MERFISH dataset onto the snRNA-seq UMAP space. We first compared the expression level of the 500 genes utilized in our panel to select features for integration with comparable expression levels. We excluded 128 genes based on the deviation of average counts per cell within the datasets, potentially resulting from faulty probes or improper assignment of RNA transcripts to segmented cells. For instance, Sox6 transcripts were not detected across the entire brain (data not shown). We successfully projected the MERFISH dataset onto the snRNA-seq UMAP space and identified 2,297 dopamine neurons with a cell similarity score > 0.5 increasing confidence in the correct assignment to snRNA-seq defined clusters. Almost all clusters were represented by more than 20 cells in our MERFISH dataset (Fig. 4A and S6A). Clusters Calb1^Ccdc192 and Calb1^Gipr were not resolved in this label transfer either due to their relatively smaller number or location in an under-sampled region. Indeed, the Vip-expressing Calb1^Gipr has been shown to be located in the caudal midbrain (PAG/DR)³⁸, a region not well covered in our MERFISH experiment. Comparing snRNA-seq clusters with the ones obtained by clustering the MERFISH dataset alone showed a strong correspondence (Fig S6B). Half of the snRNA-seq clusters were predominantly comprised of single MERFISH clusters (10/20) and 25% (5/20) were comprised of 2 MERFISH clusters, demonstrating that the integration with the snRNA-seq dataset allowed us to further refine our MERFISH-based classification (Fig S6B). We then imputed expression data from the snRNA-seq dataset onto the MERFISH cells to visualize whole genome expression levels. The imputed data across the genes used for data integration correlated strongly with MERFISH transcripts as shown for Calb1, Gad2 and Aldh1a1 genes (Fig S6C-E). Altogether, the two datasets integration allowed us to map the neuroanatomical location of 18 of the clusters identified by snRNA-seq as well as impute expression of the whole transcriptome.

Figure 4.

Locating DA neuron clusters within the midbrain

Clusters were sorted into Gad2, Sox6 and Calb1 families based on our hierarchical dendrogram (Fig 4A; Gad2 n=124 cells, Sox6 n=887 cells, Calb1 n=590 cells). Imputed Gad2, Sox6 and Calb1 gene expression within DA neurons showed a distinct medio-lateral and rostro-caudal distribution (Fig 4C-E). The spatial distribution of imputed Sox6 and Calb1 transcripts within DA neurons correlated well with previous reports^38,43,55 whereas the distribution of the Gad2 family is predominantly observed in the rostral linear/posterior hypothalamus region and CLi, and to a lesser extent in the VTA (Fig 4B, S7). Each family was overrepresented in a neuroanatomically-defined area; for instance the SNc is composed of 78.4% of neurons from the Sox6 family (Fig S7A). However, each family was also found across multiple neuroanatomical regions (Fig S7A). We also observed spatial heterogeneity within families, for example, one branch of the Calb1 family dendrogram (Fig 4, label in shades of yellow-orange defined by higher Lepr expression) was enriched in the dorsal VTA whereas the other Calb1 branch (Fig 4, label in shades of red defined by higher Ntm expression) populated mostly the ventromedial VTA. Similar observations could be made for the Sox6 family, where subtypes displayed biased distributions. This led us to explore the distribution of each individual cluster and representative markers toward the goal of generating a granular map of the DAergic system.

The hierarchical clustering of the Sox6 family reveals two major branches (Fig 2A and 4A): 1) Branch 1 defined by higher levels of Slc44a5 and composed of subtypes Sox6^Tmem132d and Sox6^March3 and 2) Branch 2 defined by Cntnap5 and composed of subtypes Sox6^Tafa1, Sox6^Arhgap28, Sox6^Kcnmb2 and Sox6^Vcan. We observed an important discrepancy between these two branches, with branch 1 being located more dorsal and medial than branch 2 (Fig 5A). Interestingly, this dorsoventral organization breaks down in the caudal SNc and RR where neurons from both branches are mostly intermixed. Only minor differences in location were observed between clusters of Branch 1, maybe reflecting their molecular similarities. Of note, whereas Sox6^Tmem132d is predominant in the rostral SNc and VTA for this branch, neurons of Sox6^March3 are more numerous in the RR (Fig 5B, Fig S7B). Within Branch 2, the molecular identity correlated with the medio-lateral distribution. For instance, neurons of cluster Sox6^Arhgap28 were found to be more lateral in the SNc compared to neurons of Sox6^Tafa1, Sox6^Kcnmb2 or Sox6^Vcan. This medio-lateral distribution of neurons of the SNc ventral tier might reflect the topographical nigrostriatal projections previously reported²⁶.

Figure 5.

The Calb1 family was the most molecularly diverse which was also reflected in its localization (Fig 6). Two major branches of Calb1 DA neurons were identified: 1) Branch #1 (Fig 6A-B) containing subtypes Calb1^Ptprt, Calb1^Sulf1, and Calb1^Sox6 and 2) Branch #2 (Fig 6C-D) containing subtypes Calb1^Chrm2, Calb1^Kctd8, Calb1^Stac, Calb1^Pde11a and Calb1^Lpar1. Broadly speaking, Calb1 subtypes were predominantly localized to the VTA and CLi with some neurons showing localization to the dorsal SNc, SNpl, or RR (Fig S7). Within the VTA, Branch #1 Calb1 neurons tended to localize to the dorsal VTA (Fig 6A-B) while Branch #2 Calb1 neurons showed distinct ventromedial localization (Fig 6C-D). In more caudal sections the difference was even more pronounced, although intermingling was always observed. Within branch #1, subtype Calb1^Ptprt and Calb1^Sox6 were largely absent from the SNc dorsal tier and RR, whereas subtype Calb1^Sulf1 was present. Some notable distinctions were observed in both groups, specifically subtype Calb1^Sulf1 and Calb1^Lpar1 which were each observed in significant numbers in the dorsolateral SNc/SNpl region. Of these, only Calb1^Lpar1 expressed significant Slc17a6 (data not shown), making it a likely candidate for the SNpl Vglut2+ DA neuron population found to innervate the tail of the striatum²⁷. Indeed, these laterally located neurons have been genetically targeted with a Slc17a6-Cre driver, and display deceleration correlated responses, and robust responses to physically aversive stimuli^{25, 27}.

Figure 6.

Within the Gad2 family we were able to resolve two subtypes. Subtypes Gad2^Egfr and Gad2^Ebf2 had very similar spatial distributions and were prominent in the posterior IF, RLi, and CLi (Fig 7, Fig S7). These regions have been associated with small DA neurons of lower Th and Slc6a3 levels. Subtype Gad2^Egfr neurons were more plentiful than cluster Gad2^Ebf2, with the latter being absent in more caudal sections. All Gad2 subtypes expressed Slc32a1, opening the possibility that these neurons co-release GABA. Further, they also expressed Slc17a6, revealing a potentially multilingual DA neuron subtype.

Figure 7.

LRRK2^G2019S KI results in pan-dopaminergic gene expression changes consistent with deficits observed in this model, without altering DA subtype proportions

With cluster identities and locations defined, we next sought to explore the downstream gene expression characteristics of DA neurons in Lrrk2 mutants. The role of Lrrk2 in PD pathophysiology remains unclear⁵⁶, and Lrrk2 mutations are often proposed as being an indirect source of dysfunction for the DA system, such as through its high expression in glial or striatal cells^{34, 57–60}. However, previous studies have shown alterations in DA signaling or gene expression even when the mutant LRRK2 was expressed specifically in DA neurons^{61, 62}. Thus, as a first step in demonstrating the possibility that mutant Lrrk2 in DA neurons themselves can contribute to dysfunction, we mapped the expression of Lrrk2 across all DA neuron subtypes (Fig 8A). We found detectable expression across most subtypes, with the notable exception of Gad2^Egfr neurons.

Figure 8.

Given the distinctive properties of different DA neuron subtypes, we next sought to test if phenotypic changes in Lrrk2 mutant mice are driven by changes in the relative proportions of said subtypes. While previous studies have shown no overt DA neuron loss, it is possible that there was a selective reduction of a subtype(s) that would be undetectable by previous methods. We first graphed the proportions of each subtype within our control and mutant datasets (Fig 8B), and found remarkable similarity across these samples, consistent with previous reports showing no DA neuron loss in these mice³⁶. To further assess for changes in cell state among clusters across conditions, we next compared the similarity of each cluster to all other clusters across conditions (Fig 8C). Meta-neighbor analysis revealed that the most similar population to each wildtype cluster was its corresponding subtype in the Lrrk2 mutant samples. This suggests that the general subtype organization is impervious to LRRK2^G2019S perturbation, and that our clustering scheme is more likely representative of cell type rather than cell state.

To assess pan-dopaminergic gene expression changes as a function of Lrrk2 genotype, we next looked at globally differentially expressed genes across our mutant and control samples. Calculating differentially expressed genes across all clusters revealed a large number of significantly enriched or diminished genes as a function of Lrrk2 mutation (1646 genes in total with BH corrected p-value <0.05, Wilcoxon rank sum test). The results of this differential expression analysis are plotted in Fig 8D, with full results and statistics available in Extended Data. Notably, most of these changes appear to be small in magnitude; among FDR-significant genes, only 51 showed increased expression and 7 showed decreased expression with log2 fold changes greater than 0.5 (Fig 8D). Two genes with high enrichment in Lrrk2 mutants, Mir124a-1hg and AC149090.1 were both of particular interest due to their potential relevance to PD pathogenesis. AC149090.1 is orthologous to the human PISD gene, which is involved in autophagy and implicated in mitochondrial dysfunction^63–65. This gene has been proposed as a sensitive marker of biological aging in the brain^{63, 66}, an interesting finding given that age is the single largest risk factor for PD. Also of note, Mir124a-1hg is a host gene (though not the exclusive gene for) miR-124, a microRNA with purported neuroprotective effects that has been found to be downregulated in PD patients or some PD mouse models^67–70.

We next performed pathway enrichment analyses using gene set enrichment analysis (GSEA) and gene ontology (GO) to compare our Lrrk2 mutant samples to our controls. Using GSEA, we found 18 gene sets with significant differential expression via GSEA (FDR adjusted p-values < 0.05). The top 13 of these enriched pathways (those with NES scores > 1.8) are shown in Fig 8E. While a broad range of biological domains are found among these results, many of the top enriched pathways were related to energy production via mitochondrial oxidative respiration, remarkably similar to prior comparisons between vulnerable and resistant DA neuron populations⁴³. Mitochondrial energy production pathways, and specifically the extraordinary bioenergetic burdens of DA neurons, have been heavily implicated in PD vulnerability^{71, 72}. Thus, the enrichment of these pathways in Lrrk2 mutants is particularly intriguing as it may reflect the structural and functional alterations reported with LRRK2 mutations in patients⁷³ and preclinical models^{74, 75} and further predispose these cells to dysfunction and ultimately degeneration. Utilizing GO, we again found many pathways displaying significant differential expression between Lrrk2 mutants and controls. Among these, the top pathways are all associated with synaptic organization and function (Fig 8F). This was particularly intriguing given that pre-synaptic dysfunction in Lrrk2 mutants such as endocytosis^76–78 or trafficking of axonal cargo impairments⁷⁹ could underlie the mechanisms for the deficits in DA signaling in these mice^{80, 81}. Of note, decreased endocytosis was observed in DA and not in cortical or hippocampal LRRK2^G2019S primary neurons⁸². These results, although in vitro, suggest cell type specific synaptic effects of dopamine neurons, consistent with their vulnerability in LRRK2-related PD.

Comparing individual subtypes across conditions allows insights into subtype-specific dysfunction

Given that global gene expression differences across our conditions appear to hold relevance for putative mechanisms of PD pathogenesis, we next sought to address the question of whether gene expression within any given subtype is specifically altered in our PD model. Although Lrrk2 mRNA seemed uniform across subtypes, downstream effects could be subtype-specific secondary to differential Lrrk2 kinase activity, distinct kinase targets, or superimposed differences in other properties of each subtype that culminate in distinct downstream transcriptional alterations. Given the proposed vulnerability of Sox6+ DA neurons and relative resilience of Calb1+ neurons in PD^{3, 41, 43}, we first compared these two cluster families across Lrrk2 mutants versus control mice to see if any changes in DA neurons are specific to vulnerable or resistant cell types. We found 729 DEGs among the Sox6 family, and 679 among the Calb1 family (BH corrected P value <0.05). Among those genes, we found 327 DEGs were shared between the Sox6 and Calb1 families (Fig S8A). The full list of DEGs and statistics is available in Extended Data. Notably, AC149090.1 was highly significant in both cohorts, implying that the pan-DA differential expression of this gene was not driven by a particular subset of DA neurons. To further compare changes within the Sox6 family, we once again performed GSEA and GO but on these groups of clusters in isolation (Fig 9A-B). In line with our previous results, GO revealed pathways that largely corresponded to synaptic and axonal function. However, the pathway with second largest enrichment now corresponded to locomotory behavior (Fig 9A). Finally, to extend the granularity our analyses, we applied the same pipelines to one individual cluster, Sox6^Tafa1, which has the highest expression of Anxa1. This SNc subtype was of particular interest to us given recent results showing that Anxa1+ ventral tier SNc neuronal activity is selectively correlated with acceleration in mice running on a treadmill, leading to the hypothesis that degeneration of these neurons may, in part, underpin PD motor dysfunction²⁵. The top enriched GO pathway in our analysis was locomotory behavior (FDR-corrected p value = 3.91E-6), with several significant pathways once again corresponding to synaptic function.

Figure 9.

Next, we compared the synaptic compartment and biological functions of DEGs (BH adjusted p values < 0.05) in Lrrk2^G2019S mutants versus controls using SynGO, a database for systematic annotation of synaptic genes⁸³. We focused primarily on the Sox6 and Calb1 families, but also on cluster Sox6^Tafa1 given that this population (SNc neurons with highest expression of Anxa1) is specifically linked to movement²⁵. Of the DEGs in the Sox6 family, Calb1 family, and Sox6^Tafa1 cluster, 24.82%, 20.19% and 29.41% of these genes had synaptic localizations, respectively, with more significant synaptic associations among populations localized to ventral tier SNc (Fig 9D-F, Fig S8B). Further compartment categorization revealed that among these synaptic DEGs annotated by SynGO, 45.83% of Sox6 family, 64.0% of Sox6^Tafa1 and 48.82% of Calb1 family DEGs have presynaptic localizations, raising the possibility of enhanced disruption of presynaptic function in acceleration-correlated DA neurons. From these presynaptic DEGs, 20/77 in the Sox6 family, 5/16in Sox6^Tafa1, and 16/62 in the Calb1 family have specific annotations for active zones, which are the site of dopamine release. Annotations based on the biological functions showed that 20.83%,28.0% and 15.74% of Sox6, Sox6^Tafa1 and Calb1 DEGs, respectively, are associated with processes in the presynapse (Fig S8C). Overall, greater proportions of DEGs are associated with presynaptic locations in cells from vulnerable DA neurons (Sox6 family, and in particular, Sox6^Tafa1), compared to less vulnerable ones (Calb1 family).

Given the differential effects of Lrrk2 mutation across subtypes, we also sought to map the association of gene expression in our dataset with co-expression of PD-associated risk loci identified by GWAS⁸⁴ using a polygenic enrichment score method, scDRS⁸⁵. We first calculated a PD GWAS risk score for each individual cell, which revealed clear bias for higher risk scores among clusters that localize to the ventral SNc (Fig 9G). Calculating the mean PD GWAS risk score for each cluster and corresponding 95% confidence intervals confirmed significant risk associations in several clusters, particularly those in the Sox6 family (Fig 9H). Interestingly, the Gad2 family of clusters all showed significant negative associations with PD risk loci (Fig 9H). Of note, the highest scDRS scores were in clusters Sox6^Tafa1 and Sox6^Vcan, the only two SNc clusters that express Anxa1. By calculating the average scores for each cluster family as a whole, as well as these two Anxa1-expressing clusters, we found that only the Sox6 family (and Anxa1-expressing SNc clusters within it) show significant associations with PD risk loci (Fig 9I).

Finally, due to the large number of potential comparisons to be made across clusters, cluster families, and genotypes, we have also developed an online tool that allows for exploration of our datasets, which we have termed Dopabase (URL: https://awatramanilab.shinyapps.io/shinyapp-2/). Utilizing this tool, users have the ability to access many analysis and plotting functions to continue to explore populations of interest.

Discussion

Our work provides five key advances. 1. We provide a large snRNAseq dataset that provides high granularity to our taxonomic schemes. 2. We provide a user-friendly portal to query this dataset. 3. We plot the spatial location of most identified DA subtypes, showing heterogeneity even within sub-domains in traditional anatomical areas. 4. We demonstrate that cluster proportions do not change in LRRK2^G2019S mice, although gene expression changes are observed, most notably in mitochondrial bioenergetics and synapse organization pathways in vulnerable DA subtypes. 5. PD GWAS risk is most prominent in Anxa1 SNc neurons.

We demonstrate the presence of 20 molecularly distinct clusters of midbrain DA neurons. These are hierarchically organized, the dendrogram being split into three levels, which we refer to as family, branch and subtype. Family (Level 1) distinctions are driven by the expression of Gad2/Kctd16 vs Slc6a3/Ddc high expressing cells. Slc6a3/Ddc high cells are further divided by Sox6/Col25a1 and Calb1/Ndst3 (Level 2) among other genes. It is likely that these fundamental divisions are, at least in part, set early in development in the mesodiencephalic floor plate, where Sox6 demarcates progenitors with a substantially higher probability of ventral SNc and lateral VTA fate, and conversely, Sox6-progenitors have a higher probability of a Calb1+ fate^{43, 55}. The Sox6 family is split into two branches with six subtypes, the Calb1 family is split into two branches with ten subtypes, whereas the Gad2 family can be further split into two subtypes. This scheme represents a more granular taxonomy than previous studies^{19–21, 23, 25, 38, 41}, and bears some resemblance to a recent study⁸⁶. Key congruencies between these studies include: 1) parsing midbrain DAergic system into broad groups of subtypes based on expression of Sox6 and Calb1, and to a lesser extent Gad2, populations, 2) greater diversity observed in Calb1+ clusters, 3) the identification of some consensus and homogenous subtypes. Key discrepancies include the extent to which each cluster further subdivided and the final number and molecular signature of DA subtypes. In the latter study⁸⁶, they refer to levels as territories and neighborhoods. For simplicity of communication, we refer to the lowest level of our dendrogram of molecularly defined DA neuron clusters, as subtypes. However, we acknowledge that because of the common developmental origin from the floor plate^87–90, their high relatedness, the blurriness of boundaries between clusters, and the polythetic nature of clustering, referring to these groups as neighborhoods is an alternative and reasonable possibility. Notwithstanding these considerations, defining markers of some of these clusters and generating Cre/Flp lines, may be useful to at very least enrich for DA neuron populations with distinctive properties as has recently been shown for the Anxa1+ SNc neurons, which at a population level are selectively correlated with acceleration, while reward responses are almost entirely absent.

Our spatial MERFISH analysis provides the most granular representation of DA subtype locations to date. This resolution was obtained by the integration of snRNA-Seq and MERFISH datasets. Our data complement current spatial transcriptomic datasets^46–48. For instance, using MERFISH, Zishen et al. identified 10 cluster in the VTA (A10), 9 clusters in the SNc (A9), and 2 clusters in the RR although little detail on the molecular identity or location of these subtypes is provided. In addition, Langlieb et al ⁴⁸, provide a molecular atlas of the adult mouse brain in which they identified 13 clusters of midbrain DA neurons and located these using Slide-Seq. By contrast, Salmani et al.⁸⁶ used established markers to locate seven DA neuron territories and 16 neighborhoods. They report a distribution of Gad2+ cells in the midline including CLi, RLi and IF, which bears similarity with our Gad2 family. They also Pdia5+/Calb1+ DA neurons in the pars lateralis, likely representing an overlapping population with our Calb1^Postn and Calb1^Sulf1. Within the SNc, our data fits also well with the location of DA neurons that could be labeled with Anxa1-Cre, although here we show that in our Seurat based clustering, Anxa1 expression is detected in Sox6^Tafa1 and to a lesser extent in Sox6^Vcan cells. Taken together, with our increasing appreciation of DA neuron heterogeneity, a consensus is bound to arise on both the identity and location of DA neuron subtypes. Integration of all these emerging datasets will be necessary to achieve a comprehensive DAergic neuron taxonomy that will ultimately encompass all morphological, connectional, physiological, molecular, and neuroanatomical parameters.

A key finding of this study is the heterogeneity observed even within sub-domains of the traditional anatomical clusters. An example of this is the lateral SNc/SNpl. Current models depict the SN heterogeneity as a mediolateral gradient⁹¹. Our data substantially refines this by showing that even within the lateral SNc/SNpl region, there are several distinct subtypes. One interesting finding is that two of the lateral subtypes from distinct families (Sox6 and Calb1), show some common gene expression signatures. For example, Asic2 or Ankfn1 are highly expressed in Calb1^Lpar1 and Sox6^Arhgap28 subtype, both of which are laterally located. One interpretation of this is that secondary gene expression programs may have been superimposed upon the early developmental sub-divisions, for instance during circuit assembly. A second example of heterogeneity within sub-domains is in the dorsolateral VTA. Again, we find enormous heterogeneity in this region, which is unaccounted for in previous models. It is likely that neurons in this region will have vastly distinct anatomical and functional properties, as exemplified by the complex intermixing of neurons with diverging axonal projections⁷. Careful intersectionally based interrogation methods will need to be developed to distinguish these populations.

We address whether our clustering schemes are more reflective of cell types or cell states. We show that the proportions and definitions of subtypes are consistent between controls and mutants, and that mutant subtypes are most closely related to their control counterpart. This suggests that our taxonomic scheme is impervious to the context of a mild perturbation such as LRRK2^G2019S, suggesting that our clusters are reflective of cell types, rather than cell states. It is likely, that with more severe perturbations, such as a toxin lesion, more substantial alterations of taxonomic schemes are observed^{86, 92}. However, we expect that for mild insults, day to day behavioral changes, or pharmacological paradigms, our clusters will be resistant to changes, although individual genes may vary.

Our work describes a rich heterogeneity of DA neurons in the murine SNc that might serve as a reference for understanding patterns of degeneration in PD. The majority of SNc neurons belong to the Sox6 family, but some are from the Calb1 family. Several studies on post-mortem human PD brain have demonstrated a relative resilience of Calb1+ DA neurons^{43, 93, 94}. Calb1+ neurons in mice SNc and SNpl project most densely to the DMS and tail of striatum (TS) respectively^{3, 27}. This could match the patterns of axon preservation observed in post-mortem PD brain, which show relative sparing of medial, ventral and caudal regions of the caudate-putamen^{95, 96}. Indeed it is plausible that patterns of degeneration could be associated with distinct clinical presentations of PD – a recent study suggests that patients with relatively spared caudate projections have a higher probability of tremor⁹⁷. Additionally, we postulate that spared DA neurons projections may be a cellular substrate mediating the side-effects of levodopa particularly at high doses. In line with this, high DA release in the TS is associated with hallucination like behaviors in mice⁹⁸.

Given the relevance of LRRK2 to genetic and idiopathic PD, we examined the effects of LRRK2 on DA neuron gene expression to start appreciating the molecular pathophysiology of relevant subpopulations in PD. We relied on a KI LRRK2 model, expressing the G2019S pathogenic mutation due to the highest prevalence of this mutation over other mutations. Although these mice do not exhibit DA neuron loss, they do have nigrostriatal synaptic alterations within a physiological range³⁵. Our transcriptomic interrogation has three main takeaways. First, the Lrrk2 gene is expressed in most DA neurons in mice; in humans Lrrk2 transcript is enriched in the Sox6 population⁴¹. While previous work has shown the LRRK2 protein expression in midbrain DA neurons^{60, 99}, its relative expression across DA neuron subtypes in mice requires further investigation. Thus, transcriptomic alterations observed within DA neurons are potentially cell-autonomous effects, albeit indirect, since LRRK2 functions as a protein kinase with phosphoregulation primarily affected. Second, while individual genes showed only modest differences, pathway analysis showed a strong dysregulation of synaptic pathways. This is conceptually aligned with prior work describing that LRRK2 phosphorylates various key presynaptic targets (not necessarily in DA axons) associated with SV cycle¹⁰⁰. Impairments of these processes at the nigrostriatal terminals may contribute to the DA release deficits reported in LRRK2 mutants^35–37. Third, alterations in pathways associated with oxidative phosphorylation and energy production are in line with the described LRRK2 impact on mitochondrial structure and mitophagy in the LRRK2^G2019S KI mice. As DA transmission changes and mitochondrial abnormalities are the two main and consistent phenotypes with the KI mouse models, our transcriptomic findings indicate that functional changes are at least partly linked with gene expression changes.

Taken together, our work complements perfectly the recent taxonomy of mouse^46–48 and human brain cell types¹⁰¹. Focusing our efforts on the DAergic system allowed us to reveal a comprehensive transcriptomic signature of midbrain DA neurons as well as ascribe with precision their spatial location. Moreover, we provide a web resource to explore this dataset through multiple angles. We also uncovered cell type specific transcriptomic changes in a prodromal model of PD revealing the molecular consequences of higher LRRK2 activity. Disentangling the molecular diversity of DA neurons, now defining 20 transcriptomic subtypes, raises important questions. For instance, what are the functional implications of this diversity and are these transcriptomic differences reflected in other cellular features? Decades of experiments might be needed to answer such questions, however, by defining genetic anchor points unique to each subtype, we have provided an avenue to target these populations in the mouse.

Methods

snRNASeq

Animals

All animals used in this study were maintained and cared following protocols approved by the Northwestern Animal Care and Use Committee. Cre mouse lines were maintained heterozygous by breeding to wild-type C57BL6 mice (RRID:IMSR_CRL:027), with the exception of DAT-IRES-Cre, CAG-Sun1/sfGFP, which was maintained homozygous and crossed to Lrrk2^G2019S mice (RRID:IMSR_JAX:030961)only for generation of experimental mice for nuclei isolation. Both males and females were used in equal number for all RNAseq experiments.

Sample Preparation

To isolate nuclei for snRNA-seq library generation, n=8 6-month old DAT-IRES-CRE (RRID:IMSR_JAX:027178), CAG-Sun1/sfGFP (RRID:IMSR_JAX:021039), Lrrk2^G2019S+/wt or littermate Lrrk2^wt/wt mice (2 female control, 2 male control, 2 female mutant, 2 male mutant) were sacrificed and rapidly decapitated for extraction of brain tissue as previously described for isolation of GFP+ (i.e. dopaminergic) nuclei. A 2mm thick block of ventral midbrain tissue was dissected out from each mouse and collected for isolation. Tissue was dounce homogenized in a nuclear extraction buffer (10mM Tris, 146mM NaCl, 1mM CaCl2, 21mM MgCl2, 0.1% NP-40, 40u/mL Protector RNAse inhibitor (Roche 3335399001). Dounce homogenizer was washed with 4mL of a washing buffer (10mM Tris, 146mM NaCl, 1mM CaCl2, 21mM MgCl2, 0.01% BSA, 40U/mL Protector RNAse inhibitor) and filtered through a 30uM cell strainer. After three rounds of washing by centrifugation at 500g for 5 minutes, nuclei pellets were resuspended resuspension buffer (10mM Tris, 146mM NaCl, 1mM CaCl2, 21mM MgCl2, 2% BSA, 0.02% Tween-20) and filtered through a 20uM strainer. This nuclei suspension was loaded onto a MACSQuant Tyto HS chip and diluted with 1x PBS. Nuclei were sorted using gates set for isolation of GFP+ singlet nuclei. Sorted nuclei were subsequently used for preparation of four 10X Genomics Chromium libraries. Protocol can be found at dx.doi.org/10.17504/protocols.io.14egn6wryl5d/v1.

Library preparations were performed by the Northwestern University NUSeq Core Facility. Nuclei number and viability were first analyzed using Nexcelom Cellometer Auto2000 with AOPI fluorescent staining method. Sixteen thousand nuclei were loaded into the Chromium Controller (10X Genomics, PN-120223) on a Chromium Next GEM Chip G (10X Genomics, PN-1000120), and processed to generate single nucleus gel beads in the emulsion (GEM) according to the manufacturer’s protocol. The cDNA and library were generated using the Chromium Next GEM Single Cell 3’ Reagent Kits v3.1 (10X Genomics, PN-1000286) and Dual Index Kit TT Set A (10X Genomics, PN-1000215) according to the manufacturer’s manual with following modification: PCR cycle used for cDNA generation was 16 and the resulting PCR products was size-selected using 0.8X SPRI beads instead of 0.6X SPRI beads as stated in protocol. Quality control for constructed library was performed by Agilent Bioanalyzer High Sensitivity DNA kit (Agilent Technologies, 5067-4626) and Qubit DNA HS assay kit for qualitative and quantitative analysis, respectively.

The multiplexed libraries were pooled and sequenced on an Illumina Novaseq6000 sequencer with paired-end 50 kits using the following read length: 28 bp Read1 for cell barcode and UMI and 91 bp Read2 for transcript. Raw sequence reads were then demultiplexed and transcript reads were aligned to mm10 genome using CellRanger (v7.0.1 – RRID:SCR_017344 – https://www.10xgenomics.com/support/software/cell-ranger/latest) with --include-introns function. snRNA datasets can be found at GEO accession GSE271781.

Data curation and analysis

Outputs from CellRanger were read into Seurat version 5.0.1 using the Read10X command for each sample. Numbers of UMIs, features and ribosomal reads and mitochondrial reads were plotted for each dataset and used to determine cutoffs for quality control pre-filtering of each sample. The following QC filtering commands were used in R: control1 <- subset(control1, subset = nFeature_RNA > 1200 & nFeature_RNA < 7800 & percent.mt < 0.5 & nCount_RNA < 29000 & percent.ribo < 0.5); lrrk1 <- subset(lrrk1, subset = nFeature_RNA > 1200 & nFeature_RNA < 6500 & percent.mt < 0.5 & nCount_RNA < 26000 & percent.ribo < 0.5); control2 <- subset(control2, subset = nFeature_RNA > 800 & nFeature_RNA < 5000 & percent.mt < 0.5 & nCount_RNA < 15000 & percent.ribo < 0.5); lrrk2 <- subset(lrrk2, subset = nFeature_RNA > 1000 & nFeature_RNA < 7500 & percent.mt < 0.5 & nCount_RNA < 20000 & percent.ribo < 0.5) (post-filtering violin plots shown in Fig S1A). The male and female datasets were then normalized and integrated using the SCTransform V2 method¹⁰² in Seurat V5¹⁰³. Cells were initially clustered using 35 principal components at a resolution of 0.8. Plotting clusters on UMAP showed one small cluster exceedingly distant from all other clusters, which was suspected to represent EW nucleus cells and was subsequently removed. Afterwards, the integration and clustering pipeline was repeated without these cells to remove their impact on PCA. Final clustering was achieved using the FindNeighbors() command with the following parameters: dims = 1:32, reduction = “cca”, n.trees = 500, k.param = 40, and with the FindClusters() command with the following parameters: resolution = 0.8, algorithm = 1, group.singletons = TRUE, graph.name = “SCT_snn”. UMAP reduction was calculated using the RunUMAP command with the following parameters: reduction = “cca”, dims = 1:32, reduction.name = “umap.cca”, n.epochs = 500, min.dist = .2, n.neighbors = 1000. In total, the integration resulted in a final dataset of 28532 nuclei, with a median UMI count of 7750.5 and median of 3056 features. All clusters were represented in both male and female samples (Fig S1B). Sex of original mice was inferred for each cell by calculating the expression ratio of genes Uty and Eif relative to Tsix and Xist for each cell, revealing roughly equal proportions of male and female cells in the final dataset (Fig S1C). All clusters were represented in each individual sample (Fig S1D), and with no overt differences in proportions across conditions (Fig S1E, Fig 8C). Clusters 21 and 16 were removed from downstream analyses due to high expression of glial markers Atp1a2 and Mbp, respectively (Fig S1F). All differential expression calculations were also performed in Seurat using the FindMarkers() command; individual parameters for differential expression at each step of our analysis are provided in the Awatramani Lab Github page (github.com/AwatramaniLab/Gaertner_Oram_etal).

ShinyCell web browser was created using the SCT assay of the Seurat Object described above¹⁰⁴. The code can be found at the Awatramani Lab Github page (github.com/AwatramaniLab/Gaertner_Oram_etal).

Mapping data query sets were done using TranserData() seurat function. The dataset and clusters from Azcorra & Gaertner et al. (Dataset can be found in the Gene Expression Omnibus GSE222558) ²⁵ was used as the reference, and the query dataset was the one described in Figure 1. Any cells clustered with a prediction.max.score <.5 were removed from the data set (2.4% of cells removed). The exact workflow can be found at the Awatramani Lab Github page (github.com/AwatramaniLab/Gaertner_Oram_etal).

Sankey plots were made using the Network3D package in R and the sankeyNetwork function. The ’source’ and ’target’ nodes are the corresponding cluster labels. ’Value’ is set as the number of cells clustered in both the source and target nodes. The exact workflow can be found at the Awatramani Lab Github page (github.com/AwatramaniLab/Gaertner_Oram_etal).

Cluster stability and homogeneity analyses

To quantify how homogenous each cluster is and uncover potential further subdivisions that may exist within our clusters, we applied two approaches highly similar to those used for our previously described dataset²⁵. First, we generated a measure of stability for each cluster through random downsampling and reclustering of the dataset using the same custom R scripts as previously described²⁵. Doing so provided cluster stability metrics (Fig S2A) representing the propensity for cells within a cluster to continue to co-cluster when data is removed.

To understand potential sources of lower cluster stability values (e.g. two clusters being grouped together as one, or additional small clusters being divided into two adjacent clusters), we utilized the ClusTree (v0.5.1 – RRID:SCR_016293 – https://github.com/lazappi/clustree) R package to plot clustering at ten incremental resolutions of 0.1 through 1.1 (Fig S2B).

Scanpy

Cells were clustered using the scanpy toolkit (Version #1.10.1 RRID: SCR_018139https://scanpy.readthedocs.io/en/stable/ on Python 3.10.13¹⁰⁵. Outputs from CellRanger were read into Scanpy using scanpy.read_10x_mtx(). Cells were filtered with the same parameters as described above in the Seurat workflow. Datasets were individually clustered using the Leiden clustering algorithm. The four datasets were then integrated using scanpy.tl.ingest.(). Each dataset was mapped onto Control 1. The exact workflow can be found at the Awatramani Lab Github page (github.com/AwatramaniLab/Gaertner_Oram_etal).

Pathway Enrichment Analyses

To evaluate pathways and gene sets enriched in subtypes or across conditions, we utilized the fgsea (v3.19; https://bioconductor.org/packages/release/bioc/html/fgsea.html; SCR_020938) and clusterProfiler (v4.8.3 – RRID:SCR_016884 – https://doi.org/10.18129/B9.bioc.clusterProfiler) R packages. Pathways for GSEA analysis (mouse canonical pathways gene sets) were obtained from MSigDB (v2023.2.Hs, 2023.2.Mm; https://www.gsea-msigdb.org/gsea/msigdb; SCR_016863). Custom R scripts were created for visualization of GSEA results and are available in the deposited code for this paper. fgsea was run using the following parameters: nperm=1000, stats=ranks, minSize = 25, maxSize = 500. For gene ontology, the biological processes gene sets were obtained using the org.Mm.eg.db R package (v3.19; https://bioconductor.org/packages/release/data/annotation/html/org.Mm.eg.db.html; SCR_002 774), and GO (part of clusterProfiler package) was run using default settings with BH adjustment of p-values and utilizing FDR-adjusted significant (p < 0.05) DEGs for the respective comparisons. Results were visualized using the enrichplot R package (v3.19; https://www.bioconductor.org/packages/release/bioc/html/enrichplot.html; SCR_006442).

scDRS Cell-level Scoring of GWAS Risk Associations

Gene expression data (‘corrected’ UMI counts) was extracted from the SCT. Expression matrices, cell metadata, and feature metadata were loaded using Pandas (v2.2.2; https://pandas.pydata.org/; SCR_018214) in Jupyter (Jupyter: v7.2.1; https://jupyter.org/install ; SCR_018315) notebook. The data were converted into an AnnData object using the anndata package (AnnData v0.10.8; https://anndata.readthedocs.io/en/latest/; SCR_018209) and stored in HDF5 format. Normalization and scaling were performed using the Scanpy package (v1.8.1) (PMID: 29409532). Total counts for each cell were normalized to 10,000 reads using the sc.pp.normalize_total function and log- transformed using sc.pp.log1p. The sc.pp.highly_variable_genes function (parameters: min_mean=0.0125, max_mean=3, and min_disp=0.5) was used to identify highly variable genes. Finally, the anndata object was scaled using the sc.pp.scale function, with a max_value parameter of 10.

GWAS summary statistics for Parkinson’s Disease risk were obtained from EBI GWAS Catalog Study ID: GCST009325 (PMID: 31701892) and processed using the MAGMA tool on FUMA’s online platform ( v1.5.2; https://fuma.ctglab.nl/snp2gene; SCR_017521) to generate a scored list (z-scores) of significant genes. Default settings were used with a 10×10 gene window. Scored gene list was processed using the munge-gs function of the scDRS package (v1.0.2; https://martinjzhang.github.io/scDRS/) to create .gs gene sets (PMID: 34431100) and loaded into the Jupyter Notebook using the scdrs.util.load_gs function. To align with the anndata object data, the gene set’s human genes were mapped to mouse orthologs and intersected with genes present in the dataset. Preprocessing was completed using the scdrs.preprocess function, binning genes by mean and variance with parameters n_mean_bin=20 and n_var_bin=20. Scoring for Parkinson’s disease was performed using the scdrs.score_cell function (parameters: ctrl_match_key=mean_var, n_ctrl=1000, weight_opt=vs, return_ctrl_raw_score=False, and return_ctrl_norm_score=True).

Bootstrapped 95% confidence intervals for mean score within clusters or cluster families were calculated in R using 10,000 permutations for each group and plotted as mean scDRS score per group with CI represented as error bars.

Parkinson’s Disease MAGMA Analysis Ranked Genes dataset can be found at 10.5281/zenodo.13076447

SynGO Analyses

SynGO analysis protocol is described in detail at dx.doi.org/10.17504/protocols.io.x54v92r94l3e/v1.

MERFISH

Animals

Adult C57Bl/6 male mice (RRID:IMSR_CRL:027) aged 2-3 months were used for the MERFISH experiments. Animals were housed at the Montreal Neurological Institute (MNI) on a 12h-12h light-dark cycle (light cycle 7:00 to 19:00) with ad libitum access to food and water. Experiments were approved by the MNI Animal Care Committee and conducted according to guidelines and regulations from the Canadian Council on Animal Care.

Gene selection and panel assembly for MERFISH

Hybridization probes were generated by Vizgen using their custom gene portal. DA neuron specific genes were selected based on differential gene expression analysis from the 20 clusters identified in the single-nuclei sequencing data. Also included were: 1) general markers of glutamatergic, GABAergic, cholinergic, serotoninergic neuron populations, 2) markers of non-neuronal cell types, 3) neuropeptide precursor genes and neuropeptide receptors and 4) general ion channels. The list of all probes on the panel can be found in Supplemental Table 1.

Tissue preparation and tissue imaging for MERFISH

Mice were perfused with 1x PBS followed by 4% paraformaldehyde (PFA). Their brains were extracted and post-fixed in PFA for 16 hours, embedded in O.C.T. (Tissue-Tek O.C.T.; 25608-930, VWR), and stored at −80 °C until sectioning. Frozen brains were sectioned at −18 °C on a cryostat (Leica CM3050S). Within 3 animals we collected 10um-thick sequential coronal brain sections every 200μm on specialized glass slides (approximately 5-6 sections per mouse). Seven sections were chosen with the best morphology across the midbrain for processing. Sample preparation was completed using Vizgen gene imaging kits provided according to manufacturer. All steps were performed with RNAse free conditions. Sections were first washed with 1x PBS three times before incubation in 70% Ethanol overnight at 4°C to permeabilize cell membranes. The following day each section was washed with sample wash buffer followed by an incubation in formamide wash buffer at 37 °C for 30 minutes. The gene panel mix containing hybridization probes was then placed on top of the tissue and placed in a humidified 37 °C chamber for 36-48h. Afterwards the tissue was incubated with formamide buffer for 30m at 47 °C twice before washing with sample wash buffer. The hybridized tissue was then embedded in a 0.05% ammonium persulfate gel and cleared overnight in a solution of Proteinase K (NEB Cat# P8107S) at 37 °C. The following morning the sample was incubated in a solution containing DAPI and PolyT stain to aid in cell segmentation. Once processed, brain sections were prepped for imaging using the commercial Merscope system developed by Vizgen. 500-gene panel cartridges were thawed in a 37°C water bath for 1h before the start of imaging. Mouse RNAse inhibitor (NEB Cat# M0314) was added to imaging activation buffer and added to the cartridge fluidic system before loading into the Merscope to prime the fluidics chamber. Slides were removed from wash buffer and placed in the imaging capsule, connected to the fluidics chamber. Sections were imaged according to the manufacturer (Vizgen). The exact protocol for tissue preparation and tissue imaging for MERFISH can be found at dx.doi.org/10.17504/protocols.io.14egn6opyl5d/v1.

Preprocessing and quality control of MERFISH data

Cell by gene matrices were processed in R using the Seurat V5 package¹⁰³. To exclude cells that underwent poor cell segmentation we removed cells that had >500μm³ and <4000μm³ volumes. To remove cells that show low transcript identity or falsely identified cells we removed any cell that had <40 total RNA transcripts detected. The remaining cells then underwent SCT normalization¹⁰² before running RPCA integration to integrate datasets from individual sections. Shared nearest neighbour clustering was completed using the first 37 principle components using a resolution of k=0.7. Genetic markers of each cluster were determined using a roc test and was used to annotate our dataset. High-level annotations for glia and neuron types were confirmed using the SCType package in R (v3.0; https://sctype.app/). The putative DA cluster was subsetted, re-normalized and subclustered using the first 19 principle components and a resolution of k=0.65. Gene markers of each subcluster was determined using a ROC test and the top three genes from each cluster were identified based on the high positive log-fold change with a low p-values (adjusted p-value<10^-50). Code can be found on the Awatramani Lab Github page (github.com/AwatramaniLab/Gaertner_Oram_etal).

Integration of MERFISH data with scRNA-seq data

To integrate MERFISH data with snRNA-seq data we utilized the weighted-nearest neighbour algorithm defined in Stuart, Butler, et al ¹⁰⁶. Briefly, anchor points were established between the two datasets using the FindTransferAnchors and MapQuery functions in Seurat. This allowed for a projection of the MERFISH data into the scRNA-seq space and a transfer of cluster identity between the two. Importantly, a label transfer prediction score was assigned to each MERFISH cell and cells with scores greater than 0.5 were classified as high confidence cells. Whole genome expression was imputed onto MERFISH cells using snRNA-seq data using the TransferData function in Seurat using the snRNA-seq normalized data as a reference. Cluster identifications were compared between transfer data and MERFISH subclustering using a Sankeyplot in R. Briefly, the number of cells that transferred between two pairs of clusters were calculated and expressed as a proportion within each MERFISH cluster. Only flows >10% were included in the visualization. Spatial visualizations of DA subtypes were generated using the ImageDimPlot function in Seurat and overlaid on the Allen Brain Atlas mouse anatomical template or the Paxinos mouse brain coronal axis using affinity designer. Gene expression values were obtained using the ImageFeaturePlot function in Seurat. Code can be found on the Awatramani Lab Github page (github.com/AwatramaniLab/Gaertner_Oram_etal).

MERFISH anatomical localization

We utilized the Vizualizer software (v2.2, Vizgen Inc. https://portal.vizgen.com/resources/software) to delineate mouse midbrain regions into 11 distinct regions: the caudal linear nucleus (CLi), interfascicular nucleus (IF), interpeduncular nucleus (IPN), periaqueductal grey region (PAG), posterior hypothalamus (PH), rostral linear nucleus (RLi), retrorubral field (RR), substantia nigra pars compacta (SNc), Substantia nigra pars reticulata (SNr), Superior Colliculus (SC) and the ventral tegmental area (VTA). Brain regions were manually sectioned based on the intensity of DAPI imaging as well as the presence of Th transcripts compared to the Allen Brain Atlas. Cells within each ROI were counted to generate the localization of each family and cluster by counts and proportions.

MERFISH data can be found at DOI #: 10.5281/zenodo.12636328.

Code Availability Statement

Code is available at Awatramani Lab Github page (github.com/AwatramaniLab/Gaertner_Oram_etal) and DOI # 10.5281/zenodo.13820931.

Acknowledgements

ZG and CO performed snRNAseq and MERFISH experiments and data analysis. CB, ES, AS and CC performed data analysis and shinycell app construction. LP and DD provided scientific input on DA classification and LRRK2 alterations. ZG, CO, LP, JFP, RA wrote the manuscript. We would like to acknowledge Dr. Kevin Petrecca, Dr. Phuong Le and Michael Luo for generously providing access to Merscope and technical advice. RA, DD, AS and LP were funded by Aligning Science Across Parkinson’s (ASAP-020600). For the purpose of open access, the authors have applied a CC-BY 4.0 public copyright license to this manuscript. RA was also funded by 1R01NS119690-01, P50 DA044121-01A1. LP was funded by R01 NS097901. ZG was funded by NINDS 1F31NS115524-01A1. JFP was funded by CIHR, HBHL, Parkinson Canada. CB was funded by Parkinson Canada and FRQS.

Supplementary Figure Legends

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Figure S8.