Explorative Combined Lipid and Transcriptomic Profiling of Substantia Nigra and Putamen in Parkinson’s Disease

Parkinson’s disease (PD) is characterized by the loss of dopaminergic neurons from the substantia nigra (SN) that project to the dorsal striatum (caudate-putamen). To better understand the molecular mechanisms underlying PD, we performed combined lipid profiling and RNA sequencing of SN and putamen samples from PD patients and age-matched controls. SN lipid analysis pointed to a neuroinflammatory component and included elevated levels of the endosomal lipid Bis (Monoacylglycero)Phosphate 42:8, while two of the three depleted putamen lipids were saturated sphingomyelin species. Remarkably, we observed gender-related differences in the SN and putamen lipid profiles. Transcriptome analysis revealed that the top-enriched pathways among the 354 differentially expressed genes (DEGs) in the SN were “protein folding” and “neurotransmitter transport”, and among the 261 DEGs from putamen “synapse organization”. Furthermore, we identified pathways, e.g., “glutamate signaling”, and genes, encoding, e.g., an angiotensin receptor subtype or a proprotein convertase, that have not been previously linked to PD. The identification of 33 genes that were common among the SN and putamen DEGs, which included the α-synuclein paralog β-synuclein, may contribute to the understanding of general PD mechanisms. Thus, our proof-of-concept data highlights new genes, pathways and lipids that have not been explored before in the context of PD.


Introduction
Parkinson's disease (PD) is the second most common neurodegenerative disease, after Alzheimer's disease. PD is rare in individuals less than 50 years old, but its prevalence increases with age, reaching up to 4% in older-age groups [1]. Clinically, PD presents with motor symptoms, e.g., rigidity, bradykinesia, tremor and postural instability [2], and non-motor symptoms, e.g., depression, hyposmia, constipation and sleep disorders [3]. The neuropathology of PD is characterized by the formation of Lewy bodies (abnormal intracellular aggregates containing metals, lipids and proteins like α-synuclein), a progressive loss of dopaminergic neurons that project from the substantia nigra (SN) to the dorsal part of the striatum (which consists of the caudate nucleus and putamen), and microgliosis [4].
Current knowledge of the molecular mechanisms underlying these pathological hallmarks has been derived from studies on 27 monogenic familial forms of PD [5], which represent only 5-10%

Human Samples
Snap-frozen post-mortem human brain samples were obtained from the Netherlands Brain Bank (NBB; Project 870; Table 1) and all material and data collected by the NBB were obtained on the basis of written informed consent. The average age for control and PD individuals was 78.3 and 77.9 years old, respectively. The average postmortem delay (time interval from death to sample-freezing) for control and PD individuals was 6.49 and 5.53 h, respectively. The cerebrospinal fluid pH for control and PD individuals was 6.54 and 6.43, respectively. None of the three variables were significantly different between the groups (i.e., age: PD and controls (p-value > 0.999), male and female controls (p-value = 0.224), male and female PD (p-value = 0.386), males PD and controls (p-value = 0.892) and female PD and controls (p-value = 0.800); postmortem delay: PD and controls (p-value = 0.184), male and female controls (p-value = 0.191), male and female PD (p-value = 0.829), males PD and controls (p-value = 0.474) and female PD and controls (p-value = 0.200) and pH: PD and controls (p-value = 0.134), male and female controls (p-value = NA), male and female PD (p-value = 0.476), males PD and controls (p-value = 0.474) and female PD and controls (p-value = NA)). The cases reported here have not been genetically analyzed for familial mutations (5-10% of Parkinson cases are thought to be familial; [28]). All procedures were in accordance with the consensus criteria established by the Netherlands Brain Bank. Tissue blocks were cold grinded in liquid nitrogen to obtain a homogenous sample for both RNA and lipid extraction. The table includes the diagnosis (control or Parkinson's disease (PD)), age at death, gender (male (M) or female (F)), post-mortem delay (PMD), pH of the cerebrospinal fluid (CSF), Braak stage and cause of death.

Lipid Profiling
Lipid profiling concerned a validated method for reproducible and high-throughput lipid analysis as described in [29][30][31]. Briefly, pellets of grinded tissue samples were mixed with UPLC-grade chloroform: methanol 1:1 (v/v) and, after 20 min, samples were centrifuged at 2000× g and the supernatant was used directly for LC-MS analysis. To this end, 10 µL was injected on a hydrophilic interaction liquid chromatography (HILIC) column (2.6 µm HILIC 100 Å, 50 mm × 4.6 mm, Phenomenex, Torrance, CA, USA), and eluted with a gradient from ACN/Acetone (9:1, v/v) to ACN/H2O (7:3, v/v) with 10 mM ammonium formate, and both with 0.1% formic acid. Flow rate was 1 mL/min. The column outlet of the LC was either connected to a heated electrospray ionization source of an LTQ-XL mass spectrometer or a Fusion mass spectrometer (both from ThermoFisher Scientific, Waltham, MA, USA). Full scan spectra were collected from m/z 450-950 at a scan speed of 3 scans/s in both the positive and negative ionization mode (LTQ-XL). On the Fusion mass spectrometer full spectra were collected in the negative ionization mode from m/z 400 to 1600 at a resolution of 120,000. Parallel data dependent MS2 was done in the linear ion trap at 30% HCD collision energy. Retention times of lipid classes were verified by comparison to authentic standards and response factors of lipid classes were calculated. Care was taken to ensure that lipid concentrations did not exceed the linear response range of the instrument.

RNA-Sequencing Analysis
RNA was extracted using the RNeasy Lipid Tissue Mini Kit (Qiagen, Cat No./ID: 74804, Venlo, The Netherlands). Total RNA samples were spectrophotometrically analyzed, and their 260/280 ration was typically above 1.9. RNA-seq directionality library preparation and sequencing with Poly(A) selection was performed by HudsonAlpha Genomic Services Laboratory (Huntsville, AL, USA). Samples underwent sequencing on the HiSeq v4 (PE, 50 bp, 25 M reads), per standard protocols. Sixteen samples were run over two lanes (3.125 M reads/sample). Data is deposited under the accession number GSE136666.
Quality control from raw data was performed with the FASTQC tool (https://www.bioinformatics. babraham.ac.uk/projects/fastqc/, Babraham Institute, Cambridge, UK). Overall, data had enough Cells 2020, 9, 1966 4 of 22 quality to follow with transcriptome alignment. Reads mapping and quantification was performed with the Salmon pipeline [32], which consists of a lightweight-mapping model, an online phase that estimates initial expression levels and model parameters, and an offline phase that refines expression estimates. Next, samples with less than 10 reads across samples were filtered out. Quality controls, including sample distances and principal component analysis, performed after variance stabilizing transformation, concluded that all samples could be included for the differential expression analysis. Differential expression analysis was performed with the DESeq2 v1.22.2 Bioconductor's package [33,34]. Finally, clusterProfiler [35] was used to analyze and visualize functional profiles.

Compilation of Transcriptomic Studies
The terms used to browse the GEO DataSets database included "Parkinson's disease" and "Substantia Nigra" or "Putamen". GEO DataSets was applied to screen the studies compiled by [25]. Thirteen studies on human postmortem SN and/or putamen samples were found ( Table 2) and they were used for the validation of our data, except for the study that distinguished between medial and lateral SN (GSE8397) was excluded. Data includes GEO Data Sets accession number, tissue, number of controls and PD patients, platform used and inclusion or exclusion criteria. The analysis and visualization of functional profiles was done with clusterProfiler [35]. Finally, we performed protein-protein interaction network analysis of the common differentially expressed genes (DEGs) with STRING, using only the highest confidence interaction score [36].

Gene Set Enrichment Analysis
The bioconductor limma parametric multivariate rotation gene set test (ROAST) [37] was applied to evaluate gene sets associated with the metabolism and/or cellular signaling pathways of the lipid classes deregulated in the SN and putamen. The terms "phosphatidylcholine", "phosphatidylethanolamine", "phosphatidylserine", "phosphatidylinositol" and "sphingomyelin" were searched at AmiGO [38] and Harmonizome [39]. From the results obtained, the gene sets with (1) only the lipid class name, (2) metabolism, biosynthesis and catabolism or (3) binding information, and containing at least 3 genes, were included in the analysis. Additionally, the specific lipid species were searched in the Human Metabolome Database [40] and retrieved the related proteins, which were also included in the analysis. More precise statistical results were obtained by running the analysis with 10,000 random rotations, and significance values were corrected for multiple testing using the Benjamini-Hochberg false discovery rate method.

Experimental Design and Statistical Analyses
In this exploratory proof-of-concept study, twenty post-mortem human SN samples (10 controls and 10 PD patients) and 18 post-mortem striatum samples (9 controls and 9 PD patients), from the same individuals, were analyzed. Age, postmortem interval and cerebrospinal fluid pH differences between groups were assessed by the non-parametric Mann-Whitney U test. Lipid data analysis was performed by combining a univariate analysis and an ensemble learning method, together with a Receiver Operating Characteristic (ROC) analysis to select optimal findings. The univariate analysis included multiple comparisons of all lipids and a selection of those with a p-value < 0.05. Random forests, an ensemble learning method, were used to identify lipids relevant for disease/gender classification. Similar to a previous study [41], the classification at each node in a tree was done with 100 randomly selected metabolites and the forest consisted of 5000 trees. The process was repeated 50 times, and the model with the highest area under the curve (AUC) was kept. Lipids with a mean decreased accuracy >2 were considered relevant. A ROC analysis was performed on all the lipids considered relevant by either method, and an AUC > 0.8 and p-value < 0.05 was established as the threshold for important lipids to distinguish between conditions. Lipids that fulfilled at least two out of the three criteria were considered as associated with the disease/gender.
For RNA-seq analysis, pools of two and three RNA samples of the same gender and condition were made for the SN and putamen samples, respectively. Differential expression analysis was performed with the DESeq2 v1.22.2 Bioconductor's package [33,34]. Since this analysis represents an exploratory study, we chose not to use multiple testing adjustments, as supported by several authors [42], but apply a threshold for effect size. Hence, genes with a p-value < 0.01 and > 0.5849-log2FC were considered as significant.
Regarding the reanalysis of published transcriptomics data, all included datasets were analyzed using GEO2R [43] to obtain a homogeneous analysis. Genes with p-value < 0.01 and log2FC < −0.5849 or log2FC > 0.5849 were considered differentially expressed genes (DEGs). To convert the obtained identifiers to gene symbols when these were not available from the GEO2R analysis, g:Profiler (version e95_eg42_p13_f6e58b9, University of Tartu, Tartu, Estonia) was used [44]. Two studies did not report any DEGs and one study had only upregulated DEGs, and the three studied were therefore excluded.

Lipid Profiling
The lipid compositions of the SN and putamen samples from PD patients and controls were analyzed by LC-MS. We identified 269 phospholipids of the following classes: Bis (Monoacylglycero)Phosphate (BMP), phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylglycerol (PG), phosphatidylinositol (PI) and phosphatidylserine (PS), and the sphingolipid sphingomyelin (SM; Spreadsheet S1). Hierarchical clustering showed a clear separation of the lipid patterns between the samples from the two brain regions, but not a clear division between the samples based on pathology or gender (Figure 1a). In the SN, we found five lipid species that have different abundances between PD patients and controls, namely

Lipid Profiling
The lipid compositions of the SN and putamen samples from PD patients and controls were analyzed by LC-MS. We identified 269 phospholipids of the following classes: Bis (Monoacylglycero)Phosphate (BMP), phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylglycerol (PG), phosphatidylinositol (PI) and phosphatidylserine (PS), and the sphingolipid sphingomyelin (SM; Spreadsheet S1). Hierarchical clustering showed a clear separation of the lipid patterns between the samples from the two brain regions, but not a clear division between the samples based on pathology or gender ( Figure 1a). In the SN, we found five lipid species that have different abundances between PD patients and controls, namely   Percentage of abundance of lipid species that differ between PD patients and controls, according to at least two of the three selection criteria, in (b-f) SN and (g-i) putamen, and between male and female controls, and also in (j) SN and (k) putamen of male PD patients compared to controls, or in (l) putamen of female PD patients compared to controls. Percentage of abundance of a lipid that (m) has different abundances in both males and female controls, and male and female PD patients, together with a difference in female PD patients compared to controls, and that is lower in females than in males, both in controls and PD patients, in (n) SN and (o) putamen. BMP, Bis(Monoacylglycero)Phosphate; PC, phosphatidylcholine; PE, phosphatidylethanolamine; PI, phosphatidylinositol; PS, phosphatidylserine and SM, sphingomyelin, where d18:1 indicates the common sphingosine base c.f. standard nomenclature [45]. Sample values are plotted individually. 'A indicates that the sn-1 position is linked via an ether linkage. Numbers (X:Y) represent the number of carbons (X) and double bonds (Y). # Differences between control and PD samples, within one gender; * differences between male and female samples, within one condition. Percentage of abundance of lipid species that differ between PD patients and controls, according to at least two of the three selection criteria, in (b-f) SN and (g-i) putamen, and between male and female controls, and also in (j) SN and (k) putamen of male PD patients compared to controls, or in (l) putamen of female PD patients compared to controls. Percentage of abundance of a lipid that (m) has different abundances in both males and female controls, and male and female PD patients, together with a difference in female PD patients compared to controls, and that is lower in females than in males, both in controls and PD patients, in (n) SN and (o) putamen. BMP, Bis(Monoacylglycero)Phosphate; PC, phosphatidylcholine; PE, phosphatidylethanolamine; PI, phosphatidylinositol; PS, phosphatidylserine and SM, sphingomyelin, where d18:1 indicates the common sphingosine base c.f. standard nomenclature [45]. Sample values are plotted individually. 'A indicates that the sn-1 position is linked via an ether linkage. Numbers (X:Y) represent the number of carbons (X) and double bonds (Y). # Differences between control and PD samples, within one gender; * differences between male and female samples, within one condition.
Since sex hormones have been found to modulate brain lipid levels [46], we performed the analyses separately for males and females (Spreadsheet S1). Interestingly, a number of lipid species displayed different abundances in female and male controls, including PC 36:4 in the SN (Figure 1j), and PI 34:1 and PS 38:3 in the putamen (Figure 1k,l), while these differences were not present in PD patients. Additionally, some lipid species had distinct levels in female and male controls and PD, e.g., PE A38:4 in the putamen (Figure 1m), whereas PS A38:4 levels were lower in females than in males, regardless of the condition and in both brain regions (Figure 1n,o).

Transcriptomic Profiling
We performed RNA-seq analysis of SN and putamen samples from PD patients and age-matched controls. Hierarchical clustering showed a clear separation of the transcriptomic profiles between the samples from the two brain regions, as well as between the male and female samples (Figure 2a).

Transcriptomic Profiling
We performed RNA-seq analysis of SN and putamen samples from PD patients and agematched controls. Hierarchical clustering showed a clear separation of the transcriptomic profiles between the samples from the two brain regions, as well as between the male and female samples ( Figure 2a). , including also the present study; interaction networks of proteins encoded by genes (h) downregulated or (i) upregulated in at least two of the eight studies included in the overlap analysis. Protein-protein networks were obtained with STRING [31]. Only the highest confidence associations were included. Non-connecting nodes were removed. , including also the present study; interaction networks of proteins encoded by genes (h) downregulated or (i) upregulated in at least two of the eight studies included in the overlap analysis. Protein-protein networks were obtained with STRING [31]. Only the highest confidence associations were included. Non-connecting nodes were removed.

DEGs in the SN of PD Patients and Controls
We then performed a restrictive analysis of the RNA-seq data (adjusted p-value < 0.05 and >0.5849-log2FC) and identified 39 upregulated and 41 downregulated transcripts (Spreadsheet S2). Among the DEGs encoding these transcripts, 7 represented unknown genes (no gene symbol associated), 2 were uncharacterized transcripts (LOC) and 3 were genes for long intergenic non-protein coding RNAs Next, we performed a less-restrictive RNA-seq data analysis (p-value < 0.01 and > 0.5849-log2FC) to determine which biological pathways were enriched in the dataset. This analysis yielded 354 DEGs in the SN of PD patients compared to controls, of which 152 were upregulated and 202 downregulated (Spreadsheet S2). To understand the biological significance of this outcome, we analyzed on the basis of gene ontology (GO) terms which biological pathways were enriched in the 354 DEGs as well as in the upand down-regulated DEGs separately (Spreadsheet S2). The top-three pathways in the all-DEGs analysis Cells 2020, 9, 1966 8 of 22 were "signal release" (GO:0023061), "neurotransmitter transport" (GO:0006836) and "amine transport" (GO:0015837; Figure 2b). The pathways "response to heat" (including GO:0009408 and GO:1900034) and "chaperone-mediated protein folding" (GO:0061077) represented the top significantly enriched pathways in the upregulated genes (Figure 2c), while the biological pathways related to "regulation of neurotransmitter levels" (GO:0001505), "signal release" (GO:0023061) and "neurotransmitter transport" (GO:0006836) were the top significantly enriched ones among the downregulated DEGs (Figure 2d).
Subsequently, we validated our results by comparison with the results of the only other RNA-seq study on PD and control SN [47]. We found that 9.8% of the DEGs identified in our study were also listed as differentially expressed in [47]. However, while 47% of the overlapping genes were differentially expressed in the same direction, the remaining 53% were differentially expressed in the opposite direction, questioning the validity of the overlap. We then aimed to validate our results and obtain more robust outcomes by comparing the publicly available transcriptomic studies performed with microarrays, homogeneously analyzed with GEO2R, and our data. Interestingly, our data partially overlaps with the results of these previous studies (Table 3). We found that 311 DEGs (112 upregulated and 199 downregulated) were present in the same direction in at least two datasets (overlapping genes; Spreadsheet S3). The top biological pathways enriched in the 311 DEGs, and also in the downregulated DEGs, were associated with "synapse organization" (e.g., GO:0050808), while the top pathways enriched in the upregulated DEGs were "regulation of cellular response to heat" (GO:1900034) and "chaperone-mediated protein folding" (GO:0061077; Figure 2e-g). The table includes the accession numbers of the studies used, the number of differentially expressed genes with p-value < 0.01 and log2FC < −0.5849 or log2FC > −0.5849 (DEGs), the number of DEGs overlapping between all studies (overlap), and the number of overlapping DEGs that are either upregulated or downregulated in the same direction in all studies where they appear as DEGs. Percentages of overlapping genes from the total DEGs, and percentages of overlapping genes in the same direction from the total number of overlapping genes, can be found between brackets in the rows "Overlap" and "Same direction", respectively.
Next, we analyzed the interactions between the proteins encoded by the 311 DEGs found in at least two studies (overlap). One of the largest clusters of proteins encoded by the downregulated DEGs (Figure 2h), which includes some of the DEGs common in four out of eight studies, such as TH, DDC, DRD2, KCNJ6 (also known as GIRK2), SLC18A2 (VMAT2) and SLC6A3 (DAT), was related to dopamine synthesis and transport, providing confidence in the effectiveness of our approach. This pertinent cluster was also associated with other proteins linked to the G-protein-coupled receptor signaling (e.g., GNG3). Additional protein clusters linked to downregulated DEGs were associated with GDNF-RET signaling (e.g., RET), and neural development (e.g., SLIT1). The protein-protein interaction networks of upregulated DEGs (Figure 2i) consisted of six main protein clusters associated with multiple cellular processes, including protein folding (e.g., HSPA1A), the aldo-keto reductase family (e.g., AKR1C1), cell migration and extracellular matrix (e.g., ITGB1), DNA (damage) (e.g., ATM) and ubiquitination (e.g., RNF130).

DEGs in the Putamen of PD Patients and Controls
We identified 15 DEGs (5 upregulated and 10 downregulated genes) in the putamen of PD patients compared to controls when using restrictive parameters for RNA-seq data analysis (adjusted p-value < 0.05 and > 0.5849-log2FC). Two of these DEGs encoded unknown transcripts (no gene symbol associated) and 1 lincRNAs RNA, respectively. The five upregulated DEGs produced to two unknown transcripts (ENSG00000244921, log2FC Next, we identified 261 DEGs in the putamen of PD patients compared to controls with the less-restrictive analysis (p-value < 0.01 and > 0.5849-log2FC), of which 68 were upregulated and 193 downregulated (Spreadsheet S4). The top three pathways enriched in these putamen DEGs were "synapse organization" (GO:0050808), "cartilage development" (GO:0051216) and "prostaglandin transport" (GO:0015732). There were no biological pathways enriched in the upregulated DEGs, while the downregulated ones were involved in "synapse organization" (GO:0050808), "cartilage development" (GO:0051216) and "synapse assembly" (GO:0007416; Spreadsheet S4).
There are less transcriptomic studies on putamen than SN from PD patients and controls, and most of them have resulted in the identification of only a very low number of DEGs. Moreover, essentially no overlap between the DEGs has resulted from these studies ( Table 4). We therefore performed no further analysis on these datasets. The table includes the accession numbers of the studies used, the number of differentially expressed genes with p-value < 0.01 and log2FC < −0.5849 or log2FC > 0.5849 (DEGs), the number of DEGs overlapping between all studies (overlap), and the number of overlapping DEGs that are either upregulated or downregulated in the same direction in all studies where they appear as DEGs. Percentages of overlapping genes from the total DEGs, and percentages of overlapping genes in the same direction from the total number of overlapping genes, can be found between brackets in the rows "Overlap" and "Same direction", respectively.

DEGs Common between SN and Putamen
Finally, we analyzed the overlap between SN and putamen mRNA expression profiles to find possible global mechanisms underlying PD. As mentioned above, since there was no robust overlap between previous transcriptomic studies on the putamen, we could only use our own data and found 33 DEGs common between the SN and putamen profiles (Table 5).

Integration of the Lipid and Transcriptomic Profiling Data
Since the lipid and transcriptomic profiling data were obtained from the same samples, we analyzed if gene sets linked to the deregulated lipids (i.e., BMP 42:8, PC 36:3, PE A36:2, PI 42:10 and PS 36:3 in the SN, and PI 34:2, SM d18:1;14:0 and SM d18:1;16:0 in the putamen) were differentially expressed. Accordingly, we performed the targeted gene set enrichment ROAST test, a hypothesis-driven analysis ( Table 6). A weak association (p-value < 0.05 and FDR-adjusted p-value < 0.25) between the upregulated SN transcripts and the dataset "phosphatidylcholines" (GeneRIF), and between the deregulated SN transcripts and the datasets "phosphatidylcholine biosynthetic process" (GO:0006656) and "phosphatidylcholine catabolic process" (GO:0034638) were found. No (weak) association between the deregulated transcripts from the SN and phosphatidylinositol-related gene sets was observed. Furthermore, a weak association between the deregulated SN transcripts and the datasets "phosphatidylethanolamine biosynthetic process" (GO:0006646) and "phosphatidylethanolamine biosynthesis" (HumanCyc Pathways) were found. Additionally, a weak association between the upregulated SN transcripts and "phosphatidylserine binding" (GO:0001786), and between the deregulated SN transcripts and the dataset "phosphatidylserine catabolic process" (GO:0006660) were found. A weak association between downregulated transcripts in the putamen and the datasets "phosphatidylinositol binding" (GO:0035091) and "Sphingomyelin Metabolism/Ceramide salvage" (HumanCyc Pathways) was also found. Unfortunately, retrieval of gene sets associated with the poorly characterized lipid BMP 42:8 was not possible.

Discussion
The data presented here represent the first lipid and transcriptome analyses performed on PD and control SN as well as putamen samples from the same individuals. We conducted a transcriptome analysis by RNA-seq rather than microarrays in order to identify not only annotated transcripts, but also currently unannotated sequences for future re-evaluation (50/354 SN DEGs and 29/261 putamen DEGs represent not-annotated transcripts) [26]. Additionally, we analyzed our data together with publicly available transcriptome datasets to draw more robust conclusions about the relevance of gene functions and pathways for PD.
Lipidomics likely represents the most informative -omics technology complementary to RNA-seq. In this context, it is indeed important to realize that lipids are structural and bioactive molecules essential for multiple brain processes such as myelination, synaptic signal transduction and acting as second-messenger precursors and, relative to water-soluble metabolites, generally have a much slower turnover and are thus less prone to post-mortem degradation. Previously, neuropathological lipidomic studies have been conducted on primary visual cortex [48], primary motor cortex [49], hippocampus [50] and cerebellum [51] from PD and/or Lewy body disease individuals. Our exploratory lipid profiling involved the analysis of seven classes of lipids and showed several lipid species modulated in the SN and putamen of PD patients compared to age-matched controls. Two other studies have been performed on the lipid profile of the SN [52,53] and one on the putamen [51] of PD patients, but they have reported only differences in lipid classes rather than lipid species and focused on sphingolipid and gangliosides, respectively, hindering a comparative analysis. The differences we found in the putamen include PI 34:2 and two saturated SM species, namely SM d18:1;14:0 and SM d18:1;16:0. Unfortunately, the limited current knowledge of putamen lipids does not allow a functional interpretation of the putamen lipid changes we detected. Lipids BMP 42:8 and PI 42:10, both with increased abundance in the SN of PD patients compared to controls, are thought to have arachidonic acid (20:4) as one of the two side chains, while the three decreased lipid species, PC 36:3, PE A36:2 and PS 36:3, are likely to have linoleic acid (18:1) as one of their two side chains. Although such specific changes in lipid composition are difficult to directly link to distinct enzymatic activities in specific metabolic conversion steps, it is tempting to speculate on their wider biological implications. The presence of the five differentially expressed SN lipid species could point towards a neuroinflammatory component in disease etiology, since increased arachidonic acid has been demonstrated in acute neuroinflammation [54], while linoleic acid is one of the sources of arachidonic acid [55]. A neuroinflammatory involvement is also in line with the findings of our transcriptomic analysis revealing DEGs associated with the transport of prostaglandins, which regulate neuroinflammatory pathways. Lipid BMP 42:8 belongs to a lipid class highly enriched in lysosomal membranes [56] and is key for the proper functioning of lysosomal hydrolases [57]. Interestingly, BMP accumulation also occurs in macrophages from patients with Gaucher disease [58], a rare genetic disorder characterized by deposition of the lipid glucocerebroside in macrophages and that highly increases the life-time risk to develop PD [59]. Nevertheless, we realize that it is difficult to draw functional conclusions concerning individual lipids since the number of samples employed in this proof-of-concept study is relatively small and only the roles of lipid classes have been characterized.
Even though our sample size was small, it is important to note that we found gender-related differences in the abundance of lipids both in PD patients and controls. This observation is consistent with the gender-related differences that have been found in the SN lipid profiles of PD patients [53] and the "female pregnancy" pathway detected in our transcriptional analysis. Current data indicates that gender differences exist not only in PD incidence and prevalence [60], but also in its clinical manifestation [61][62][63], with estrogens as a possible mediator for these differences [64], although the underlying molecular mechanisms are unknown. Since estrogens are steroid hormones and thus lipids modulating lipid metabolism in the brain [65], and gender dimorphism has been observed in the fatty acid profiles of control mouse brains [66] as well as in our data, brain lipidome modulation by sex hormones may be one of the underlying protective mechanisms in PD. Clearly, the possible link between sex hormones, brain lipids and PD deserves further attention.
In our exploratory transcriptomic study, the top upregulated DEG in the SN was ORC7C1, which encodes an odorant receptor. Mesencephalic dopaminergic neurons express olfactory receptors and they respond to odorant-like molecules in mice [67]. Several odorant receptors have been found to be downregulated in the frontal cortex of PD patients [68], but the difference in the directionality of expression (up-or down-regulation) could be brain region related. The fact that the top upregulated DEG as well as the 7th upregulated DEG (i.e., OR7A5) represent odorant receptors highlights the importance of a poorly studied group of receptors that could modulate the cell behavior and fate via binding of small molecules and, thus, could play a crucial role in PD. The relationship between the upregulated MTRNR2L8 and the mitochondrial MT-RNR2 is at present unclear [69] and its role in PD has not been analyzed yet. The protocadherin encoded by PCDH20, which is thought to be involved in the establishment and function of specific cell-cell connections in the brain, has been found to be also upregulated in a meta-analysis of transcriptomic studies from SN of PD patients [70]. Nevertheless, the role of this protein has thus far not been studied in the context of PD. The transcription factor KLF5 binds to GC box promoter elements and activates their transcription, mediating multiple cellular processes, such as proliferation, migration and differentiation, among others [71], but its significance for PD is currently unknown.
The top-downregulated DEG in the SN is TPH2, which encodes the rate-limiting enzyme in the synthesis of serotonin. Interestingly, serotoninergic neurons are progressively lost in PD, which is implicated both in motor and non-motor manifestations of PD [72]. Moreover, polymorphisms in TPH2 are associated with addictive behaviors [73] and depression [74] in PD patients. Additionally, Tph2 KO mice, which have serotonin deficiency, present with systemic oxidative stress and lipidomic abnormalities [75], and swallowing dysfunction [76]. The alpha-D1 adrenergic receptor encoded by ADRA1D is also downregulated in the hippocampus of Alzheimer's disease and dementia with Lewy bodies patients [77], and noradrenergic impairment is present in PD patients as well [78]. Both findings highlight the importance of neurotransmitters other than dopamine in the pathology of PD. It is therefore imperative to do further studies into their role to obtain a better understanding of the complexity of disease. The observed downregulated expression of IL1B is in contrast with earlier findings showing increased IL1B expression in striatum [79] and cerebrospinal fluid [80] of PD patients. However, the complex effects of IL1B seem to depend on the time, place and level of IL1B expression, e.g., infusion of IL1B in the striatum of rats 5 days before 6-OHDA injection had a protective effect [81]. Finally, DOK7 encodes a protein involved in neuromuscular synaptogenesis, but its role in PD pathogenesis remains to be established.
PD is characterized by the loss of dopaminergic neurons from the SN [4], which reaches a reduction of at least 70% by the time the disease is diagnosed. Thus, our top finding of decreased expression of genes associated with neurotransmitter levels, and dopamine synthesis and transport in postmortem PD SN samples is likely linked to the decreased number of dopaminergic neurons and fibers in the SN and putamen, respectively. One of these was the PD-associated RET gene [82], encoding a GDNF receptor required for the preventive and compensatory mechanisms linked to dopaminergic system degeneration that is triggered by the neurotrophic factor [83,84]. Remarkably, RET and some of its interaction partners (GFRA1, DOK4 and DOK6) are downregulated in multiple SN transcriptomic studies, highlighting the importance of the process of dopaminergic neurodegeneration in PD brains. The two most-enriched upregulated pathways were "cellular response to heat" and "protein folding", which have an overlap of 80% among their DEGs. Since misfolded proteins are a hallmark of PD [85], upregulation of genes associated with these pathways, such as the molecular chaperone HSPA1A, may represent a compensatory cellular mechanism to handle the accumulation of misfolded proteins. Furthermore, we confirmed other previously reported PD mechanisms, such as a dysfunction of protein folding. Interestingly, other SN DEGs comprise novel PD-related pathways that may provide a better understanding of the disease pathology, including "positive regulation of acute inflammatory response" (GO:0002675), "female pregnancy" (GO:0007565) and "G-protein-coupled glutamate receptor signaling pathway" (GO:0007216).
The overlap of our RNA-seq data with results from earlier microarray studies on postmortem PD and control SN confirms the importance of the degeneration of dopaminergic neurons in this brain region, since 10 out of the 12 DEGs that were found in at least half of the studies are associated with dopaminergic neurons or PD (ALDH1A1, DDC, DLK1, DRD2, GAP43, KCNJ6, SLC18A2, SLC6A3, SV2C and TH). The other two genes are AGTR1, the angiotensin II receptor subtype AT1, which has been shown to be downregulated in 5 out of the 8 SN studies, and the neuroendocrine proprotein convertase 1, PCSK1, which was downregulated in 4 out of the 8 SN studies. AGTR1 and PCSK1 do not have a known role in dopaminergic neurons or PD and are therefore interesting candidates for future functional studies. The reason for the restricted overlap between the results of our study and those of the only other PD/control SN RNS-seq analysis is at present unclear.
The top upregulated protein-encoding transcript in the putamen is MYOT, which codes for a component of a complex of actin cross-linking proteins. Transcript levels of MYOT have been found to be upregulated and downregulated in the SN and blood of PD patients, respectively [86]. Furthermore, PD patients present a higher prevalence of MYOT autoantibodies than controls [87]. Thus, MYOT may play distinct roles in the pathology of PD, but the nature of its involvement has not been elucidated. Neither have the function of the transcription factor ZNF646 or its link to PD been widely studied. The upregulated gene HP1BP3 mediates chromatin condensation and modulates cognitive aging [88], but not in known association with PD.
The top downregulated transcript in the putamen is SLC30A3, which encodes a protein involved in the accumulation of zinc in synaptic vesicles, rather than in the cytosol. Decreased levels of the protein have been observed in the (pre)frontal cortex of DLB and PDD patients [89,90] and is thus associated with dementia. Moreover, PD patients present changes in the cellular distribution of zinc in the SN [91] and treatment of midbrain primary cultures with MPP+ leads to zinc ion accumulation [92]. These findings highlight the relevance of bivalent metal ions in the pathology of PD and encourage further studies in this direction. The second most downregulated transcript was DUP2, a phosphatase involved in the negative regulation of ERK1 and ERK2. Members of the ERK family of proteins are known to control cellular proliferation and differentiation, but association of this function with PD has not yet been found. Similarly, GNG2 encodes a gamma subunit of heterotrimeric G proteins, which are involved in cellular responses to external signals, but the significance of this function for PD still needs to be clarified. PPL codes for a component of desmosomes and is thus involved in filament binding. Although PPL has not been associated with PD, genetic variants of another Plakin family member, MACF1, lead to decreased expression and confer risk for PD [93]. Finally, GRM2 encodes a metabotropic glutamate receptor and, interestingly, glutamate-induced excitotoxicity has been suggested to result in the loss of dopaminergic cells in PD [94]. Moreover, activation of both metabotropic glutamate receptors 2 and 3, which have inhibitory functions, gives rise to striatal protection in PD rodent models [95]. Hence, this finding can be considered further support for the involvement of multiple neurotransmitter systems in PD.
Analysis of data from various transcriptomic studies on putamen showed no overlapping genes. While SN is considered to be one of the primary sites of degeneration in PD, the neurodegeneration in the dorsal striatum may be mainly a consequence of the depletion of dopaminergic innervation [96]. Thus, the lack of overlap among putamen DEGs could point to more heterogeneous pathology in the putamen than in the SN. We found an enrichment of downregulated putamen genes that play a role in the control of neurotransmitter levels and transport, which may well be associated with the loss of dopaminergic innervation [96]. Furthermore, we observed DEGs associated with prostaglandin transport, which may argue for a neuroinflammatory response in this brain region [97]. Indeed, neuroinflammatory events may play a role in nearly all known neurodegenerative diseases [98].
Interestingly, we found 33 genes that were differentially expressed in both the PD SN and putamen, hinting at mechanisms shared by the two brain regions. The use of pathway analysis software nor the findings from extensive literature searches led to the identification of a specific pathway within the set of 33 genes. Of note is the downregulation of β-synuclein, the paralog of α-synuclein, that inhibits the assembly [99], aggregation [100] and toxicity [101,102] of α-synuclein. Since α-synuclein aggregation is one of the neuropathological hallmarks of PD [103], β-synuclein upregulation may represent a potential target to ameliorate the disease. Hence, overexpression of β-synuclein in cellular and animal PD models would be central to better understand its role in the disease. Other potentially interesting examples of the set of common SN and putamen DEGs include DDIT4, which mediates apoptosis in cellular PD-models [104], and P2RX7, the antagonist of which is neuroprotective in cellular and animal models for PD [105]. Additionally, SNPs in the common genes encoding EDN1 and MYO5B have been associated with multiple systems atrophy (which presents with parkinsonism [106]) and loss of the sense of smell (an early non-motor PD-symptom [107]), respectively. Therefore, functional studies on the genes that are differentially expressed in both brain regions, e.g., by manipulating their expression in cellular and animal PD models, may provide clues for further understanding of the pathobiological complexity of PD.
The parallel analysis of the lipid and transcriptomic profiles of the same samples allowed us to integrate the two datasets. In this connection, one has to realize that the characterization of lipid species is limited, while the proteins involved in the metabolism of lipid classes have been well studied. Yet, targeted gene set enrichment analysis revealed that gene sets associated with the lipid classes differentially expressed in the SN and putamen were deregulated in the respective brain region. Thus, the lipidome and the transcriptome appear to be tightly connected, and transcriptomic changes leading to differences in the lipid profile may be one of the molecular mechanisms underlying neurodegeneration in PD. Still, additional studies involving multiple -omics approaches are necessary to allow a better characterization of the role of lipids in PD. Such studies will be of particular interest since several familial PD genes, such as LRRK2 [108], PINK1 [109], DJ1 [110], ATP13A2 [111], PLA2G6 [112] and ATXN2 [113], are known to modulate various lipid metabolism pathways.
This study has several limitations. First, relative to control samples, the SN samples from PD patients are characterized by an extensive loss of dopaminergic neurons and gliosis, which hinders making a distinction between lipid/mRNA changes caused by the disease-generated differences in the composition of cellular populations or by the molecular mechanisms leading to neurodegeneration (i.e., the "consequence" or the "cause"). Second, our lipid analysis was focused on phospholipids and sphingolipids, and did not include other lipid classes that could be relevant for PD, such as sterols and glycolipids [24]. Third, we pooled samples for RNA-seq analysis, which hampered a gender-dependent analysis as performed with the lipid profiling, although the hierarchical clustering showed a clear separation by gender. Moreover, the pooling of samples might have masked minor differences in mRNA expression levels due to a small sample size and thus less statistical power.

Conclusions
We here reported the first combined lipid and transcriptomic analyses of SN and putamen samples from the same PD patient and control individual. This exploratory proof-of-concept study allowed us to unravel changes in the PD brain lipid profile, some of which were gender dependent and should be studied further to evaluate their relevance. Additionally, we confirmed molecular pathways associated with PD and identified several potentially interesting novel genes and pathways that may help to better understand PD pathology.