Study designIdentifying causal genes that can be further considered as potential therapeutic targets remains a complex undertaking24. Characterizing causal variants that underpin genetic associations with a disease, as well as elucidating the mechanistic impacts of genetic aberrations on the human proteome and normal protein interactions necessitates integrative approaches to combine genetics with other omics data types such as proteomics and epigenomics25,26. Moreover, since majority of drug targets are proteins, then proteome-based Mendelian randomization (MR) as well as modeling alterations in the human interactome caused by genetic abnormalities is of great importance in terms of identifying novel disease causal genes/proteins27,28. In this study, we designed a three-phase framework (Fig. 1) to characterize the impact of PD genetic variants on proteins followed by a series of network-based analyses, epigenomic and transcriptomic analyses aimed at creating a comprehensive mapping of interconnections between genome and proteome and systemically integrate such data-types to identify proteins that might be associated with PD. Such a systems-level design using orthogonal data from the same cohort of patients allowed us to detect not only direct genetic risk factors, but potential indirect mediators for which no sentinel genetic variants were identified that could play critical roles in disease etiology and pathology. Such mediators may either be proteins identified through causal inference analysis or through network analysis. The designed workflow consists of (1) proteogenomic analysis; (2) network analysis; (3) transcriptome/epigenome data analysis. The outcome of the proteogenomic analysis is subsequently utilized in the following phases.Fig. 1: An overview of the study design.The study consists of three stages, including proteogenomic analysis, network analysis, and epigenomic data analysis.The proteogenomic analysis consists of six steps as follows: (1) genome-wide association analysis (GWAS) in UKB on PD; (2) meta analyzing the summary statistics of the UKB GWAS with the GWAS results from FinnGen; (3) protein quantitative trait loci (pQTL) mapping across 54,219 individuals and 2923 proteins in UKB; (4) causal inference; (5) characterizing the enriched biological pathways by the identified proteins followed by determining the pathway specificity of the PD risk variants; (6) aggregating the identified risk variants and generating gene-based PD association scores to be used in network analysis.Using the gene-based scores from the proteogenomic analysis phase, we leveraged a publicly available protein-protein interaction network (PPI)29 to run a series of network analyzes including: (1) a hypothesis-free approach, where we mapped the gene-based scores onto the PPI, conducted an exhaustive search to identify PD-associated risk modules, and paired the most significant modules with pQTLs to pinpoint dysregulated proteins; (2) a hypothesis-driven approach, where we mapped the genetic and proteomic signals to identify pQTL-associated proteins interacting with the genetics-derived risk proteins; (3) characterizing the enriched biological pathways by the identified PD-associated protein complexes. Finally, in the third phase, we conducted epigenomic and single-cell transcriptomic analyzes to investigate the regulatory roles of the identified genetic signals as well as determining the cell-specific transcriptional patterns of PD risk factors in the human brain.In the following sections, we will describe the details of each step of this study design along with the findings.Proteogenomic analysisUsing UKB clinical and imputed array data, we created a cohort of 2,864 clinically defined PD patients and 158,876 controls of European (EUR) ancestry and conducted GWAS using REGENIE v2.2.430 (Methods, Supplementary Data 1). A list of top identified significant associations are provided in Fig. 2a and Table 1. Details of the conducted GWAS are provided in the Supplementary Methods.Fig. 2: Manhattan plots for PD case-control GWAS using the EHR-based clinically defined cohorts from the UKB and the UKB-FinnGen meta-analysis in European samples.(a) UKB case-control GWAS, (b) UKB-FinnGen GWAS meta-analysis. A solid blue line indicates the significance threshold P < 5×10-8. The dashed line shows P = 10-6.Table 1 Genome-wide significant loci from GWAS performed on the PD cohorts in UKB EUR*We next conducted a meta-analysis (Methods, Supplementary Figs.) of the UKB-curated cohort and FinnGen R9 (Table 2, Fig. 2b). FinnGen R9 contains 4,235 PD cases and 373,042 healthy controls leading to a total 7,099 PD cases and 531,918 healthy controls in this study. Overall, 11,114,080 SNPs shared between the two datasets were evaluated for meta-analysis. Details of the performed meta-analysis is provided in the Supplementary Methods. Meta-analyzing UKB and FinnGen has added significant power to replicating many of known PD signals. The FinnGen GWAS signals itself shows two significant loci UBQLN4 and ARL17B in chromosomes 1 and 17, respectively. Upon meta analyzing it with the UKB GWAS, four other significant loci reached genome-wide significance including HIP1R, HLA-DRB1, GBA, and TMEM175. These established loci have been identified in the literature mainly in EUR population. Therefore, meta-analyzing UKB with FinnGen not only did not lead to spurious associations, but increased the number of loci reaching genome-wide significance. This can be attributed to similarities in the genetic architecture of the samples from the UKB and participants of FinnGen. We should note that all the samples in the UKB GWAS were selected from British ancestry.Table 2 Genome-wide significant loci from the meta-analysis of the UKB and FinnGen PD cohortsaTo investigate if any of the identified significant loci are supported by other studies, we used the largest available meta-analysis GWAS by Nalls et al.1. Focusing on the top significant variants in UKB GWAS and the meta-analysis GWAS, seven variants have been replicated, including (a) UKB GWAS: rs1474055-T (intergenic to STK39), rs356203-T (intronic to SNCA), rs35749011-A (intergenic to KRTCAP2), and rs34637584-A (missense to LRRK2); (b) UKB/FinnGen meta-analysis: rs4630591-T (intronic to KANSL1), rs34311866-C (missense to TMEM175), and rs34637584-A (missense to LRRK2). Given the low allele frequency of rs34637584 in the imputed array data (MAF EUR < 0.001), we also checked Gene Bass31 and verified findings through analysis of whole exome sequencing data (OR = 3.37, P = 3.81E-8). The 95% credible set is provided in Supplementary Data 1.We acknowledge that the conducted meta-analysis sample size is modest compared to existing studies. However, our objective was to align genetic associations using orthogonal data sets generated from the same samples and proteome. Moreover, in order to generate gene-based scores to be used in the network analysis, meta-analysis was performed so that we later used the generated scores to weight the nodes of the PPIs. This became possible by leveraging high-throughput proteomic data from the UKB-PPP. Within our created case/control cohort used for running GWAS in the UKB, 638 PD patients and 15,522 healthy controls had proteomic data. In order to integrate PD genetic association with protein interactome for downstream candidate target identification, we used MAGMA32 to generate per-gene scores based on the GWAS meta-analysis results (Supplementary Data 2, Methods). We used the SNPs in the body of each gene to calculate a gene-based p-value associated with PD, which was then used in the network analysis.Using results generated by the UKB-PPP consortium, we evaluated if any of the Parkinson’s disease susceptibility variants also act as pQTLs in blood plasma (Methods). Briefly, the UKB-PPP results comprise 14,287 primary genetic associations derived from 2,923 plasma proteins measured across 54,219 individuals14. Details of proteomic data generation and downstream quality control are outlined in [14]. We intersected all GWAS significant variants from our meta-analysis of UKB and FinnGen with all UKB-PPP pQTL results. In total, 103,118 associations between genetic signals and protein abundances were found, translating to a total of 577 proteins (Supplementary Data 1) where 10% of the observed associations were cis (Fig. 3a). These proteins may not necessarily be disease-associated, thus we further conducted causal inference analyses in the MR and colocalization section to nominate putative causal proteins.Fig. 3: Chromosomal positions of the identified pQTLs relative to their associated proteins.a Position of pQTLs and their associated proteins; b PD-associated SNPs and their associated proteins.Focusing on the GWAS meta-analysis top signals (Table 2), seven PD-associated SNPs were associated with a total of 163 proteins, including rs2230288, rs34311866, rs2760980, rs10847839, rs2183082, rs4630591, and rs6857, among which 4 SNPs were only associated with one protein including rs2230288-T (Ephrin, EFNA; P = 2.4E-10), rs34311866-C (Alpha-L-Iduronidase, IDUA; P = 4.6E-8), rs10847839 (Huntingtin-interacting protein 1-related protein, HIP1R; P = 6.7E-30), and rs2183082-T (Galectin-3, LGALS3; -log10(P) = 311.6). Three other SNPS including rs2760980-A, rs4630591-T, and rs6857-T were associated with 95, 21, 50 proteins, respectively (Fig. 3b). All pQTLs associating with a single protein were in cis.We sought to investigate the enrichment(s) of those proteins across distinct biological pathways. Accordingly, we conducted pathway enrichment analysis (Methods) on 577 proteins associated with GWAS signals and focused on the top 10 significantly enriched pathway at false discovery rate (FDR) of less than 0.05 (Supplementary Fig. 1). We found several PD-related mechanisms to be enriched in the identified proteins, including lysosome (FDR = 1.6E-6), cytokine-cytokine receptor interaction (FDR = 2.2E-16), and chemokine signaling pathway (FDR = 1.2E-6). Details of the enriched pathways are provided in Supplementary Data 3. Among the observed pathways, recent studies have been increasingly focusing on the role of lysosomal genes as risk factors for idiopathic PD33,34. We identified 21 proteins enriched for lysosome pathway including cathepsins (CTSC, CTSD, CTSE, CTSF, CTSL, CTSO, CTSV, CTSZ), ARSA, GLB, GUSB, IDUA, IGF2R, LAMP3, NAGA, NAGPA, PLA2G15, PSAPL1, SCARB2, SGSH, and SMPD1. In addition to known PD-related pathways, majority of the other significant pathways were enriched for inflammatory markers including chemokines such as chemokine ligand (CCL) and CXCL family as well as cytokines such as tumor necrosis factor (TNF) and interleukin (IL). ‘Rheumatoid arthritis’-related pathways were also predominantly enriched for inflammatory markers. This is in line with the existing knowledge of multiple shared genetic variants between PD and rheumatoid arthritis and several other autoimmune diseases35,36. We point out that the UKB-PPP is a population-based study with proteomic evaluation based on a commercial panel of 2923 proteins, therefore many known PD-related proteins were not profiled in the consortium. The proteins in the Olink panel include low-abundant inflammatory proteins, proteins actively secreted into blood circulation, approved and ongoing drug targets, organ-specific proteins leaked into blood circulation, and proteins representing exploratory biomarkers. As a result, some known PD-related pathways may not show significant enrichment as the relevant proteins were not assayed. Moreover, we conducted an enrichment analysis on the SNPs that were found to be pQTLs in the UKB-PPP using MAGMA32 to account for potential confounding factors caused by linkage disequilibrium (LD). We observed a few enriched pathways including cytokine-cytokine receptor interaction (FDR = 1.7E-5), chemokine signaling pathway (FDR = 2.7E-7), lysosome (FDR = 1.5E-7), and rheumatoid arthritis (FDR = 2.3E-5). These observations are in accordance with our pathway enrichment results applied to the pQTL-associated proteins.The UKB-FinnGen meta-analysis is restricted in terms of the sample size. Therefore, we sought to leverage the summary statistics provided by Nalls et al.1 and evaluate their potential associations with plasma protein concentrations in the UKB-PPP. Among the 107 SNPs provided by Nalls et al., 51 SNPs were pQTLs in the UKB-PPP associated with a total of 204 proteins from 287 pQTL-protein associations (Supplementary Data 1). Out of 204 identified proteins, 149 proteins were shared with the 577 PD pQTL-associated proteins from UKB-FinnGen, while 55 proteins were not identified based on UKB-FinnGen GWAS. pQTL-associated proteins are the proteins whose abundance is significantly associated with the SNPs identified in the UKB-FinnGen meta-analysis. We conducted a pathway enrichment analysis on the 204 proteins leading to a total of 6 significant pathways (Supplementary Data 3). Among these, three pathways were shared with the enriched pathways from the PD associated proteins from UKB-FinnGen, including lysosome (FDR = 4.62E-4), graft-versus-host disease (FDR = 0.003), and cell adhesion molecules (CAMs) (FDR = 0.03). The other three pathways that are unique to this analysis include antigen processing and presentation (FDR = 2.88E-5), staphylococcus aureus infection (FDR = 3.94E-4), and natural killer cell mediated cytotoxicity (FDR = 0.03).We conducted a Mendelian Randomization analysis to identify potential PD causally-linked proteins. We used 6381 identified cis-pQTLs from the UKB-PPP as the instrumental variables and protein abundance as the Exposure in this analysis (Supplementary Data 1). We took advantage of summary statistics from Nalls et al.1 to use in the MR analysis as the outcome (Methods). We identified a total of 122 proteins showing a significant association at FDR < 0.05. We also tested for horizontal pleiotropy on the MR results, where no significant pleiotropic effect was observed. Following MR, we conducted a colocalization analysis between the UKB-FinnGen meta-GWAS and pQTLs from the UKB-PPP to characterize if meta-analysis signals colocalize with pQTL hits in the MR-derived proteins aimed at assessing whether the genetic factors from meta-GWAS also colocalize with the genetic predictor of the protein abundances for supporting the MR findings. We found strong colocalization support (H4.PP > 0.9) for a shared causal variant for abundance and PD risk for 17 proteins identified through MR. These 17 proteins are shown in Table 3 and the full list of all MR results as well as sensitivity analysis results are provided in Supplementary Data 1. In addition, we ran a colocalization analysis on the UKB-FinnGen meta-analysis genes (Table 2) and Nalls et al. signals. Among which HIP1R (H4 posterior probability, H4.PP = 1), GCH1 (H4.PP = 0.99), SNCA (H4.PP = 1), TMEM175 (H4.PP = 0.95) show a significant colocalization signal. We performed another colocalization analysis between the pQTLs from the UKB-PPP and GWAS results by Nalls et al.1. We found three loci BST1 (PP.H4 = 0.95), GPNMB (PP.H4 = 0.94) and CTSB (PP.H4 = 0.96) with showing a colocalization support (H4.PP > 0.9). These three genes were also found as significant MR evidences in our analysis. Integrating the MR results with the colocalization signals from the UKB-FinnGen meta-analysis and Nalls et al. signals as well as the pQTLs from the UKB-PPP and Nalls et al., we observe four proteins including GPNMB, BST1, CTSB, and HIP1R to be shared. Among which, GPNMB and BST1 have also been reported in CSF as significant MR signals by Kaiser et al.21. Given that our findings are based on plasma, such an overlap shows the great potential of these two proteins to be investigated as therapeutic targets.Table 3 The identified plasma proteins putatively causal for PD with significant MR and colocalization evidencesTo gain a better insight into the SNP-protein distribution, we performed a specificity analysis of the PD-associated SNPs identified as pQTLs for the proteins profiled in UKB. Initially, we included all genome-wide significant variants from the UKB-FinnGen meta-analysis GWAS, then we narrowed down the analysis to the top significant hits reported in Table 2 (n = 10). Among all the significant variants, the number of proteins per pQTL ranged from 1 to 125 while the number of pQTLs per protein ranged from 1 to 2913 (Figs. 4a, b). The frequencies of cis associations (defined as ±1 kb flanking regions of gene boundaries) as well as trans associations per pQTL were significantly different, where the number of former associations per pQTL ranged from 1 to 16. On the other hand, trans associations ranged from 1 to 109 and 798 pQTLs were only associated in trans (Fig. 4d). Next, we checked the frequency of variants per proteins in both cis as trans associations. We found that 540 proteins had only trans associations and 37 proteins had both types of associations (Fig. 4c).Fig. 4: Overview of blood plasma pQTL architecture in PD.a histogram of the number of PD-associated pQTLs per protein; b histogram of the number of associated proteins per pQTL; c–d number of associations in panels a, b separated by the association type.We investigated the specificity of the SNPs reported in Table 2 to proteins profiled in the UKB-PPP. rs2230288 (missense variant in GBA) shows one significant association in cis with EFNA1. rs34311866 (missense variant in TMEM175) was observed to be significantly associate with IDUA in cis (Beta = -0.05, P = 4.6E-8). Of note, rs2760980 (an intergenic variant close to HLA-DRB1) shows 98 significant associations among which five associations were in cis while the rest were in trans. The cis-associated proteins included AGER, AIF1, DXO, HLA-DRA, and TNXB. We found that rs10847839 (an intronic variant in HIP1R) is in fact associated with HIP1R in the UKB-PPP (Beta=0.09, P = 6.7E-30). rs11158026 (intronic to GCH1) was significantly associated with decreased levels of LGALS3 (Beta = -0.28, -Log10(P) = 316). rs4630591 (intronic to KANSL1) shows 21 significant association among which only one association is in cis with LRRC37A2 (Beta=1.07, -Log10(P) = 2751) while the remaining association show -Log10(P) < 13. Finally, rs6857 (NECTIN2) was found to be highly pleiotropic with 50 significant associations among which only one association in cis with APOE was identified (Beta = -0.85, -Log10(P) = 1573). We can see that these variants are mostly specific to significant cis associations and a few of them are highly pleiotropic with larger p-values in trans associations compared to cis associations.Circulating plasma proteins to study neurological diseases is constrained by the absence of proteins that do not cross the blood-brain barrier (BBB). One alternative tissue for neurological studies is Cerebrospinal Fluid (CSF). Therefore, a systematic comparison between plasma-derived protein measurements versus CSF can potentially shed light on proteins that cannot be captured in plasma either because they cannot make it through the BBB or can make it through at lower concentrations leading to their dilution to noise when preprocessing the data. Leveraging pQTL associations derived from CSF of Parkinson’s disease patients by Kaiser et al.21, we made a systematic evaluation using our identified cis-pQTLs and their associated proteins. Please note that there are systematic differences between the proteomic data used in our study as opposed to the data used by Kaiser et al.21. UKB-PPP proteomics data is generated by the Olink Explore 3072 antibody-based proximity extension assay (PEA), measuring 2,941 protein analytes capturing 2,923 unique proteins. On the other hand, the data used by Kaiser et al.21. is generated by SomaScan aptamer-based proteomics platform37 which captures 4,135 proteins. The two platforms share a total of 1,583 unique proteins. Please note that due to technical differences of data processing on the two platforms as well as the protocols used for characterizing pQTLs, not all the pQTLs might have been shared between the two studies.The CSF-based study reports significant cis pQTLs for 856 Off-rate Modified Aptamers (SOMAmers) corresponding to 744 unique proteins21. Among these proteins, 474 proteins are also profiled in the UKB-PPP data. Comparing the CSF-derived pQTL-associated proteins with our cis-pQTL associated proteins, 151 proteins are shared between the two (Supplementary Data 1) leading to a total of 323 proteins that are found to show significant cis-pQTL associations only in CSF. We conducted a pathway enrichment analysis on significant proteins unique to CSF using KEGG gene sets38. Only one pathway ‘Complement and coagulation cascades’ (FDR = 0.002) passed the significance threshold of FDR < 0.05. We had not identified this pathway to be significant in our analysis based on the pQTL-associated proteins in plasma. The complement system is a tightly regulated innate immune system playing a key role in regulating the normal function and development of the central nervous system (CNS)39. According to Fatoba et al.39., accumulating evidence from human postmortem studies demonstrate that neuronal and glia cells are capable of expressing a majority of the complement molecules in CNS. This pathway has been reported in several CNS-related studies such as multiple system atrophy40, brain proteome profiling of PD in humans and mice41,42,43, and psychosis44.In addition to evaluating the pQTL-associated proteins in both datasets, we also compared the cis-pQTLs that were identified in both studies. Two pQTLs were shared on both studies including rs429358 and rs557011. In CSF, the missense variant rs429358 was significantly associated with three proteins APOC2 (Beta=0.89, P = 3.56E-39), APOE (Beta=0.82, P = 2.49E-45) and CEACAM19 (Beta=0.46, P = 4.05E-17). However, this variant was only significantly associated with APOE in plasma with a negative direction of effect (Beta = -1.01, -log10(P) = 2069.74). rs557011 was associated with HLA-DQA2 in CSF (Beta=0.53, P = 1.72E-42) while being significantly associated with HLA-DRA in plasma (Beta=0.37, -log10(P) = 455.6). The main difference observed is the opposite association of rs429358 with APOE in plasma versus CSF.Our observation indicates the intrinsic limitations of measuring circulating plasma proteins to study brain diseases and need for leveraging CNS-specific tissues to shed light on genome-proteome aberrations in PD that may not be captured in plasma.Network modeling reveals PD-associated protein complexesThe complex etiology of PD and its polygenic nature requires employing a systems-level approach for joint modeling of distinct data modalities to identify the interactions between the genetic variants with encoded proteins as well as protein-protein interactions that might be disrupted as the outcome of genetic aberrations. Network modeling is among the most efficient tools for modeling complex interactions and analyzing their overall contribution to the disease onset and progression19,45. Using the gene-based scores from GWAS meta-analysis of UKB-FinnGen, the identified proteins, and existing protein-protein interaction (PPI) networks, we performed hypothesis-free and hypothesis-driven network analyses to identify protein complexes associated with PD and to characterize potential mediators that indirectly contribute to the disease. In the following, we will discuss the details of each experiment.To conduct a hypothesis-free network analysis with the goal of identifying potential targets and mediators through the integration of PD-associated genetic signals and PPI networks, we used the GWAS meta-analysis results of UKB-FinnGen and calculated a gene-based score based on the p-values of all the SNPs localizing within the gene boundaries using MAGMA32. We employed the PPI network from PICKLE29, which includes 217,963 interactions containing 16,420 proteins. Then, we used gene-based scores to assign weights to each of the corresponding nodes in a PPI network (Methods). Using the generated weighted PPI network, we identified PD-associated modules using dmGWAS46. Here, a ‘module’ is defined as a subgraph of the larger PPI network with a local maximum proportion of low p-value genes. Upon completing the analysis, each module will have a score based on the gene-based p-values of its member nodes. Overall, 13,897 modules with a varying number of nodes were identified. We sorted the modules based on their scores, focusing on the top 1% significantly PD-associated modules. In total, the top-identified modules contained 175 proteins (Supplementary Data 4). Next, we compared the list of 577 proteins, that were associated with the identified pQTLs, with the identified the proteins in the top 1% of modules. We found 14 proteins (P = 0.003) to be shared between the two groups, including: SNCA, LGALS3, CSF1R, STX4, APOC1, APOE, APOA2, SMPD3, PXN, PAFAH1B3, LDLR, HSPB1, BRK1, and CSNK2A1 (Fig. 5a). Please note that the red nodes shown in Fig. 5a include all the genome-wide significant loci irrespective of being the top locus in the meta-analysis GWAS. Notably, two established familial PD risk loci appeared among the top disease modules only through this PPI analysis – PINK1 and PARK7.Fig. 5: Network analysis results.a hypothesis-free analysis: aggregated view of the top 1% of the identified PD-associated protein modules. Red nodes represent GWAS signals, green nodes represent pQTL-associated proteins in PD, dark blue nodes are pQTL-associated proteins which are PD-associated GWAS loci, and the size of the nodes is proportional to the gene-based scores generated by MAGMA; b hypothesis-driven analysis: aggregated view of the identified loci from UKB-FinnGen meta-analysis GWAS (shown in red) and their first-degree neighbors in the PPI network. Green nodes represent pQTL-associated proteins in PD.Genetic correlations between PD, Alzheimer’s disease, and Lewy body dementia (LBD) have been investigated in the literature47,48 and a few loci have been reported to be shared between these diseases such as TMEM175, MAPT, KANSL1. We sought to identify potential overlaps between the 175 proteins in the top 1% module of the hypothesis-free network analysis and risk genes identified in Alzheimer’s Disease. Using GWAS Catalog49, we found 23 overlapping genes, including APOE, MAPT, PLEKHM1, KANSL1, NSF, KAT8, BCKDK, APOC1, TOMM40, LRIG1, ICA1L, USP8, UBXN11, FAM210B, AVL9, DISC1, PRKAA1, BIN1, TEAD1, LDLR, CSNK2A1, L3MBTL3, and RAC1. Moreover, we investigated shared genes among the top identified modules with LBD50 where we observed a significant overlap at P = 1.8E-63. These findings can be further explored for examining the shared molecular signatures among neurodegenerative disorders.We next conducted a hypothesis-driven network analysis, mapping the UKB-Finngen PD meta-analysis GWAS significant loci onto the human PPI network, and then extracting their first-degree neighbors in the network (Fig. 5b). Among the first-degree neighbors, pQTL-associated proteins identified through proteogenomic analysis are shown in green. In total, 53 proteins were found to interact with PD risk nodes. We ran a set of permutation tests to check if the overlap between the pQTL-associated proteins and the direct neighbors of PD-associated risk loci was not by chance. We randomly selected proteins from the entire PPI network with a size equal to the network presented in Fig. 5b and tested their overlap with the pQTL-associated proteins for 10,000 iterations. In each iteration, we tested the overlap using Fisher’s exact test and created an empirical distribution of the p-values. An empirical p-value of 0.014 was observed, denoting significant enrichment of pQTL-associated proteins within the direct neighbors of PD GWAS loci in the human PPI network. The extracted network in this section was then compared with the top 1% modules in the hypothesis-free network analysis. Considering all the nodes in both networks, we observed a significant overlap between the two (Fisher’s exact test P = 7.6E-56). Such observation indicates the importance of considering not only the PD-associated loci from GWAS, but also their interacting proteins whose dysregulation may contribute to the pathways that contribute to the disease. One of the limitations of this analysis is the constrained number of proteins profiled in the UKB-PPP. Therefore, hypothesis-driven analysis may not reflect the full spectrum of the human proteome and as a result, this might have led to overlooking some parts of the human interactome that might be involved in the disease.Transcriptomic and Epigenomic data analysisIn the third phase of our analytical pipeline, we conducted epigenomic and transcriptomic analyses of genes identified from the first two phases. We first studied how GWAS-significant PD risk variants potentially disrupt the expression of linked genes through physical chromatin loops. Second, we investigated how the expression of the identified GWAS loci, as well as their associated proteins in the brain, may be disrupted in PD.We investigated how the identified PD variants impact disease risk via potential influences on cell type-specific enhancer-promoter interactions using PLAC-seq data by Nott et al.51. First, we mapped the identified GWAS SNPs (Fig. 2b) onto the developed cell type-specific interactome, examining whether any variants localized to regulatory regions of their corresponding genes or impacted the regulatory regions of distal genes through physical interactions with chromatin loops.We found 6 PD variants to be linked with distal genes in cell types essential for central nervous system function (Supplementary Data 5). Two variants were found to exert regulatory effects only in neuronal cells, including rs356219 (intronic variant in SNCA) and rs113100008 (intronic variant in CRHR1). One variant, rs7542186 (intronic to LMNA), was found to act only in oligodendrocytes. Two variants impacted both neuronal cells and oligodendrocytes, including rs10847839 (intronic variant in HIP1R) and rs8067056 (intronic to STH). Finally, only rs429358 (one of the two APOE e4-defining missense variants) was found to be active in microglia. We checked the localization patterns of the linked genes to the PD variants in each cell type. For all the variants, the linked genes localized in the same region. Since the proteomic data in the UKB-PPP is plasma-derived, we sought to investigate if any of the identified cell-specific variants in Supplementary Data 5 localize in open chromatin regions in blood-derived data, to investigate if these variants are brain-specific. We used the assay for transposase-accessible chromatin using sequencing (ATAC-seq) data52 processed by ATACdb53 generated from healthy individuals. None of the brain-specific variants were found to be localizing in open chromatin regions in blood implicating that such variants have cell-specific regulatory roles in brain.To complement the analysis, we conducted a fine-mapping (Methods) on the UKB-FinnGen meta-analysis results aimed at identifying the putative causal variants. Then we integrated them with the data by Nott et al. We consider a variant causal if the fine-mapping posterior probability (PP) > 0.9. We found two variants rs356219 (SNCA) and rs10847839 (HIP1R) (Supplementary Data 5) to be causal at PP = 1 which were both identified to be acting in neuronal cells. These variants have also been reported by Schilder and Raj54 to be causal in PD. The remaining variants reported in Supplementary Data 5 did not have a PP > 0.9. Moreover, no other significant variants with PP > 0.9 were found to fall in the chromatin interaction loops provided by Nott et al.51.We examined the cell-specificity of the identified PD risk variants as well as their dysregulations in PD patients compared to normal controls. First, we leveraged single-nucleus droplet-based sequencing data from Lake et al.55 to evaluate if the identified PD risk genes are markers of specific cell-types in distinct brain regions. These data cover >60,000 single nuclei from the human adult visual cortex, frontal cortex, and cerebellum. Using the summary statistics provided by Lake et al.55, we extracted significantly upregulated genes for each cluster of cells and used them as the marker for that cell type. Among all the PD GWAS risk genes, we found three genes to act as specific markers of certain cell types, including SNCA and LRRK2 in excitatory neurons and KANSL1 in Purkinje cells.To further investigate cell specificity in more disease-relevant brain tissues, we leveraged another set of single-cell RNA-seq (scRNA-seq) data from Kamath et al.56. These data56 contain >22,000 dopamine neurons extracted from substantia nigra pars compacta, where ten populations of dopamine neurons were identified. We evaluated how the identified PD loci in Tables 1 and 2 were differentially expressed (DE) between PD, LBD, and healthy control individuals (Fig. 6a). Among the investigated genes, 8 genes were DE in PD compared to healthy controls. Moreover, we observed different dysregulation of several genes in PD versus LBD including SNCA, STK39, HIP1R, GCH1, KANSL1, which were overexpressed in >75% of cells in PD. However, three other genes GBA, MMRN1, and CRHR1 were expressed in less than 50% of cells. Next, we focused on the shared proteins between the hypothesis-free and hypothesis-driven network analyses results and examined how their encoding genes are expressed in PD, LBD, and healthy control individuals in dopamine neurons, microglia, and astrocytes (Figs. 6b–d). These proteins included: SNCA, APOE, CSNK2A1, SMPD3, STX4, APOA2, PAFAH1B3, LDLR, HSPB1, BRK1, and LGALS3. Among the genes encoding these proteins, HSPB1 was the only gene to be upregulated in >75% of the cells in dopamine neurons, microglia, and astrocytes. APOE was upregulated in >75% of cells in microglia within PD and LBD cohorts while it was downregulated in astrocytes in the PD cohort. SNCA was significantly upregulated in dopamine neurons and microglia in >75% of cells while it showed upregulation only in ~38% of cells in astrocytes. SMPD3 was significantly upregulated in dopamine neurons in >75% of cells while showing the same expression patterns in microglia and astrocytes at just a small fraction of cells. STX4, APOA2, and PAFAH1B3 showed a stable number of expressing cells across the three cell types with a varying degree of expression levels. BRK1 shows a larger number of expressing cells in dopamine neurons while less expressing cells in microglia and astrocytes. We found that BRK1 is significantly downregulated in PD patients compared to controls while being significantly upregulated in LBD versus controls. Finally, we found LGALS3 to be significantly upregulated in PD and LBD in ~40% of cells in microglia while showing downregulation in LBD in astrocytes.Fig. 6: Expression profiles of PD-associated genes among PD, LBD, and healthy individuals.a Expression of PD-associated genes identified through UKB GWAS and UKB-FinnGen meta-analysis GWAS; b–d Expression of the gene encoding the identified proteins from the network analyzes in dopamine (DA) neurons, microglia, and astrocytes, respectively. Source data are provided as a Source Data file.Converging lines of evidence suggests 9 proteins as potential mediators of PDWe sought to identify ‘indirect’ mediators of PD risk–i.e., proteins that may not harbor a PD susceptibility risk variant but may contribute to disease-relevant mechanisms such as neuroinflammation. Based on the proteogenomic analysis and different network modeling results, we defined three sets of targets as follows: (a) pQTL-associated proteins in blood plasma derived from UKB-PPP; (b) proteins that were identified among the top 1% of PD associated modules from the hypothesis-free network analysis whose abundance was also associated with PD pQTLs; (c) PD pQTL-associated proteins which directly interact with the PD meta-analysis GWAS loci in the hypothesis-driven network analysis. Then, we mined the literature to collect additional independent causal inference evidence on each protein. Overall, 9 proteins with converging lines of evidence from the three target sets were identified, including: LGALS3, CSNK2A1, SMPD3, STX4, APOA2, PAFAH1B3, LDLR, HSPB1, BRK1. We screened these proteins for previous peer-reviewed studies linking them to PD and/or other neurodegenerative diseases.We identified a number of studies reporting links between PD and LDLR57,58 and CSNK2A159. Two prior studies identified associations between PD and STX4. In the first, Tewhey et al.60 employed multiplexed reporter assays to investigate the downstream effects of GWAS risk variants across several traits including PD, reporting an association between the PD risk variant, rs11865038, and expression of STX4. In the second study, STX4 was identified as a common risk gene for AD and PD through a transcriptome-wide association study (TWAS)61. We further identified three prior studies on the impact of SMPD3 on neurodegeneration and cognitive impairments62,63,64. Several lines of evidence were found about LGALS3. For the remaining nominated proteins, we did not observe published studies to demonstrate their links with neurodegenerative disorders.We found a series of research studies highlighting the role of LGALS3 in mediating microglial-mediated neuroinflammation in AD65,66, cognitive function, and neurodegeneration67,68,69, and PD70,71. Two recent causal inference studies, employing Mendelian randomization (MR), have reported a significant causal association between LGALS3 and PD21,72. In the next section, we further explore the potential impact of LGALS3 on the development of PD.Potential role of galectin-3 (LGALS3) in PDAs an evolutionarily conserved protein, LGALS3 drives inflammation73. The role of LGALS3 in the immune system has been extensively investigated, but its role in the central nervous system is less explored. LGALS3 is observed to be upregulated in several central nervous system diseases associated with inflammation, including AD, hypoxia, and stroke65,74. Such an observation is in line with our finding in analyzing the PD scRNA-seq data (Fig. 6). LGALS3 activates microglia and inflammation in diseases such as stroke75, Huntington’s disease76, and multiple sclerosis (MS) and it is found that it can induce microglial activation only when coupled with tissue damage such as infection or neurodegeneration74. Our initial hypothesis about LGALS3 was based on the PD protective variant rs11158026 (beta = −0.1) identified through the meta-analysis of UKB-FinnGen GWAS. Given that rs11158026 was a pQTL associated with downregulation of LGALS3, we hypothesized that this variant is exerting its effect through downregulation of LGALS3. So far, LGALS3 has been explored as a target by a number of teams and companies for diseases such as heart disease, fibrosis, and AD where a few drug compounds have made it to the phase II and III clinical trials77,78,79. However, to our knowledge, no clinical efforts are in place for targeting LGALS3 in brain or inflammatory diseases. Therefore, we believe further investigations to study how inhibition of LGALS3 can ameliorate the footprints of PD is warranted.LGALS3 is a member of family of galectins which are a set of proteins that share a carbohydrate recognition domain motif which interacts with β-galactoside glycan80. Galectins play a role in a wide array of cellular functions including signaling, inflammation, autophagy, and immune response80. We queried the GWAS significant SNPs in the proteogenomic analyzes to identify if LGALS3 abundance is associated with any the identified PD-associated variants. We had identified rs11158026-T as a significant SNP with a protective effect (beta = −0.1, P = 3.4E-8) on PD in UKB-FinnGen meta-analysis GWAS (Table 2). We hypothesized that the identified protective association can potentially be explained through the results of proteogenomic analysis and found that rs11158026-T in GCH1 is a pQTL associated with reduced LGALS3 expression (beta = -0.286, -log10(P) = 316). Therefore, we hypothesized that rs11158026 exerts protective effects against PD by downregulating LGALS3 (Supplementary Fig. 2). A recent study by Burbidge et al.70 identified that LGALS3 mediates unconventional secretion of SNCA/alpha-synuclein in response to lysosomal membrane damage in human midbrain dopamine neurons. Based on this finding, protective association of rs11158026-T with PD, and its association with reduced expression of LGALS3, we can hypothesize that rs11158026 has a protective role against PD through downregulation of LGALS3. Through MR analysis performed in the Proteogenomic analysis section, we found LGALS3 to significantly associate with PD (beta=0.22, FDR = 0.001) indicating that its overexpression is putatively causal to increased risk of PD. In our MR analysis, rs11158026 was the instrumental variable used in the Wald Ratio test. Testing for potential pleiotropic associations of this SNP, we found no significant pleiotropic association for this variant (Horizontal pleiotropy P-value = 0.94). To further validate this hypothesis, we used two MR analysis results generated from human brain and CSF. In the first analysis, Storm et al.72 performed an MR using eQTL data from human brain tissue to identify genes with potential causal associations with PD. They identified LGALS3 (OR = 2.11, FDR = 0.004) as a potential causal gene whose overexpression is associated with increased risk of PD. In a second independent study, Kaiser et al.21 performed an MR analysis using the existing GWAS summary statistics on PD along with pQTL data from human CSF. They have also identified a causal association between LGALS3 and PD (beta=0.07, FDR = 7.1E-4), indicating that overexpression of LGALS3 is significantly associated with increased risk of PD. In fact, the direction of effect between our findings and observations and the developed hypothesis, match both independent MR studies. Convergence of all these independent lines of evidence had led us to nominate LGALS3 as a potential mediator of PD which can indirectly contribute to the PD pathogenesis.We investigated the abundance of LGALS3 in the UKB-PPP data and analyzed its genotype-specific abundance among PD patients. To gain a deeper insight into the genotype-specific abundance of LGALS3 in pre-diagnosis versus diagnosed PD cases, we divided the PD cohort into ‘prevalent’ and ‘incident’ cases. Prevalent cases are patients who have been diagnosed with PD when donating blood samples, while incident cases are individuals who have not been diagnosed with PD when donating blood and were diagnosed later, sometimes years later. Three sets of tests were performed. Initially, we used all the incident and prevalent cases (Fig. 7a) and ran an ANCOVA (analysis of covariance) to evaluate the abundance of LGALS3 between patients with different genotypes at rs11158026. Then, we ran a similar experiment on prevalent and incident cases only (Figs. 7b, c). The combined set of cohorts shows a significant decrease in LGALS3 abundance in homozygous carriers of rs11158026-T compared to the heterozygous carriers. Both were also showing a lower abundance compared to non-carrier patients. Focusing on prevalent cases, we found a significant downregulation of LGALS3 in TT carriers versus non-carriers as well as a marginally significant downregulation in CT compared to CC. However, no significant difference between homozygous and heterozygous carriers was observed (Fig. 7b). Finally, considering only incident cases, we found a significant downregulation in the whole three sets of comparisons (Fig. 7c). Moreover, the number of non-carriers (n = 363) is larger than the carriers among incident cases (CT, n = 307; TT, n = 83). So we compared the number of carriers and non-carriers between the incident cases and the control cohort in the UK Biobank. We found heterozygous carriers of this variant had a marginally significantly lower proportion of PD incident cases compared to non-carriers (Fisher’s exact test P = 0.04, OR = 0.8). Similarly, homozygous carriers who later developed PD were significantly less than non-carriers (Fisher’s exact test P = 0.02, OR = 0.75). Genotype-specific abundance of LGALS3 among the incident cases, prevalent cases, and control individuals were conducted (Fig. 7d). In each genotype, we did not observe a significant difference between the three sets of individuals.Fig. 7: Genotype-specific abundance of LGALS3 among carriers of the protective SNP rs11158026 among.a incident+prevalent PD cases; b prevalent PD cases only; c incident PD cases only; d cases and controls. NPX: Normalized Protein Expression. In (d) the p-values indicated at the bottom of the figure represent the differences between the NPX values once combining incident, prevalent and controls in each genotype. All data panels are represented as median values plus the first and third quartiles. Two-sided ANCOVA test was performed to evaluate the genotype-specific abundance of LGALS3 adjusted for age and sex. Source data are provided as a Source Data file.
