Fresh snap-frozen human post-mortem tissue was provided by the Banner Sun Health Research Institute Brain and Body Donation Program of Sun City, Arizona, United States, where all human subjects donate after signed informed written consent [Western Institutional Review Committee (WIRB), Seattle, Washington] (Beach et al., 2015). Sample data used in the preparation of this article were obtained from the Arizona Study of Aging and Neurodegenerative Disorders (AZSAND) database.¹ Detailed information about the samples analyzed is provided in Supplementary Figure 1. Controls (n = 10) didn’t have any neurodegenerative disease and were selected to match the age and the gender of the PD cases. Since loss of dopaminergic markers occurs rapidly and is completed by 4 years post-diagnosis (Kordower et al., 2013), we selected two groups of PD cases: early post-diagnosis 2-4 years (n = 6) and late post-diagnosis 5-24 years (n = 6). ILBD cases (n = 7) were also included. Only samples with a post-mortem index < 5h were included. To get robust results, we used ≥six cases per condition (Schurch et al., 2016).

Cryosectioning was used to obtain 10 μm thick sections of the ventral midbrain, which were mounted on membrane glass slides (Zeiss, 415190-9041-000) that were kept at −20∘C during sectioning and subsequently stored at −80∘C. Prior to LCM, the tissue was dehydrated by a series of ethanol washes (75% 30 s, 95% 30 s, 100% 30 s). The Leica LMD7000 system was used to capture NM+ neurons from the SNc. 10 controls, 12 PD and 7 individuals with ILBD were analyzed. Cells were cut at ×40 magnification with a minimum laser power and collected into lysis buffer (Picelli et al., 2013). From each SNc, 4–12 pooled samples of 100–200 micro-dissected cells were collected. As a reference, 18 samples were also collected from the cortex of control brains to provide a non-DA source of RNA. Cortex neurons were visualized using a quick Histogene (Arcturus) staining.

Libraries for RNA-seq were generated using the Smartseq2 protocol (Picelli et al., 2013). cDNA libraries were tagmented using home-made Tn5 enzyme (Picelli et al., 2014) and Nextera dual indexes (Illumina). The cDNA quality was checked on high-sensitivity DNA chip (Agilent Bioanalyzer). Illumina HiSeq 2000 was used for sequencing, giving 43 base pair reads after de-multiplexing.

Following de-multiplexing, reads were aligned to the hg19 human genome using STAR v2.3.0 and filtered for uniquely mapped reads (Dobin et al., 2013). Gene expression was calculated as reads per kilobase gene model and million mappable reads (RPKMs) for each transcript in Ensembl release 69 using rpkmforgenes (Ramsköld et al., 2009). The low-quality libraries were filtered out based on following criterias: <10.3% uniquely mapping reads, >68% fraction mismatches, <10.5% exon mapping reads, >57% 3′mapping, < 26% of all genes detected, <100000 normalization reads. 270 samples passed the quality control.

Principle Component Analysis (PCA) was used to analyze how individual samples relate to each other on overall gene expression profiles and DESeq2 package (Love et al., 2014) to determine the differentially expressed genes between sets of samples.

Pathway enrichment analysis was performed using the Toppgene suite (Chen et al., 2009) on differentially expressed genes detected with DESeq2 that fulfilled the following criteria: padj < 0.05, mean RPKM value in control samples > 1 and 1.5 fold change.

Snap-frozen tissue blocks, received from the Banner Sun Health Research Institute Brain and Body donation program, were used to obtain cryo-sections of the ventral midbrain. Laser capture microdissection (LCM) was used to isolate mDA neurons from postmortem cases of sporadic PD, ILBD and age-matched controls (Figure 1A and Supplementary Figure 1). mDA neurons were identified by the presence of neuromelanin (NM) and dissected using LCM (Figure 1B). Several replicates of 100–200 individually isolated mDA neurons per individual (either PD, ILBD or age-matched controls) were collected. Cortical neurons from control brains (n = 18 pools), visualized by Histogene quick staining, were also collected by LCM in order to provide a reference source of mRNA from a brain region devoid of mDA neurons. Libraries for RNA sequencing were generated by the Smart-seq2 protocol (Picelli et al., 2013). Figure 1C summarizes cases used in the study and the number of analyzed samples per case, while Supplementary Figure 1 provides more information on all cases included in the study.

Libraries were sequenced (see Methods for details) and analyzed by various bioinformatic methods. The Spearman coefficient between biological replicates of neuromelanin-positive (NM+) cell samples averaged 0.78 and was not below 0.57 for any sample, revealing a very high level of overall correlation (Supplementary Figure 2A). Principal component analysis (PCA) based on all variable genes showed that NM+ samples from healthy controls clustered separately from cortex samples along the second principal component (PC2) (Figure 2A). Further comparison of cortex and NM+ samples revealed that 8845 genes were differentially expressed in the two populations (padj < 0.05, RPKM in controls > 1, fold change > 1.5; Figure 2B and Supplementary Table 1). Known dopaminergic markers such as TH, SLC6A3, EN1 and PITX3 were identified among top-ranked genes enriched in NM+ samples demonstrating that these samples provide a good representation of the transcriptome of human mDA neurons (Figures 2B–E). Moreover, the comparison demonstrated that pan-neuronal markers such as GAP43 and MAP2 were enriched in the NM+ cells (Figure 2C) while glial cell markers were identified in cortex samples where dissection did not discriminate between neurons and other cell types including astrocytes and oligodendrocytes (Supplementary Figure 2B). This difference was also apparent from gene ontology analysis (Supplementary Figure 2C).

Interestingly, samples visualized by PCA showed that ILBD and control samples clustered mostly together while a large proportion of PD samples clustered separately along one of the principle components (PC5; Figures 3A,B). Along other PCs, separation was not evident (Supplementary Figures 3A,B). Thus, the analysis detected gene expression profiles that differed in PD samples as compared to controls and ILBD samples. In addition, disease duration affected the transcriptional profile of NM+ cells from PD samples along PC2 as evident from different clustering when comparing samples derived from patients with a disease duration of 2-4 years (PDearly) with those with disease duration of 5 to 24 years (PDlate); (Figures 3C,D). We next considered how additional parameters affected the transcriptional profile of PD samples and found that samples did not cluster separately across any principal component based on either gender, post-mortem intervals (PMI; < 2.5 h versus > 2.5 h) or patient age (< 80 years versus > 80 years) (not shown); however, consistent with the observed different clustering of PDearly vs. PDlate samples, a trend for separation along PC2 was also seen when comparing Braak stages (I, II versus III, IV; Supplementary Figure 3E). Thus, disease duration appears to have an influence on transcriptional profiles in mDA neurons remaining in PD patients at the time of death.

Consistent with PCA demonstrating differences between control and PD samples, DESeq2 statistical analysis (Love et al., 2014) identified hundreds of differentially expressed genes. In PDearly samples, 220 genes were expressed at a higher level and 1029 genes were expressed at a lower level compared to healthy controls (Supplementary Table 2); in PDlate samples, 448 genes were more highly expressed and 565 genes were expressed at a lower level (Supplementary Table 3); in ILBD samples, 277 genes were more highly expressed and 397 genes were expressed at a lower level (Supplementary Table 4), (padj < 0.05, RPKM in controls > 1, fold change > 1.5; Figure 4A). Figure 4B displays a heat map visualizing selected top enriched genes in NM+ neurons from PDearly (orange) and PDlate (red) and Supplementary Figure 4B from ILBD. Top downregulated genes are shown in Supplementary Figure 5.

We focused on the top 100 up- or down-regulated genes (padj < 0.05, RPKM in controls > 1, fold change > 1.5) in PDearly, PDlate and ILBD samples and searched literature in order to identify genes that have been previously associated to PD. A striking observation in PDlate samples was the higher expression of several genes shown previously to confer neuroprotection, such as the pituitary adenylate cyclase activating polypeptide ADCYAP1 (Reglõdi et al., 2004) the calcium-binding protein S100A4 (Pedersen et al., 2004), the transcription factor NFIL3 (Tamai et al., 2014), the regulator of G-protein signaling RGS10 (Tansey and Goldberg, 2010). Neuroprotective gene such as PRDX5 (Plaisant et al., 2003) was also expressed at a higher level in PDearly cases suggesting that some neuroprotective mechanisms are quite active also during early stages of the disease. Genes previously shown to be more highly expressed in SNc from PD patients, including DDIT4 (Malagelada et al., 2006) and NPTX2 (Moran et al., 2008), were also found to be enriched in PDlate cases. Of interest is also that variants of several regulated genes in both PDearly and PDlate samples, including GSTT1 (Wang et al., 2014), CXCR4 (Bonham et al., 2018), BAG3 (Cao et al., 2017) and HLA-DRB1 (Hollenbach et al., 2019), were previously associated with PD.

Relatively few genes expressed at lower levels in PD samples were found to be associated with PD. However, LGALS3 (Galectin-3), a gene linked to neuroinflammation, was expressed at a lower level in PDearly samples and a BCL2 antagonist, was also expressed at lower levels in PDearly samples [see e.g., (Bernardini et al., 2019)]. BAK1, the gene encoding a presumptive target of the E3 ubiquitin ligase Parkin, was expressed at lower levels in PDlate samples. Notably, quite many genes behaved in an opposite manner being enriched in PDearly samples but expressed at lower levels in PDlate samples. This further underscores the different profiles of gene expression seen in samples derived from patients at either early or late stages of disease (Supplementary Tables 2, 3).

We next focused on known mDA neuron markers and genes that cause familial PD (padj < 0.05, RPKM in controls > 1). The expression of known mDA neuron markers RET, GFRA1, LMX1B, and FOXA2 were not significantly affected in either PD group. Expression of DDC and LMX1A were the only genes that were significantly increased in PDearly but not in PDlate samples (Figure 5A). SLC18A2, NR4A2, EN2, and OTX2 were expressed at lower level only in PDearly, while TH, SLC6A3, DRD2, EN1, and PITX3 were low only in PDlate samples (Figure 5A). SOX6, encoding a transcription factor specifically expressed in SNc but not in ventral tegmental area (VTA), was previously shown to be expressed at lower levels in NM+ cells in PD patients (Panman et al., 2014). Consistent with this result SOX6 was the only analyzed mDA neuron marker that was significantly expressed at lower level in both PDearly and PDlate cases (Figure 5A). Most genes causing familial forms of PD did not show differential expression in any PD cohort as compared to controls. UCHL1 was significantly enriched only in PDearly cases, ATP13A2 was enriched in PDearly but lower than controls in PDlate cases, while GIGYF2, EIF4G1, and PLA2G6 were lower in PDlate but not in PDearly cases. PINK1 and NOTCH2NL were the only familial PD-associated genes that was expressed at lower levels in both early and late PD cases (Figure 5B).

The diversity of mDA neurons within the ventral midbrain is substantial. In PD, mDA neurons of the ventral-lateral SNc are more vulnerable than mDA neurons localized at a dorsal-medial region of the SNc and within the VTA (Gibb and Lees, 1991). Although tissue samples were dissected from the substantia nigra it is likely that NM+ cells isolated by LCM in this study represents mDA neurons belonging to different subtypes and with variable vulnerability (Poulin et al., 2020). A recent study identified markers specific for human SNc and VTA mDA neurons, respectively (bioRxiv 334417; doi: https://doi.org/10.1101/334417). The study defined 33 markers that can distinguish between human VTA and SNc mDA neurons. We assessed if any of these genes were enriched in either PDearly, PDlate or incidental Lewy body disease. Interestingly, a clear trend was observed when comparing the different samples (Figure 5C). Accordingly, the majority of VTA markers were significantly enriched in PDlate samples but not in PDearly samples indicating that surviving cells isolated from PDlate samples are enriched for cells with VTA or “VTA-like” gene expression profiles. In contrast, although loss of neurons within the SNc has been reported in ILBD (40%), VTA markers were not enriched in ILBD samples (Iacono et al., 2015) (Figure 5C).

The overlap between differentially expressed genes in PD and ILBD samples was limited (Supplementary Figure 4A). Of note, among genes that are differentially expressed in both PD and ILBD samples, several are among those that previously have been associated with PD. Thus, GPNMB, PRDX5, GSTT1, HLA-DRB1, and DDIT4 were enriched in both ILBD and PD samples and linked to PD as described above. Moreover, LGALS3, a gene associated with neuroinflammation, was expressed at reduced levels in both PD and ILBD samples.

Gene ontology analysis was performed to identify pathways that are enriched among genes that are differentially expressed at higher levels in PD and ILBD samples (padj < 0.05, RPKM in controls > 1, fold change > 1.5; Figure 6 and Supplementary Table 5). Terms related to negative regulation of apoptosis appears to be upregulated in PDearly samples. However, it should be noted that almost all genes contributing to these GO-terms belong to the same family of genes with homology to the mitochondrial encoded MT-RNR2 (humanin) gene (Supplementary Figure 6). Humanin has been shown to have anti-apoptotic and neuroprotective effects and has been associated with Alzheimer’s disease (Hashimoto et al., 2001). Several GO-terms relating to regulation of apoptosis and programmed cell death were also identified when analyzing enriched genes in PDlate (Figure 6 and Supplementary Table 5). Another set of GO terms among genes enriched in late PDlate samples were dominated by terms related to the response to unfolded proteins and chaperone-mediated protein folding. Genes associated with the endoplasmic reticulum (ER) unfolded protein response pathways (UPR) are included among these genes, e.g., XBP1, CHAC1, and DDIT3 (CHOP), but also many other genes encoding proteins that regulate cytoplasmic response to unfolded proteins (Supplementary Figure 6 and Supplementary Table 5). Enriched pathways identified from analysis of downregulated genes demonstrated that there was a significant downregulation of synaptic signaling and neuronal projection organization in PDlate samples suggesting defects in synaptic transmission and axon/dendrite integrity, which are hallmarks of PD pathology (Supplementary Figure 7). Finally, in ILBD samples, downregulated genes were associated with GO terms relating to genes involved in cell junction. Strikingly, enriched genes were strongly linked to terms associated with oxidative phosphorylation (Figure 6 and Supplementary Table 5).