Single-nucleus transcriptome analysis reveals cell type-specific molecular signatures across reward circuitry in the human brain

Single cell/nucleus technologies are powerful tools to study cell type-specific expression in the human brain, but most large-scale efforts have focused on characterizing cortical brain regions and their constituent cell types. However, additional brain regions - particularly those embedded in basal ganglia and limbic circuits - play important roles in neuropsychiatric disorders and addiction, suggesting a critical need to better understand their molecular characteristics. We therefore created a single-nucleus RNA-sequencing (snRNA-seq) resource across five human brain regions (hippocampus, HPC; dorsolateral prefrontal cortex, DLPFC; subgenual anterior cingulate cortex, sACC; nucleus accumbens, NAc; and amygdala, AMY), with emphasis on the NAc and AMY, given their involvement in reward signaling and emotional processing. We identified distinct and potentially novel neuronal subpopulations, which we validated by smFISH for various subclasses of NAc interneurons and medium spiny neurons (MSNs). We additionally benchmarked these datasets against published datasets for corresponding regions in rodent models to define cross-species convergence and divergence across analogous cell subclasses. We characterized the transcriptomic architecture of regionally-defined neuronal subpopulations, which revealed strong patterns of similarities in specific neuronal subclasses across the five profiled regions. Finally, we measured genetic associations between risk for psychiatric disease and substance use behaviors with each of the regionally-defined cell types. This analysis further supported NAc and AMY involvement in risk for psychiatric illness by implicating specific neuronal subpopulations, and highlighted potential involvement of an MSN population associated with stress signaling in genetic risk for substance use.


Introduction
Recent advances in single-cell and single-nucleus RNA-sequencing (scRNAseq/snRNA-seq) technologies have facilitated the molecular characterization of diverse cell types in the postmortem human brain during development (Darmanis et al., 2015;Li et al., 2018a;Zhong et al., 2018Zhong et al., , 2020, and have been used to assess cell type-specific gene expression differences in the context of several brain disorders, including Alzheimer's disease, autism spectrum disorder, multiple sclerosis, and major depressive disorder (Mathys et al., 2019;Nagy et al., 2020;Schirmer et al., 2019;Velmeshev et al., 2019). Identification of cell type-specific gene expression signatures has contributed to the understanding of the relationship between molecular identity and cell function as it relates to brain health, neurological disease, and genetic risk for neuropsychiatric disorders, such as schizophrenia (Skene et al., 2018).
While substantial advancements have been made in understanding cell type heterogeneity within regions and across the human brain, the majority of snRNA-seq reports are limited to a small number of brain areas. These primarily include the hippocampus (HPC) (Franjic et al., 2020;Habib et al., 2017) and several heavily studied sub-regions of the cortex (Lake et al., 2016), including the dorsolateral prefrontal cortex (DLPFC) (Li et al., 2018a;Nagy et al., 2020), medial temporal cortex (Darmanis et al., 2015;Hodge et al., 2019), entorhinal cortex (Grubman et al., 2019), and anterior cingulate cortex . Molecular profiling of less studied cortical subregions including the subgenual anterior cingulate cortex (sACC), as well as striatal and limbic brain regions, including the nucleus accumbens (NAc) and the amygdala (AMY), is lacking in the human brain. The sACC, NAc, and AMY are interconnected within well-established circuit loops that mediate important behavioral and neurobiological functions, including signaling for reward and motivation as well as processing emotional valence, particularly for fearful and threatening stimuli (Haber and Knutson, 2010;Janak and Tye, 2015;Russo and Nestler, 2013). Importantly, the cellular composition of individual neuronal subtypes in these regions substantially differs from previously well-profiled cortical and hippocampal regions (Saunders et al., 2018;Zeisel et al., 2018).
For example, the NAc contains dopaminoceptive populations of GABAergic medium spiny neurons (MSNs) -the principal projecting cell type comprising up to 95% of neurons in rodent -that harbor unique physiological and cellular properties (Gerfen et al., 1990;Kawaguchi, 1997;Kronman et al., 2019;Russo and Nestler, 2013). Early functional characterization of MSNs revealed two distinct classes of MSNs based on expression of D1 versus the D2 dopamine receptors (D1-MSNs and D2-MSNs, respectively) (Lobo, 2009;Lobo et al., 2006). However, recent sc/sn-RNAseq studies in the rodent striatum, and in the NAc specifically, revealed more complex transcriptional diversity within broader D1 and D2-MSN subclasses than was previously appreciated (Gokce et al., 2016;Saunders et al., 2018;Stanley et al., 2020;Zeisel et al., 2018). Moreover, subpopulations of MSNs are differentially recruited in response to cocaine exposure, and mediate divergent functional effects on behavioral responses to drugs of abuse (Savell et al., 2020). Similarly, single-cell profiling studies in the rodent AMY identified specialized populations of Cck-expressing neurons that are preferentially activated by behavioral experience, including exposure to acute stress .
However, whether and to what extent this transcriptional diversity is conserved in these areas of the human NAc and amygdala has not yet been fully explored. Given evidence for the functional importance of unique cell types in these areas of the rodent brain, profiling these regions in human by snRNA-seq may identify analogous cell populations, which can then be analyzed in the context of neurobiological dysfunction in human brain disorders.
Further, while unique transcriptomic profiles have been identified for specialized cell types that are localized to specific brain regions (e.g. granule cells of the dentate gyrus; MSNs of the striatum), it remains unclear to what extent broad populations of "common" cell types, such as glutamatergic excitatory neurons, are transcriptionally similar both within and across different brain regions. Given that many snRNA-seq studies have used different nuclei isolation protocols (sucrose vs iodixanol gradients), cell enrichment techniques (NeuN or DAPI sorting), library/sequencing technologies (10x Genomics Chromium, SMARTseq, Fluidigm, DroNc-seq), and data analysis workflows (dimensionality reduction techniques, handling batch effects, cluster annotation), comparing cell type-specific signatures across studies and between brain regions has been computationally and methodologically challenging. While efforts to generate a human cell type atlas at the single-cell level are underway (Han et al., 2020), the landscape of specialized molecular cell types across the human brain remains largely unexplored.
Here we defined the molecular taxonomy of distinct cell types in subcortical regions (NAc and AMY), which act as key nodes within circuits that mediate critical brain and behavioral functions including reward signaling and emotional processing. We also validated molecular profiles for previously identified cell types in the HPC and DLPFC, and identify similar cell types in the sACC, an additional cortical region central to limbic system function that has been implicated in affective disorders. Furthermore, we evaluate cross-species conservation of NAc and AMY cell types between human and rodent, specifically focusing on comparisons of MSN sub-populations identified as playing key roles in reward-processing and addiction. Finally, by integrating genetic studies for substance use and neuropsychiatric disorders, we show differential cell type expression of genes associated with schizophrenia, autism spectrum disorder, major depressive disorder, bipolar disorder, posttraumatic stress disorder, and features of addiction, highlighting the clinical relevance of understanding cell type-and regionspecific expression in the human brain.

Results
We profiled 5 brain regions (DLPFC, HPC, sACC, NAc, and AMY) across up to 5 neurotypical, adult male subjects of European Ancestry using 10x Genomics Chromium technology. To minimize potential batch effects, regions/donors were split across Chromium runs, for a total of 14 samples (sample/demographic information found in Table S1). Nuclear preparations were generated and purified by flow cytometry using DAPI staining (and NeuN selection for a subset of samples) to obtain nuclei from all cell types in a brain region. After sequencing, data processing and QC (Methods ; Table S2), we analyzed a total of 42,308 highquality nuclei across these 14 samples, which were then analyzed in respective region-specific, in addition to pan-brain, analyses.
We next evaluated the conservation of NAc cell types across species by comparing our transcriptional profiles with those generated in a previous study of acute exposure to cocaine, which analyzed a total of 15,631 rat NAc nuclei (Savell et al., 2020). Correlation analyses between our NAc subpopulations with those derived from rat NAc revealed that glial populations, including astrocytes, microglia, oligodendrocytes, and oligodendrocyte progenitor cells, were highly conserved ( Figure 1F). Inhibitory interneuron populations were also wellcorrelated across species as rat Sst-expressing and likely-Pvalb-expressing clusters overlapped with human Inhib.1 and Inhib.3 clusters, respectively. We also observed substantial correlation between rat and human D1 and D2-MSNs, especially between rat Drd1-expressing MSNs and human D1.4 MSNs. Further, beyond the expected overlap of rat Drd2-expressing MSNs in the human D2.2 MSN subcluster, we additionally saw positive correlations across D1 and D2 MSN subtypes, such that rat Drd2-expressing MSNs also showed enrichment in our human D1.4 MSN subcluster. This result is not likely fully explained by co-expression of DRD1 and DRD2 in the same nucleus because, while we did find that ~9% of all MSNs expressed some degree of both DRD1 and DRD2, these dual-expressing nuclei were mainly restricted to the D1.3 subcluster, as opposed to D1.4 ( Figure 1B). Additionally, many of the top markers for either the D1.4 or D2.2 subcluster were highly expressed in both MSN clusters (Figure S4), suggesting that the majority of canonically dichotomous D1 or D2 MSNs may be more molecularly similar than previously believed. We did not observe strong enrichment for rat Drd3-and Grm8expressing MSNs in any of our human MSN subclusters. Likewise, D1.1 and D1.2 human MSNs did not appear to be well represented in rat MSN subtypes (see Discussion). Taken together, while these data suggest strong overall conservation between rat and human NAc cell types, there appear to be transcriptional features that are unique among specialized subpopulations of rodent and human MSNs.

Atlas of molecularly-defined cell types in amygdala
The amygdala, a medial structure of the temporal lobe is noted for its role in processing emotional valence, particularly for both fear and reward (Janak and Tye, 2015;Wassum and Izquierdo, 2015). Dysfunction in amygdalar signaling is implicated in major depressive disorder, bipolar disorder and posttraumatic stress disorder (PTSD) (Fenster et al., 2018;Garrett and Chang, 2008;Murray et al., 2011). The human amygdala can be subdivided into a number of distinct regions based on histology, immunohistochemical classifications, connectivity, and neural activation patterns as revealed by functional magnetic resonance imaging (fMRI) of the brain (Barger et al., 2012;Schumann and Amaral, 2005;Sorvari et al., 1995;Tyszka and Pauli, 2016;Zhang et al., 2018). Studies in the rodent and non-human primate amygdala have identified different cell compositions across the amygdala, which likely correspond to differential patterns of synaptic connections between cell types across amygdalar subregions, and with extra-amygdalar brain regions (Chareyron et al., 2011). Hence, it is likely that various cell types with unique molecular signatures also exist within the human amygdala, which can be surveyed by snRNA-seq. We therefore analyzed 6,582 nuclei from the greater amygdala of two adult neurotypical donors to create a molecular taxonomy of cell types in this brain region. We identified 12 clusters that corresponded to four broad glial cell types and eight neuronal subclusters (Figure 2A; Figure S6). Glial cell populations were present at similar proportions between the donors (52% and 53% Oligo; 15% and 11% Astro; 13% and 10% Micro; 10% and 9% OPC), but we observed a varied distribution of neuronal subclusters between donors (see Discussion; Table S3). Despite this, after correcting for donor batch effects, we identified hundreds of genes enriched in each broad glial and neuronal subcluster at FDR < 1e-6 (top markers shown in Table S4).
Within the eight neuronal subclusters expressing the pan neuronal marker gene SNAP25, three clusters were enriched for excitatory neuronal markers (SLC17A7, SLC17A6) and five clusters were enriched for inhibitory GABAergic markers (GAD1, GAD2; Figure S6).
The GABAergic subcluster Inhib.3 was specific for expression of CALB1 (calbindin), on the other hand, and both Inhib.3 and Inhib.4 were enriched for CRH (corticotropin release hormone/factor). CRH is a key regulator of the hypothalamic-pituitary-adrenal (HPA) axis, which is critical for both the acute stress response and adaptation to chronic stress. Finally, NPFFR2 and TLL1, additional genes associated with HPA axis regulation, were selectively expressed in Inhib.5 (Lin et al., 2016;Tamura et al., 2005).
We then compared our subcluster-level transcriptomic profiles to those of a previously published single-cell dataset derived from the mouse medial amygdala (MeA)  to evaluate conservation of amygdalar cell types between humans and rodents ( Figure 2C). Across all shared homologous genes, we observed substantial correlation between several mouse and human amygdala cell types. For example, our human glutamatergic subcluster Excit.1 (SLC17A6+, SLC17A7+) most closely correlated with the mouse MeA glutamatergic subcluster 'N.11' (Pearson correlation: r = 0.455). Indeed the marker genes that were most highly conserved between these subclusters included SLC30A3, NPTX1, and NRN1.
Another notable pair of clusters conserved between species was human inhibitory neuronal subcluster, Inhib.5, and mouse inhibitory subcluster MeA 'N.8' (r = 0.465). The top shared genes between these clusters included NPFFR2, GRM8, and FOXP2. Though we observed selective co-expression of NPFFR2 and TLL1 in human Inhib.5, we note absence of orthologous Tll1 expression in all mouse MeA neuronal subclusters ( Figure S6), including the corresponding cluster 'N.8', suggesting species differences in the molecular characteristics of neuronal subpopulations. Importantly, several neuronal subpopulations in the mouse and human datasets lacked strong correlation with each other (e.g. human Excit.2, mouse 'N.3' and N.7'), either suggesting possible divergence between species, or unique differences between the cell-type makeup of amygdalar subregions, such that all subpopulations may not be fully represented in our human amygdala sample compared to mouse MeA samples. Our crossspecies analysis demonstrates the conservation of neuronal subtypes between human amygdala and mouse MeA, but highlights potential differences in the cellular distribution and transcriptomic profiles across neuronal subtypes.

Convergent cell classes with unique molecular signatures across brain regions
We lastly analyzed our 12 homogenate samples (a total of 34,005 nuclei) together at the pan-brain level to identify broad patterns of transcriptional dynamics across the human brain, assessing intra-and inter-regional similarities. We applied the same clustering, annotation, and marker-defining strategies as in region-specific analyses, allowing us to compare the higher-resolution cell type/subpopulation annotation, within each region, to clustering and annotations performed across brain regions. This analysis yielded 17 robust cell type clusters broadly classified into excitatory (Excit) and inhibitory (Inhib) neurons, oligodendrocytes (Oligo), oligodendrocyte precursor cells (OPC), astrocytes (Astro), and microglia (Micro; Figure S7). These pan-brain clusters showed some apparent regional specificity ( Figure 3A; Table S5). For example, excitatory neuronal subclusters Excit.1, Excit.2, Excit.5, and Excit.7 were nearly entirely comprised of nuclei from both cortical regions (DLPFC and sACC), likely reflecting the similar laminar distribution of neuronal subpopulations in the cerebral cortex. Alternatively, Excit.6 was driven by HPC nuclei, whereas Excit.8 had nuclei from both HPC and AMY. Ultimately, these 17 clusters defined across our five regions, analyzed together, provide only a broad overview of similarities between nuclei in a region-agnostic approach, as opposed to better-resolved regionally-defined populations (see Discussion).
To complement the subcluster populations described in the previous sections for the NAc and AMY, we additionally defined the catalog of cell type clusters and cluster-specific genes within all other brain regions (sACC, DLPFC, and HPC), separately ( Figure S7), applying methods using spatial transcriptomics data (Maynard et al., 2020a) to provide layer-specific annotations to our many excitatory neuronal subpopulations for the DLPFC (see Methods). We further benchmarked our transcriptomic profiles against other published datasets that profiled similar regions in the postmortem human brain. Overall, our HPC subpopulations correlated well with the broad cell classes as reported in (Habib et al., 2017); Figure S8). We additionally observed strong overlap between our DLPFC to the reported PFC profiles from ; Figure S9), or similarly, sACC to the ACC set ( Figure S10). Interestingly, our sACC subpopulations did not correlate more strongly with the ACC subpopulation profiles than their corresponding PFC profiles from , and our DLPFC subclusters generally correlated only slightly more strongly to PFC than ACC subpopulations. This suggests that these cortical regions share a high degree of overlap in their nuclear transcriptomic profiles.
The strength of correlation to these benchmark datasets demonstrates the robustness and utility of our pipeline, and the presented data significantly expand the existing repository of postmortem brain snRNA-seq datasets.
We next compared all regionally-defined subcluster expression patterns to generate a global view of the transcriptomic architecture across the five brain regions. Each glial cell subpopulation (Oligos, Astros, OPCs, and Micros) had highly consistent gene expression patterns across all five brain regions, suggesting a lack of overall regional specificity ( Figure S11), in line with previous analyses of broad non-neuronal cell populations using DNA methylation data (Rizzardi et al., 2019). Within the neuronal set of region-specific annotations, totaling 47 neuronal subpopulations, most inhibitory or excitatory populations preferentially expressed genes that clustered these classes together across brain regions ( Figure 3B).
Notably, all neuronal subtypes in the NAc (interneuron or MSNs) correlated poorly or lessstrongly with either excitatory or inhibitory subtypes from all other regions, reflecting divergent nuclear RNA expression profiles and the unique functional specializations of the neuronal populations in this brain region. These comparisons also suggested a potential excitatory function in the NAc-specific MSN.D1.4 subcluster, as this profile correlated more strongly with excitatory subclusters across other brain regions than inhibitory subclusters. This particular D1 MSN subcluster most highly correlated with one D2 MSN subcluster (MSN.D2.2; Pearson correlation: r = 0.50), reflecting their sharing of top marker genes ( Figure S4), and strikingly these two NAc subpopulations negatively correlate with all other MSN or inhibitory interneuron populations within the NAc, further highlighting the specializations across NAc neuronal subpopulations. We also observed strong similarities between unique pairs of neuronal subpopulations across regions, such as between AMY and DLPFC ('In.3_amy' and 'In.3_dlpfc'; r = 0.74). Indeed, this AMY inhibitory subpopulation shares many top markers with its DLPFC counterpart (Table S4), such as GABRD and MME, the latter encoding for the neuropeptidecleaving protein, Neprilysin. Further, we observed analogous pairings between each of layerspecific excitatory DLPFC and spatially-undefined sACC subpopulations (e.g. 'Ex.L5_dlpfc' and 'Ex.4_sacc'; r = 0.86). This suggests that spatially-registered snRNA-seq information, as generated in the DLPFC (Maynard et al., 2020a), can be projected into regions where cell type architecture is expected to be similar. Such an approach is useful since spatial transcriptomic data generation has not yet reached the pace of sc/sn-RNAseq technologies, and remains unavailable for most brain regions.

Enrichment of region-specific cell subtypes in psychiatric disease and substance use
Genome-wide association studies (GWAS) have identified a plethora of genetic risk variants or loci (segregating variants in linkage disequilibrium, or LD) for common psychiatric disorders, including schizophrenia (SCZ: (Pardiñas et al., 2018; Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014)), autism spectrum disorder (ASD: (Grove et al., 2019)), bipolar disorder (BIP: (Stahl et al., 2019)), major depressive disorder (MDD: (Howard et al., 2019;Wray et al., 2018)), and posttraumatic stress disorder (PTSD: (Nievergelt et al., 2019)). Additionally, a large GWAS was recently performed with 1.2 million individuals to identify the genetic risk and correlations for alcohol and tobacco use (Liu et al., 2019). Approaches have been developed to identify the biological context or relevance of the hundreds of risk loci that are often identified for a given disorder or phenotype, such as LD score regression (Finucane et al., 2015), which assesses the heritability of complex phenotypes across input categories/genomic regions and their measured LD with single nucleotide polymorphism (SNP)level variants. Multi-marker Analysis of GenoMic Annotation (MAGMA) (de Leeuw et al., 2015) is an alternative approach that defines gene-level localization of GWAS risk, then integrates this with gene set observations, affording flexibility to assess a variety of marker lists, such as for our brain region-specific snRNA-seq subcluster profiles, in two separable analyses.
We used MAGMA to identify which cell subtypes in this study harbored aggregated genetic risk for psychiatric disorders, and found robust signals across many nuclear profiles in each of the five profiled brain regions. As expected, many DLPFC and HPC neuronal subtypes exhibited significant effect sizes for both SCZ and BIP GWAS, spanning both excitatory and inhibitory subpopulations (Figure S12), extending and strengthening previous findings in (Bryois et al., 2020;Skene et al., 2018). However significant association signals with BIP genetic risk only remain in the excitatory, layer-specific DLPFC subpopulations, after controlling with the strict Bonferroni multiple test correction across all MAGMA gene set tests (threshold p-value < 6.8e-5). Conversely, subpopulations of both excitatory and inhibitory classes associated with BIP risk at this threshold in HPC. This was additionally the case for sACC subclusters, which suggests potential regional differences in inhibitory subpopulations between the two cortical brain regions, in their manifestation of genetic risk for bipolar disorder. None of the subcluster profiles in these cortical or hippocampal regions retained significant signal for aggregated risk for GWAS in ASD, MDD, or PTSD, after Bonferroni correction, though there were some region subpopulations with less-stringent FDR-significant signals (controlling for false discovery rate < 0.05) across all tests for these disorders (Table S6).
Previously, it has been shown that broad mouse striatal neuronal populations (interneurons, Drd1, and Drd2-expressing medium spiny neurons, or MSNs) additionally associated with SCZ (Skene et al., 2018) and BIP (Bryois et al., 2020) genetic risk. We demonstrate that some of our refined subpopulations in the human NAc, including MSN.D1.4, and MSN.D2.1, and D2.2 exhibit strong associations to schizophrenia with variable effect sizes, after Bonferroni correction ( Figure 4A). MSN.D1.3 and MSN.D2.1 both additionally associated with BIP at this threshold, and interestingly, D1.3 was the strongest-associating subpopulation to BIP, whereas D2.1 had the strongest association to SCZ. Within the amygdala, we observed significant associations in each of our neuronal subpopulations to SCZ ( Figure 4B). Notably, the strongest signal and effect sizes were observed across each of Inhib.1 through Inhib.5, in comparison to AMY excitatory subpopulations, and Inhib.4 exhibited the strongest effect size across all regionally-defined subpopulations tested in any GWAS phenotype ( = 0.21). Thus, our analysis here with the NAc show not only complementary findings of Drd1 and Drd2expressing striatal MSN associations with schizophrenia and bipolar disorder, but we dissect these mouse association signals with more relevant human interneuron and MSN subpopulations, and further extend this analysis to human AMY snRNA-seq-defined subpopulations.
We further tested for alcohol and tobacco use GWAS (Liu et al., 2019) genetic risk associations across subcluster profiles from each of our brain regions, focusing on the subcortical regions centered in reward circuitry, the NAc and AMY, and their subcluster profiles described above. This highlighted various MSN and inhibitory subpopulations in the NAc associating with genetic risk for regular smoking behavior ('SmkInit') at FDR < 0.05 (Figure 4A), along with associations to multiple AMY excitatory and inhibitory neuronal subpopulations ( Figure 4B). Additionally, both Oligo and OPC clusters from both regions exhibited aggregated risk for regular smoking behavior. However after applying the above Bonferroni threshold across all regions and phenotypes tested, only the AMY 'OPC' population retained significant association with this behavior, whereas no NAc subpopulations met this significance threshold.
Conversely, only the NAc contained a subpopulation harboring Bonferroni-significant risk for any of the other phenotypes in the addiction GWAS. This was the GABRQ and TAC1-expressing subcluster MSN.D1.3, which associated with heaviness of drinking ('DrinkWk'). Collectively, these results provide complementary human findings for genetic risk associations to those previously described for psychiatric disease, further identifying subpopulations in the NAc and AMY harboring aggregated risk for substance use behaviors.

Discussion
In this study, we used snRNA-seq to profile five human brain regions across the ventral striatum (NAc), limbic system (AMY and HPC), and two cortical subregions (sACC, DLPFC).
While single-nucleus transcriptomic profiling in the postmortem human brain has rapidly accelerated, most efforts to date (Mathys et al., 2019;Nagy et al., 2020;Velmeshev et al., 2019) have focused on cortical regions and the hippocampus. This is the first study, to our knowledge, to systematically profile and compare across additional human brain areas selected for their function associated with risk for neuropsychiatric disorders and addiction. Specifically, we focused on the NAc and the AMY given their roles in emotional processing and reward signaling. While this study was performed only among neurotypical donors, the strong cell typespecific associations to genetic risk for these disorders provide important information about disease etiology. This link to genetic risk is important, given that differential gene expression identified in case-control studies of postmortem tissue are difficult to interpret and may represent consequences (and not causes) of these disorders (Collado-Torres et al., 2019;Jaffe et al., 2020). More generally, understanding the transcriptomic architecture and cell type composition across the normal human brain is crucial to understanding the etiology of disease, the molecular pathology observed in postmortem tissues, and novel potential disease targets.
Our study is a significant contribution as it demonstrates differential enrichment of disease risk in snRNA-seq-defined cell populations across multiple brain regions, including the NAc and AMY, which have not yet been profiled at the single-nucleus level in the human brain.
The NAc is a central hub for reward signaling, and altered function in circuits encompassing the NAc is implicated in a number of psychiatric disorders as well as drug addiction. Hence, we sought to define molecular profiles for NAc cell types, with a specific focus on functionally dichotomous subtypes of DRD1-and DRD2-expressing MSNs. Consistent with studies that used single-cell sequencing to profile the mouse striatum, including ~1000 striatal cells in each study (Gokce et al., 2016;Stanley et al., 2020), we identified several discrete subpopulations of D1 and D2-expressing MSNs in human NAc. Similar to Gokce et al., we also identified 2 discrete subpopulations of D2 MSNS, including a subtype expressing the serotonin receptor gene HTR7. However, distinct from Gokce et al. and Stanley et al., we identified 4 discrete subpopulations of D1-MSNs, which we validated by smFISH. Several reasons may explain why we identified different discrete D1 cell types, including differences in species (human vs. mouse), region (NAc-specific vs. striatum-wide), sample preparation (whole cells vs. nuclei), number of MSNs profiled (about 8x greater in our dataset) and single cell technology employed (10x Genomics Chromium vs. SMART-Seq v2). However, in agreement with these studies, we also observed co-expression of DRD1 and DRD2 in a small subset of MSNs. While these dual-expressing neurons did not emerge as their own cluster, they were largely found in the D1.3 subpopulation ( Figure 1B). Interestingly, this cluster showed the strongest enrichment of genes associated with psychiatric disease and addiction, indicating that this particular subpopulation might be especially vulnerable to dysfunction in these disorders. Indeed, among D1 subtypes, D1.3 MSNs show selective expression of CRHR2, a gene encoding corticotropin releasing hormone receptor 2, suggesting that they may be particularly susceptible to the effects of corticotropin-releasing hormone (CRH), which is released and mediates the physiological and behavioral response to stress, modulating several neurotransmitter systems, including dopamine release (Bonfiglio et al., 2011;Payer et al., 2017). Given that dysfunction of the CRH system has been associated with many psychiatric disorders, including depression, anxiety, and PTSD (Claes, 2004), understanding which cell types express CRH receptors may aid in more specific targeting of the stress axis for therapeutic developments.
Similar to Gokce et al., we also observed promiscuous expression of "typical" D1 and D2 neuropeptide marker genes (TAC1 and PENK, respectively) in both D1 and D2 MSN subpopulations, providing further evidence that these classic markers may not be as selectively expressed as previously understood. Future studies using spatial transcriptomic approaches will be important to clarify whether TAC1-expressing D1 and D2 MSN subpopulations show topographical organization in the NAc core vs. shell. Anatomical location may explain differences in TAC1 and PENK expression in specific MSN subpopulations, as it is well established that specific neuropeptides are expressed in a spatial gradient across the core and shell (Prensa et al., 2003;Salgado and Kaplitt, 2015;Stanley et al., 2020;Voorn et al., 1989).
To better interpret clinical implications of studies that focus on circuitry encompassing the NAc in animal models, further understanding of similarities and differences across species for cell types that contribute to NAc function are important. While many cell populations were conserved between rat and human NAc (Savell et al., 2020), we did observe differences in specific MSN subpopulations, which may indicate unique molecular features between analogous MSN subpopulations and/or the existence of divergent MSN subclasses, as exemplified by the lack of a corresponding human MSN subpopulations with the rat 'Grm8-MSN' subpopulation ( Figure 1E). However, given the small positive correlations measured with human D1.1 (r = 0.19) and D1.3 (r = 0.26) subtypes, it is possible that this Grm8-expressing population encompasses the species-equivalent of these less abundant D1 subtypes, since D1.1 expresses abundant GRM8, while D1.3 expresses little (Figure S4), even though the latter correlates more strongly with the rat 'Grm8-MSN'. We also were unable to identify a population of cholinergic interneurons. While cholinergic interneurons are thought to be more abundant in the human neostriatum compared to the rodent, where they only account for ~0.3% of neurons (Graveland et al., 1985;Rymar et al., 2004;Tepper and Bolam, 2004), we think it is likely that the low rate of sampling and this population's relative rarity accounts for this lack of observation.
In addition to profiling NAc cell types, we also generated a molecular taxonomy of human amygdala cell types. We identified eight distinct neuronal subpopulations as well as accompanying gene marker annotations, including NRN1 (neuritin) and NPTX1 (neuronal pentraxin 1) for the AMY Excit.1 subcluster. Neuritin is a neurotrophic factor which modulates neurite outgrowth and plasticity (Yao et al., 2018), whereas neuronal pentraxin 1 regulates neuron excitability via synapse density (Figueiro-Silva et al., 2015). Additionally, the highest SLC17A6 (VGLUT2)-expressing subcluster, Excit.2, specifically expressed high levels of VCAN (Versican) amongst other neuronal subpopulations, which has multiple isoforms exhibiting different mechanisms for synaptic regulation (Horii-Hayashi et al., 2008). Among the diverse set of inhibitory subpopulations in the AMY, the stress modulator CRH was specifically enriched in Inhib.3 and Inhib.4. Top markers in the AMY Inhib.5 subcluster included NPFFR2 (Neuropeptide FF Receptor 2) and TLL1 (Tolloid-Like 1), which are both associated with glucocorticoid signaling and the response to stress (Lin et al., 2016(Lin et al., , 2017Tamura et al., 2005). Interestingly, TLL1 harbors top intronic genetic variant associations for PTSD in multiple European American cohorts (Xie et al., 2013), suggesting that dysfunction of neurons in Inhib.5 subcluster might mediate genetic risk for PTSD. Using MAGMA, we assessed whether the top nuclear marker genes for AMY subclusters associated with PTSD genetic risk (Nievergelt et al., 2019), and found that Inhib.5 was indeed one of two AMY neuronal subclusters that reached significance (FDR < 0.05) across all gene set tests ( Figure 4B). This result did not maintain significance using a stricter Bonferroni threshold; however this GWAS is underpowered, and it is anticipated that with increased sample size in future iterations, stronger association signals for PTSD genetic risk may be revealed. Comparing human AMY subcluster profiles to data from the mouse medial amygdala (MeA; ), we found that Inhib.5 and its corresponding population in mouse (MeA 'N.8' subcluster, Figure 2C) were the most strongly correlated neuronal subpopulations. While Tll1 expression was notably absent in mouse MeA, Npffr2 and other top MeA 'N.8' marker genes were shared with Inhib.5 ( Figure S6). These insights highlight the importance of deriving reference snRNA-seq datasets across the human brain, as molecular gene markers may not be shared across species between analogous neuronal subpopulations.
We next used the snRNA-seq data from the five profiled regions to ask whether identified cell subpopulations harbored aggregate genetic risk for various neuropsychiatric disorders and/or features of substance use. We confirmed previous findings by identifying strong associations for neuronal subpopulations in the DLPFC and HPC with both schizophrenia (SCZ) and bipolar disorder (BIP) (Bryois et al., 2020;Skene et al., 2018), and significantly extended these findings by providing associations with specific sACC excitatory and inhibitory populations ( Figure S12). Additionally, we not only confirmed previously observed associations to broad striatal populations defined in the mouse, but used our human NAc subcluster profiles to further refine these findings by demonstrating that SCZ most strongly associated with subpopulation MSN.D2.1, whereas D1.3 had the strongest association for BIP ( Figure 4A). This suggests, for the first time, that individual populations of dopaminoceptive (DRD1/2) neurons in the human NAc may be differentially associated with SCZ and BIP. We also found that specific subpopulations of inhibitory interneurons in the human AMY were preferentially associated with SCZ, with AMY Inhib.4 exhibiting the strongest effect across the five profiled brain regions.
These observations highlight a potential role for these subcortical brain regions in mediating genetic risk for SCZ and BIP.
As both the NAc and AMY play critical roles in reward signaling, we also evaluated enrichment of genetic risk for addiction or substance use behaviors (Liu et al., 2019).
Intriguingly, the genetic risk for adopting regular smoking associated more broadly across most neuronal populations, whereas other phenotypes assessed in this addiction GWAS showed more preferential associations to certain subpopulations. This suggests that the risk for adopting addictive-like behaviors might affect these brain regions more broadly than specific features of addiction ( Figure 4A/B). With regard to the other features, the MSN.D1.3 subpopulation significantly associated with genetic risk for heaviness of drinking ('DrnkWk') and smoking ('CigDay'), after Bonferroni and FDR correction, respectively. As a top marker for this subpopulation was CRHR2, this might be a key population in understanding these features of addiction. Indeed, many rodent studies have implicated CRH receptors in alcohol consumption and alcohol dependence (Heilig and Koob, 2007;Yong et al., 2014). Finally, though no neuronal AMY subpopulations met our strict Bonferroni threshold for association, the potentially HPA axis-involved 'Inhib.5' population exhibited more FDR-significant associations to addiction phenotypes in this GWAS than any other neuronal subpopulation, suggesting that this NPFFR2/TLL1-expressing AMY subpopulation might be of interest in understanding the reward circuitry underlying substance use. From these analyses, we surveyed our diversity of neuronal subpopulations profiled in the NAc and AMY for their clinical relevance in psychiatric disease and addiction behaviors. Additionally, we have extended such analyses for these regions, which have formerly only been performed on cell-type profiles defined in murine models (Bryois et al., 2020;Skene et al., 2018) to their relevant human context, and with increased resolution of molecularly-defined subpopulations. Finally, we narrowed down on those subpopulations manifesting the greatest genetic risk, potentially highlighting some neuronal subclasses mediating certain substance use behaviors.
Future studies conducted in multiple brain regions warrant important considerations in study and analysis design for interpretability. Here, we performed: (1) a 'pan-brain' analysis with all of our 12 homogenate (non-NeuN-selected) samples, combined in a region-agnostic manner ( Figure 3A); and contrasted this against (2) performing the above cell type clustering and analysis within each of the regions, separately. There were few benefits of the pan-brain approach, as identified clusters showed little regional variation (e.g. glial and broad neuronal populations) or near-complete segregation by brain region (Table S5). We therefore recommend prioritizing region-specific analyses even when data from multiple brain regions are collected, and then integrating information from these analyses across brain regions, over performing combined-region analyses. Of course, the optimal approach may depend on the research question being asked.
While we identified and characterized a diversity of robust neuronal subpopulations with our analytical pipelines for this study, we recognize that our sample sizes will not capture all cell types or subpopulations--even for our most-sampled brain region, the NAc, such as cholinergic interneurons mentioned above. The most direct evidence for this is that there remains some bias in donor makeup of certain subpopulations (Table S3). However, despite steps to mitigate the impact of the small input for our sample processing protocol (see Methods), we expect some degree of sampling bias since cell type makeup is not expected to be homogeneous within a single region. For example, the NAc core or shell have different functional properties, and differ in regards to their afferent and efferent connections, and thus differences in cell composition across these two subregions is expected (Heimer et al., 1991;Li et al., 2018b;Zahm and Heimer, 1993). Integration of spatial transcriptomic technologies with snRNA-seq data in these regions (Maynard et al., 2020a) will help resolve expected heterogeneity across these adjacent subregions. To similar effect, we acknowledge that we may have missed some expected nonneuronal, or potentially even neuronal, subpopulations in our approaches to define these subcluster profiles (see Methods). For example, we expected some number of endothelial cells and pericytes in each of our regions, making up the vasculature of the brain; however, these likely clustered with the more abundant microglia, due to a shared cell lineage (Bailey et al., 2006). Ultimately, our goal was to better characterize neuronal subpopulations across these five brain regions, and with these sample sizes being still limited, we did not pursue any nested or further cluster refinement of the presented populations. With all of this in mind, we were able to validate the presence and relative ratios of subpopulations defined by our snRNA-seq analytical pipeline via smFISH.
Another caveat to these snRNA-seq data is the lack of gene expression information from the cytosolic compartment, such as the neuropil. This is an important caveat given that synaptic signaling is implicated in neuropsychiatric disorders, and gene products localized to the synapse are enriched for SCZ genetic risk (Skene et al., 2018). In addition, mRNA from some expected marker genes, e.g. PVALB, may be preferentially localized to the cytosol, as demonstrated with smFISH for the GAD1+ interneuron 'Inhib.3' population in the NAc (Figure S5). However, this seems to be cell population-specific, as PVALB was highly expressed in some DLPFC subpopulations (data not shown; see Data and Software Availability). These and observations by others thus emphasize that snRNA-seq will not capture the full transcriptomic profile of cell populations, including activation-induced or disease-associated molecular changes restricted to the cytosol (Thrupp et al., 2020). However, as we have previously demonstrated (Maynard et al., 2020a), snRNA-seq-defined cell populations can be registered to spatial transcriptomic data, which does retain such information, for further characterization of transcriptomic profiles.
In summary, we used snRNA-seq to profile five human brain regions. We defined transcriptomic profiles for 68 regionally-defined cell type subpopulations and characterized the architecture of molecular relationships across brain regions. We finally identified associations with genetic risk for neuropsychiatric disorders and addiction in unique neuronal subpopulations in the NAc and AMY. This study takes a large step towards mapping the single-nucleus transcriptomic atlas of the human brain, further demonstrating such a utility in understanding the diversity of cell populations and their roles in biology and disease.

Post-mortem human tissue
Post-mortem human brain tissue from five neurotypical donors of European ancestry from age 40 to 62 (Table S1) was obtained by autopsy from the Office of the Chief Medical

Examiner for the State of Maryland under State of Maryland Department of Health and Mental
Hygiene Protocol 12-24. Clinical characterization, diagnoses, and macro-and micro-scopic neuropathological examinations were performed on all samples using a standardized paradigm, and subjects with evidence of macro-or micro-scopic neuropathology were excluded. Details of tissue acquisition, handling, processing, dissection, clinical characterization, diagnoses, neuropathological examinations, RNA extraction and quality control measures have been described previously (Lipska et al., 2006). Dorsolateral prefrontal cortex (DLPFC) and hippocampus (HPC) tissue was microdissected using a hand-held dental drill as previously described (Collado-Torres et al., 2019). The subgenual Anterior Cingulate Cortex (sACC) was dissected under visual guidance from the medial aspect of the forebrain at the level of the rostrum of the corpus callosum. Dissections were performed ventral to the corpus callosum, and dorsal to the orbital frontal cortex (BA11). Medially it was bounded by the interhemispheric fissure, while laterally it was bounded by the corona radiata/centrum semiovale. For the amygdala, a block containing the structure was dissected under visual guidance at the level of its maximal size, taken from a 1 cm thick slab of one hemisphere, and sectioned in the coronal plane. The amygdala block was chosen by visual inspection at a level that contained the maximal number of subnuclei. Landmarks for selection of the amygdala block included presence of the internal and external segments of the globus pallidus, the anterior commissure, and optic tract. The block containing the nucleus accumbens was taken from a 1 cm thick slab of one hemisphere, and sectioned in the coronal plane. The NAc block was chosen at a level where the putamen and caudate are joined by the accumbens at the ventral aspect of the striatum, with clear striations separating the putamen from the caudate. Additional landmarks include the presence of the anterior aspect of the temporal lobe and the claustrum.
Nuclei were resuspended in the EZ lysis buffer, centrifuged again, and equilibrated to nuclei wash/resuspension buffer (1x PBS, 1% BSA, 0.2U/μL RNase Inhibitor). Nuclei were washed and centrifuged in this nuclei wash/resuspension buffer three times, before labeling with DAPI (10μg/mL). For 2 NAc samples from individual donors, nuclei were additionally labeled with Alex snRNAseq raw data processing We processed the sequencing data with the 10x Genomics' Cell Ranger pipeline, aligning to the human reference genome GRCh38, with a reconfigured GTF such that intronic alignments were additionally counted given the nuclear context, to generate UMI/featurebarcode matrices (https://support.10xgenomics.com/single-cell-geneexpression/software/pipelines/latest/advanced/references). Per the output metrics of Cell Ranger, each sample was sequenced to a median depth of 253.0M reads (IQR: 148.7M-274.9M). We started with raw feature-barcode matrices from this output for analysis with the Bioconductor suite of R packages for single-cell RNA-seq analysis (Amezquita et al., 2020) using Bioconductor (Huber et al., 2015) versions 3.10 and 3.11. For quality control (QC) and nuclei calling, we first used a Monte Carlo simulation-based approach to assess and exclude empty droplets or those with random ambient transcriptional noise, such as from debris (Griffiths et al., 2018;Lun et al., 2019). This was then followed by mitochondrial rate adaptive thresholding, which, though expected to be near-zero in this nuclear context, we applied a 3x median absolute deviation (MAD) threshold, to allow for flexibility in output/purity of nuclear enrichment by FACS using scater's isOutlier (Lun et al., 2016). This QC pipeline yielded 5,399 high-quality nuclei from the DLPFC, 10,444 nuclei from HPC, 6,632 nuclei from AMY, 7,047 nuclei from sACC, and 13,241 nuclei from NAc. Collectively, these exhibited a median unique molecular identifier (UMI) count of 8,747 (interquartile range, IQR: 5,895 UMIs) per nucleus, and a median detected gene count of 3,047 (IQR: 2,224-5,359) genes captured per nucleus. These feature-barcode gene counts were then rescaled across all nuclear libraries, using scater's librarySizeFactors (Lun et al., 2016). Finally, these rescaled counts were log2-transformed for identification of highly-variable genes (HVGs) with scran's modelGeneVar (Lun et al., 2016), taking all genes with a greater variance than the fitted trend.

Dimensionality reduction and clustering
Principal components analysis (PCA) was then performed on the HVGs to reduce the high dimensionality of nuclear transcriptomic data, both in region-specific analyses and panbrain analyses. The optimal principal component (PC) space was defined with iterative graphbased clustering to determine the d PCs where resulting n clusters stabilize, with the constraint that n clusters </= (d + 1) PCs (Lun et al., 2016), resulting in a chosen d between 45-96 PCs, in the region-specific analyses, or 204 PCs for the pan-brain analysis (all homogenate, DAPI-sorted samples). In this PCA-reduced space, graph-based clustering was performed to identify what we classified as preliminary clusters; specifically, k-nearest neighbors with k=20 neighbors and the Walktrap method from R package igraph (Csardi and Nepusz, 2006) for community detection. We used an increased k neighbors (from the default k=10) as a means to increase the connectivity in the kNN graph, as an alternate approach for handling potentially donor-driven preliminary clusters, instead of manually identifying batch-correlated PCs and removing these.
We then took all feature counts for these assignments and pseudo-bulked counts (Crowell et al., 2019;Kang et al., 2018;Lun and Marioni, 2017) across these preliminary nuclear clusters, rescaling for combined library size and log-transformed normalized counts, as above. With the pseudo-bulked count profiles, we then performed hierarchical clustering to identify preliminary cluster relationships, and finally merged with the cutreeDynamic function of R package dynamicTreeCut (Langfelder et al., 2016), or keeping split clusters at the preliminary resolution, if generally well-represented across donors, as this suggested biologically valid subpopulations (for example, neuronal subtypes) as opposed to more likely batch-driven preliminary clusters.
However, in some cases, cluster marker identification (see below) suggested sample bias in true, biological subpopulations (see Discussion). The final clusters merged at the appropriate tree height were then annotated for broad cell type identity with well-established cell type markers (Mathys et al., 2019), and with a numeric suffix where multiple broad cell class populations were defined ('Excit.1','Excit.2',etc.). We also used Bioconductor package scater's (McCarthy et al., 2017) implementation of non-linear dimensionality reduction techniques, t-SNE (van der Maaten and Hinton, 2008) and UMAP (McInnes et al., 2018), with default parameters and within the aforementioned optimal PC space (or a reduced dimensional space in the case of the NAc and AMY, for t-SNE coordinates which would better reflect the final subcluster distribution), simply for visualization of the high-dimensional structure in the data, which complemented the clustering results. Additionally, in each of our five within-region analyses, a small cluster driven by low transcript capture would remain even after hierarchical cluster-merging of preliminary clusters, but these were removed prior to downstream analyses and from the t-SNE display, resulting in a final n nuclei analyzed per region of: 5,231 from the DLPFC, 10,343 nuclei from HPC, 6,582 nuclei from AMY, 7,004 nuclei from sACC, and 13,148 nuclei from NAc (an average of 98.8% nuclei kept post-QC, above). These final numbers of nuclei analyzed per regionally-defined subcluster and by donor can be found in Table S3.

Cluster marker identification
For marker identification with our final clusters defined in each brain region or at the panbrain-level analysis, we utilized scran's findMarkers (Lun et al., 2016) function for two sets of statistics: 1) Pairwise t-tests, to identify differences between each cluster, or 2) Implementing the function findMarkers to perform a cluster-vs-all-other-nuclei t-test iteration In both cases, we included a donor/processing date covariate to model (in the design= parameter) on these expected and unwanted batch effects. The latter approach, 2), we consider a less-stringent marker test for enriched genes in a given cluster, but which would not necessarily differentiate between said cluster and all others. We used the results from both tests to interpret cell type identity beyond the broad classes (excitatory vs. inhibitory neuron), and to identify markers to probe via smFISH (below). The top 40 markers from each test result are provided for each regionally-defined subpopulation in Table S4 (regions separated by worksheet), where the '_pw' suffix corresponds to the pairwise tests (set 1), and '_1vAll' to the enriched expression test (set 2).
Importantly, 2) can be used to return a statistic, Cohen's D, or the standardized log-fold change, which we used to back-compute a single t-statistic for each cluster per gene, using: t = std.logFC * sqrt(n), where n = the total n nuclei (per region/dataset) * Back-computing a single t-statistic cannot be generated with the result of 1) due to pairwise testing.
These t-statistics can finally be used to compare such 'transcriptomic profiles' to those we computed for publicly-available datasets, using the provided cell type annotations (or across our 5 regions), and compute the Pearson correlation coefficient, as was done in the spatial registration approaches in spatialLIBD (Maynard et al., 2020a). To perform cross-species conservation analyses, we generated these t-statistics per gene per reported cell annotation, subsetting on shared homologous genes between our human data and rat or mouse, using the 'HomoloGene.ID' identifier provided by (http://www.informatics.jax.org/downloads/reports/HOM_AllOrganism.rpt), before computing the pairwise correlations. In the case of "many-to-many" scenarios, we took the highest-expressing paralog as the surrogate for each homologous pair, though these were small sets of genes in both rat and mouse cases.

GWAS association analyses with MAGMA
The latest version (v1.08) of Multi-marker Analysis of GenoMic Annotation (MAGMA; (de Leeuw et al., 2015) was used to test for genetic risk association of our 68 regionally-defined subpopulations with schizophrenia (SCZ: (Pardiñas et al., 2018; Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014)), autism spectrum disorder (ASD: (Grove et al., 2019)), bipolar disorder (BIP: (Stahl et al., 2019)), major depressive disorder (MDD: (Wray et al., 2018)), posttraumatic stress disorder (PTSD: (Nievergelt et al., 2019)); and for alcohol and tobacco use (Liu et al., 2019). For the marker gene sets, we used any genes defined as enriched per subpopulation (using marker test set 2, from above), at the Benjamini & Hochberg false discovery rate (FDR) < 1e-12 (Benjamini and Hochberg, 1995). SNPs were first annotated to genes, using window sizes from -10kb to +35kb of each gene, with the 1000 Genomes EUR reference panel, and gene-level analyses were performed, using provided summary statistics from each of the above listed GWAS (via https://www.med.unc.edu/pgc/download-results/) and the snp-wise=mean model, to test whether there was enrichment of genetic risk for disease/phenotype in each gene. Following this, we performed the default competitive gene set analysis with the 68 regionally-defined marker sets, testing for association of gene-level risk and whether genes were enriched/specific to each subpopulation. From the empirical p-value of the gene set analysis, we performed multiple test correction with both false-discovery rate (FDR) and the stricter Bonferroni procedure (threshold p < 6.8e-5) across all 680 (68 regionallydefined subpopulations and 10 GWAS phenotypes tested) tests. All genetic association test results are provided in Table S6.

RNAscope single molecule fluorescent in situ hybridization (smFISH)
Fresh frozen NAc from two independent donors was sectioned at 10μm and stored at - were fluorescently labeled with Opal Dyes (Perkin Elmer, Waltham, MA; 1:500) and stained with DAPI (4′,6-diamidino-2-phenylindole) to label the nucleus. Lambda stacks were acquired in zseries using a Zeiss LSM780 confocal microscope equipped with a 63x x 1.4NA objective, a GaAsP spectral detector, and 405, 488, 555, and 647 lasers as previously described (Maynard et al., 2020b). All lambda stacks were acquired with the same imaging settings and laser power intensities. For each subject, high magnification 63x images were randomly acquired in the NAc (n= 2 subjects, n=2 sections per subject, n=12 images per section). Following image acquisition, lambda stacks in z-series were linearly unmixed in Zen software (weighted; no autoscale) using reference emission spectral profiles previously created in Zen (Maynard et al., 2020b) and saved as Carl Zeiss Image ".czi" files. Images were segmented and quantitatively analyzed in MATLAB using dotdotdot software (Maynard et al., 2020b) and statistical analyses were performed in R v3.6.3.
For each of the five experiments, we combined ROI-level data from all respective images, and used data-driven cutoffs based on distributional overlap to determine binary expression levels (i.e. expressed or unexpressed) for each gene/channel. For each experiment, we calculated the Euclidean distance of the vector of expressed targeted genes in each ROI to the cell type-specific targeted designs, where distance = 0 meant that ROI matched a predicted snRNA-seq cell type cluster (shown as larger points in Figures E and S1

Data and Software Availability
All the code for processing and analyzing the data is available on the GitHub repository https://github.com/LieberInstitute/10xPilot_snRNAseq-human. The processed results files are hosted on Amazon S3 and the links are available on the README.md in the GitHub repository.
For each of the five brain regions, we created an interactive website with the data using iSEE (Rue-Albrecht et al., 2018) and deployed at the LIBD shinyapps.io account at URLs such as https://libd.shinyapps.io/tran2020_Amyg/ (tran2020_sACC, tran2020_DLPFC, tran2020_NAc,  (Savell et al., 2020) derived from rat NAc with data-driven clusters provided in their processed data.   Heatmap is colored by empirical -log10(p-value) for each association test. Displayed numbers are the effect size ( ) for significant associations (controlled for false discovery rate, FDR < 0.05), on a Z (standard normal) distribution. Bolded numbers are those that additionally satisfy a strict Bonferroni correction threshold of p < 6.8e-5.

Supplementary Material
Supplementary Tables   Table S1. Demographic/Sample Info   Table S2. Cell Ranger QC and Mapping Statistics Table S3. Regionally-defined Nuclear Subpopulations, Donor-stratified Table S4a-e. Top 40 Marker Lists (Both Pairwise and Cluster-vs-All-Others Results) Table S5. Pan-brain Cluster Annotations, Region-stratified Table S6. MAGMA association statistics across all regionally-defined subpopulations