Immune system disruptions implicated in whole blood epigenome-wide association study of depression among Parkinson's disease patients

Although Parkinson's Disease (PD) is typically described in terms of motor symptoms, depression is a common feature. We explored whether depression influences blood-based genome-wide DNA methylation (DNAm) in 692 subjects from a population-based PD case-control study, using both a history of clinically diagnosed depression and current depressive symptoms measured by the geriatric depression scale (GDS). While PD patients in general had more immune activation and more accelerated epigenetic immune system aging than controls, the patients experiencing current depressive symptoms (GDS≥5) showed even higher levels of both markers than patients without current depressive symptoms (GDS<5). For PD patients with a history of clinical depression compared to those without, we found no differences in immune cell composition. However, a history of clinical depression among patients was associated with differentially methylated CpGs. Epigenome-wide association analysis (EWAS) revealed 35 CpGs associated at an FDR≤0.05 (569 CpGs at FDR≤0.10, 1718 CpGs at FDR≤0.15). Gene set enrichment analysis implicated immune system pathways, including immunoregulatory interactions between lymphoid and non-lymphoid cells (p-adj = 0.003) and cytokine-cytokine receptor interaction (p-adj = 0.004). Based on functional genomics, 25 (71%) of the FDR≤0.05 CpGs were associated with genetic variation at 45 different methylation quantitative trait loci (meQTL). Twenty-six of the meQTLs were also expression QTLs (eQTLs) associated with the abundance of 53 transcripts in blood and 22 transcripts in brain (substantia nigra, putamen basal ganglia, or frontal cortex). Notably, cg15199181 was strongly related to rs823114 (SNP-CpG p-value = 3.27E-310), a SNP identified in a PD meta-GWAS and related to differential expression of PM20D1, RAB29, SLC41A1, and NUCKS1. The entire set of genes detected through functional genomics was most strongly overrepresented for interferon-gamma-mediated signaling pathway (enrichment ratio = 18.8, FDR = 4.4e-03) and T cell receptor signaling pathway (enrichment ratio = 13.2, FDR = 4.4e-03). Overall, the current study provides evidence of immune system involvement in depression among Parkinson's patients.


Introduction
Parkinson's disease (PD) is a progressive neurodegenerative disorder characterized by gradual loss of dopaminergic neurons in the substantia nigra region of the brain. While PD is typically described in terms of motor symptoms, neuropsychiatric symptoms and especially depression are common, both preceding and subsequent to motor symptom onset (Aarsland et al., 2012). Depression in PD is strongly related to diminished health-related quality of life, reduced function, and cognitive decline (Aarsland et al., 2012;Lieberman, 2006). Although often underdiagnosed, clinically relevant depressive symptoms occur in an estimated 35% of PD patients and epidemiologic evidence indicates a strong association between a history of major depression and subsequent development of PD(1).
The question remains whether depression is a prodromal symptom of PD or a risk factor for PD incidence. However, whichever is the case, there is mounting evidence for biologic underpinnings of depression in PD, including structural changes in the brain of depressed PD patients relative to non-depressed patients (e.g. loss of striatal dopamine transporter availability (Weintraub et al., 2005); loss of white matter within the cortico-limbic network (Kostić et al., 2010)), alterations in neurotransmitter systems (e.g. dopaminergic and serotonergic circuits (Grosch et al., 2016;Lu et al., 2019;Mayberg et al., 1990)), and an involvement of inflammatory and neurotrophic factors (e.g. extended stress-induced activation of the brain via cytokines and glucocorticoids (Pace and Miller, 2009); downregulation of the MAPK-MEK pathway (Wang and Mao, 2019)).
Depression likely arises due to a combination of environmental and genetic factors (Aarsland et al., 2012). Currently, we have a limited understanding of etiologic mechanisms that contribute to depression in PD. Here we are concentrating on DNA methylation (DNAm) as it may offer a readout for biologic pathways that are impacted. DNAm also regulates gene expression/repression (Moore et al., 2013) and it also reflects both environmental and genetic determinants (Law and Holland, 2019). Specifically, DNAm alterations may reflect key mechanisms through which exposures interact with genetic predisposition to determine an individual's susceptibility. Furthermore, depression has also long been linked to immune dysregulation (Maes, 1995). For instance, PD, depression, and related risk factors, such as psychosocial stressors, can induce proinflammatory cytokines and have been associated with sustained epigenetic changes (Slavich and Irwin, 2014;Reale et al., 2009;Cunliffe, 2016;Zannas, 2019). Thus, blood may be an informative tissue for assessing depression related DNAm changes.
To investigate blood-based DNA methylation related to depression in PD, we performed an epigenome-wide association study (EWAS) of both having a history of clinical depression and current depressive symptoms in 692 participants of a population-based case control study of PD with longitudinal data for patients. We further used publicly available data (GTEx and BIOS) to map the associated CpGs to different biologic layers, using methylation and expression quantitative trait loci (QTL) to assess networks of CpGs, SNPs, and transcripts related to depression in PD.

Results
Our analysis of blood-based DNAm draws from participants recruited as part of the Parkinson's disease, Environment, and Genes (PEG) study (Ritz et al., 2016). Data are available on Gene Expression Omnibus (GEO), accession numbers GSE72774 and GSE72776(18,19). PEG is a population-based study of Parkinson's disease, which enrolled patients and controls from California's Central Valley (2001-2007 & 2010-2016). Primary analysis was restricted to 465 PD patients of European ancestry (based on AIMs ancestry markers) for whom we had DNAm data available and 227 controls of European ancestry. Given the limited sample size, we used the control population only to assess whether depression associations seen among the PD patients were also seen among the controls (i.e., related to depression in general, or related to depression in PD specifically). Characteristics of the study population can be found in Table 1.
We assessed two indicators related to depression, first, having a history of clinical depression based on self-reported physician diagnosis and, second, having current depressive symptoms measured by the geriatric depression scale (GDS) at the time of blood draw. We dichotomized the GDS into groups with no or low depressive symptoms (GDS<5) and with depressive symptoms (GDS≥5). We selected this GDS cut-point as we have previously validated the indicator (GDS≥5) as having high specificity and sensitivity in distinguishing minor and major depression in a subset of our patient population (Thompson et al., 2011). This validity has also been shown in other studies of PD(21). An overview of the analysis plan is shown in Supplemental Fig. 1.
Overall, 130 PD patients (28%) reported a history of clinical depression and 130 PD patients (28%) had current depressive symptoms at baseline. The two depression indicators were moderately correlated among the PD patients (Pearson's R = 0.32, 95% CI = 0.25, 0.39), with 70 of the patients having both a history of clinical depression and experiencing current depressive symptoms at the time of blood draw.
However, when examining patients by a history of clinical depression instead of based on current depressive symptoms (GDS≥5), we found these differences were limited to current symptoms. Immune cell composition was similar among PD patients with and without a history of clinical depression (Supplemental Fig. 2). Furthermore, the observed differences were also specific to PD patients. Non-PD controls with higher levels of depressive symptoms (GDS≥5) and a history of clinical depression had a similar cell composition and immune system age acceleration as controls without depressive symptoms or a history of clinical depression (Supplemental Fig. 3).

Depression epigenome-wide association study
In our primary analysis, we related the two depression indicators to genome-wide DNAm levels among PD patients, adjusting for cell composition (proportion of neutrophils and CD4Ts), age, sex, smoking, PD duration at baseline, European fractional ancestry, and study wave. Without controlling for cell composition, the GDS≥5 indicator for current depressive symptoms was related to 129 CpGs at an FDR<0.05 and 675 CpGs at an FDR<0.10 among PD patients. However, after controlling for cell composition, the current depressive symptoms indicator was not related to any CpG levels in site-by-site analysis (p > 7.1e-6; FDR>0.99).

Fig. 1. Epigenetic Immune System Markers and Current Depressive Symptoms in PD.
Mean comparisons of A) Houseman DNAm estimated leukocyte proportions, B) Neutrophil-to-lymphocyte ratio (based on DNAm leukocyte proportions), and C) epigenetic immune system age acceleration (EEAA) across three groups: controls (no depressive symptoms), PD patients without depressive symptoms (GDS<5), and PD patients with depressive symptoms (GDS>5).
On the other hand, having a history of clinical depression prior to PD onset was associated with several differentially methylated positions (DMP). An overview of the results is displayed in a Manhattan plot (Supplemental Fig. 4). Overall, having a history of clinical depression was associated with 35 CpGs at an FDR<0.05, 569 CpGs at an FDR<0.10, 1718 CpGs at an FDR<0.15, and 5270 CpGs at an FDR<0.25 (Supplemental Table 1). The top CpGs based on FDR are detailed in Table 2. The most significantly associated CpG was cg18774195 (FDR = 0.007) in the 5 ′ UTR region of the SLC7A14 gene. Other notable CpGs included cg23426156 and cg11042505, linked to schizophrenia in EWAS(25,26), 32 CpGs linked to aging, and 6 CpGs linked to smoking. These CpGs are listed in the supplemental materials and were determined via query of the MRC-IEU catalog of epigenome-wide association studies (Battram et al., 2022). In sensitivity analysis, we assessed whether the CpGs associated with clinical depression among PD patients were also associated with a history of clinical depression among controls (Supplemental Table 1). From the 35 CpGs associated at FDR<0.05 among the PD patients, three CpGs, cg22801913, cg21276379, and cg25290938 were also associated with depression among controls at p < 0.05. The three CpGs were located in body of gene C11orf49, the transcription start site of EPAS1, and the transcription start site of ARHGAP22, respectively, and all three were also associated with aging in a prior EWAS(27). Otherwise, of the 5270 CpGs associated at FDR<0.25 among the patients, 242 in total were also associated with depression among controls (p < 0.05), indicating that among our population most of the depression CpGs were specific to depression among PD.
Among the patients, having a history of clinical depression prior to PD onset was also associated with differentially methylation regions (DMR) identified via Bumphunter in the ChAMP R package. Supplemental Table 2 details the DMR results. The most significant DMR was Chr1:205818668-205819609 (p = 2.0e-5, FWER = 0.007), a 941 base pair region with 9 CpGs. This region is part of the PM20D1 gene. PM20D1 has been associated with both PD and Alzheimer's via meta-GWAS (top PD SNP in gene has GWAS meta p-value = 6.45e-08 (28,29)) and is near NUCKS1, another PD gene identified in the meta-GWAS (top SNP in gene meta p-value = 1.96e-16(28,29)), which we identified as important in the quantitative trait loci (QTL) analysis described below. Chr6:32145470-32146232 DMR was also suggestively associated with having a history of clinical depression prior to PD onset after multiple-testing correction (p = 2.7e-4, FWER = 0.09). This is a 762 base pair region with 28 CpGs in the AGPAT1 gene of the major histone compatibility (MHC) class III region, which per NCBI is involved in phospholipid metabolism and metabolism gene pathways. Differential expression of the AGPAT1 protein has previously been linked to major depression (Scifo et al., 2018).
As the PD patient population used here is unique in that they all live in a highly agricultural region on California with high levels of pesticide use, we also conducted a sensitivity analysis controlling for pesticide exposure which we have previously linked to methylation levels (Paul et al., 2018a;Furlong et al., 2020). The results were very similar with a correlation coefficient between model predicted betas at R = 0.9 (p < 2.2e-16; Supplemental Fig. 5). CpG associations for different sensitivity models can be found in Supplemental Table 1.

Blood-brain CpG level correlation
Given the neurologic pathogenesis of PD and depression as well as the enrichment of neural function pathways, we next assessed the correlation of methylation levels in blood and brain for the depression EWAS CpGs using the Iowa Methylation Array Graphing Experiment for Comparison of Peripheral Tissue and Grey matter (IMAGE-CpG) data (n = 37; two studies) (Braun et al., 2019). Limiting this analysis to the 35 CpGs at FDR≤0.05, multiple CpGs showed both relatively high positive and negative correlations for methylation in blood and brain ( Fig. 2A).

Methylation and expression quantitative trait loci
Finally, we assessed how the EWAS-associated CpGs were related to other omic layers. First, using BIOS and our study data, we determined methylation quantitative trait loci (meQTLs, e.g. SNPs which are significantly associated with CpG levels). Second, using GTEx, we determined whether the meQTLs were also related to transcript abundance measured in both whole blood and brain. BIOS is a public database of meQTLs based on whole blood methylation, determined in 3,841 samples from five Dutch biobanks (generally age>45 years, including one longevity study with n =~2600 participants >89 years of age) (Bonder et al., 2017). GTEx provides characterization of genetic associations with gene expression based on 838 individuals [mean age 53.4 (range 21-70)], 52 tissues, and two cell lines (Aguet et al., 2017;Lonsdale et al., 2013).
We determined that amongst the 35 CpGs at FDR≤0.05, 25 were associated with 45 different meQTLs (Supplemental Table 3 and  Table 4). This indicates that 45 SNPs (i.e. meQTLs) predicted the levels of 25 of the depression EWAS CpGs. Using GTEx, we determined 26 of these meQTLs were also expression QTLs (eQTLs), as the SNP was also associated with blood-based transcript abundance of 53 different gene transcripts (Table 4). For instance, methylation levels of cg26297819 in the transcription start site (TSS200) of the TEF gene was significantly associated with rs202637, a nearby intergenic SNP, based on both BIOS and PEG data (SNP-CpG p-value of 7.9e-7). rs202637 was also significantly related to the transcript abundance of TEF measured in blood (SNP-transcript p-value = 7.98e-28), as well as three other transcripts (MEI1, DESI1, and PMM1; Table 4). Interestingly, rs202637 is associated with PD (meta-GWAS p < 0.05 (Lill et al., 2012;Nalls et al., 2014)) and TEF polymorphisms have been linked with depression and depression in PD specifically (Kripke et al., 2009;Hua et al., 2012). cg15199181, an intergenic CpG, was associated with four meQTLs, most strongly rs823114 (SNP-CpG p-value = 3.27E-310). Not only was rs823114 also picked up in the PD meta-GWAS (p = 1.78E-13 (Lill et al., 2012;Nalls et al., 2014)), but this SNP was also related to blood-based transcript levels of four different genes, PM20D1, RAB29, SLC41A1, and NUCKS1 based on GTEx. One final notable CpG was cg21769117 in the CLIC1 gene transcription start site (TSS1500), which is located in the MHC class III region. This CpG was significantly associated with three meQTLs, SNPs which were also eQTLs associated with the abundance of 18 different transcripts from genes in the region, many immune system related, including HLA genes and complement system protein genes. cg21769117 and its associated SNP and transcript network are mapped onto the chromosome in Fig. 2B. Supplemental Fig. 6 shows all genes in the dense chromosomal region where cg21769117 and the me/eQTLs are located. This figure demonstrates that the me/eQTLs were not associated with transcript abundance of every gene in the region surrounding cg21769117, but instead selectively associated with the 18 specific transcripts.
To synthesize the network of associated CpG, SNP, and transcripts, we further took the entire set of implicated genes (79 genes with an associated CpG from the EWAS, SNP via meQTL, or transcript associated with the me/eQTLs in blood) and performed gene set overrepresentation analysis to test for overrepresentation of biologic processes. The overrepresented (FDR<0.05) gene ontology (GO) biologic processes are shown in Supplemental Table 4 and Fig. 2C. Given the large number of transcripts from the MHC class III region associated with the three me/ eQTLs linked to cg21769117, immune system processes were strongly overrepresented, most significantly: interferon-gamma-mediated signaling pathway (enrichment ratio = 18.8, FDR = 4.4e-03); T cell receptor signaling pathway (enrichment ratio = 13.2, FDR = 4.4e-03); antigen receptor-mediated signaling pathway (enrichment ratio = 10.9, FDR = 0.01). Many of the associated biologic processes are related      NA: Not a applicable (CpG not associated to meQTL, or meQTL is not associated with transcript adundance and not also an eQTL. PDGene (meta-GWAS): ƚ meta p < 0.05; ƚƚ meta p = 1.78E-13. a Ensembl_id when HGNC unavailable.
processes, shown in a graph of enriched GO terms that displays the connections between the GO terms (Fig. 2D). Interestingly, of 45 me/eQTLs that were associated with abundance of 53 different transcripts in blood, 14 were also associated with the abundance of 22 different transcripts in one of three brain regions of interest: the substantia nigra, putamen basal ganglia, and frontal cortex (Table 5). Notably, fewer transcripts in the MHC class III region were related to the eQTLs (8 in the brain versus 18 in blood). There was also differential abundance of three transcripts, CYP21A1P, RAB29, and CPNE1, in the substantia nigra specifically related to three eQTLs. The PD GWAS SNP rs823114, which was associated with cg15199181, was also linked to expression of RAB29 in all three regions (substantia nigra, putamen basal ganglia, and frontal cortex) and PM20D1 in the frontal cortex.

Discussion
Depression is highly prevalent in Parkinson's disease (PD) (Kasten et al., 2010). Yet, the etiology of this important non-motor feature, while likely multifactorial, is not well understood (Simuni and Sethi, 2008). To assess biologic disruptions associated with depression in PD, we used a population-based study of patients early in disease course and employed blood-based methylation to relate leukocyte composition and genome-wide white blood cell methylation levels to two indicators of depression: having a history of clinical depression and current depressive symptoms. Although the two depression measures were moderately correlated (R = 0.32), indicating some patients had both a history of clinical depression and current depressive symptoms, each indicator was associated with different methylation patterns. Nevertheless, both implicated immune system function involvement.
We first characterized PD patient immune cell profiles via immunomethylomics and methylation profiling (Kelsey and Wiencke, 2018). Levels of several leukocyte indicators suggested more immune activation in patients with current depressive symptoms than in either PD patients or controls without current depressive symptoms. This was perhaps most comprehensively shown by the neutrophil-to-lymphocyte ratio (NLR). The NLR, which is simply a ratio of lymphocyte to neutrophil counts, here measured with epigenetic surrogates as leukocyte markers, is often assessed as a measure of systemic inflammation and immune activation (Song et al., 2021). It reflects the balance between inflammation (acute or chronic), as measured by neutrophil levels, and adaptive immunity, measured by lymphocyte levels (Song et al., 2021). Interestingly, even our controls were showing an average NLR of 2.1, which is just below the "grey zone" of 2.3-3.0 that indicates immune activation and perhaps early signs of pathologic states (Zahorec, 2021). This likely reflects our aged study population, as age strongly predicts decreased immune function, increased immune activation, and a higher NLR(44,45). Yet, the PD patients experiencing current depressive symptoms showed a markedly higher NLR on average (3.4) than either controls (2.1) or PD patients without depressive symptoms (2.8). Furthermore, the patients with current depressive symptoms also showed more accelerated epigenetic immune system aging, which is an indicator of biologic aging and inflamm-aging. We have previously linked other immune system aging markers to PD as well (Paul et al., 2021). However, the current study demonstrates even more pronounced increases for PD patients experiencing current depressive symptoms. Interestingly, PD patients with a history of clinical depression prior to PD diagnosis had immune cell profiles similar to PD patients without a history of clinical depression. This result seems to indicate that the differences in leukocyte composition of the blood are related to current depressive symptoms and health states rather than reflecting a lasting pathology of clinical depression without current symptoms.
Although current depressive symptoms were necessary to see activated immune cell composition profiles, PD patients with a history of clinical depression showed differences in CpG levels across the genome that were unrelated to blood cell composition. We identified 35 specific CpGs associated at an FDR<0.05 and 569 at a relaxed significance threshold for discovery (FDR<0.10). Several of the CpGs are noteworthy, including cg23426156, cg11042505, and cg21769117. Two, cg23426156 and cg11042505, have been linked to the important neuropsychiatric disorder schizophrenia in three different study populations (Hannon et al., 2016(Hannon et al., , 2021. While cg23426156 is intergenic, it is in the transcription factor binding site for several factors (i.e., CEBPB, STAT3, JunD, GR, RFX5_(N-494), Rad21, p300_(N-15), SMC3_(ab9263), Pol2-4H8, c-Fos). The CpG cg21769117 on the other hand is in the transcription start site region of the gene CLIC1. CLIC1 and the chloride intracellular channel 1 (CLIC1) protein have been compellingly linked to neurodegeneration. For instance, during chronic inflammatory states in the central nervous system, CLIC1 increasingly accumulates in peripheral blood mononuclear cells, as shown in Alzheimer's patients (Carlini et al., 2020). Proteomics analysis of plasma from PD patients and controls also found higher levels of CLIC1 in PD patients (Dong et al., 2019). CLIC1 has also been implicated in microglia-mediated β-amyloid peptide neurotoxicity (Skaper et al., 2013) and IL-1β biology (Domingo-Fernández et al., 2017), connecting the gene and its protein with inflammatory neurodegenerative processes. Our pathway and enrichment analyses further implicated several immune function pathways with a history of clinical depression among the patients, including among KEGG pathways, cytokine-cytokine receptor interaction and leukocyte transendothelial migration, according to Reactome pathways, immunoregulatory interactions between a lymphoid and a non-lymphoid cell, interleukin-4 and interleukin-13 signaling, and neutrophil degranulation, and from Panther pathways, FAS (death receptor signaling on cytotoxic T cells and NK cells) and Integrin (principal cell adhesion receptors that mediate leukocyte migration and activation) signaling pathways.
To assess how these CpGs identified as associated with a history of depression in PD may be related to other biologic layers, we further linked CpGs to quantitative trait loci and determined whether meQTLs were eQTLs based on GTEx. We found that 45 single nucleotide polymorphisms (SNPs) were associated with 25 of the 35 CpGs associated at FDR<0.05. Furthermore, many of these methylation QTLs were also eQTLs such that the SNP variant was related to both methylation levels of one of the EWAS CpGs and expression-based abundance of different transcripts measured in both blood and brain regions. This is expected, as many if not most cis-eQTLs occur at the same genomic location as a cis-meQTL(50). One of the best described functions of methylation is regulating gene expression (Moore et al., 2013) and our QTL analysis seems to support this, closely linking genetic variation, methylation, and expression. For instance, one of the top hits, cg26297819, is in the transcription start site region of TEF. One meQTL was associated with the CpG, rs202637. Not only has this SNP previously been associated with PD in a meta-GWAS (p < 0.05 (Lill et al., 2012;Nalls et al., 2014)), but TEF variants have been associated with depression in PD as well (Kripke et al., 2009;Hua et al., 2012). The meQTL was also strongly related to transcript abundance of not only TEF measured in blood, but also CSDC2, POLR3H, and MEI1 in brain, transcripts which have been linked to depression in a large meta-GWAS(51).
The differentially methylated region (DMR) and QTL analysis also converged on two regions of relevance for our outcome. First, having a clinical history of depression among the PD patients was related to a DMR in chromosome 1, in a genetic region near PM20D1, SLC41A1, RAB29, and NUCKS1. The most significant individual CpG related to depression in the region we identified was cg15199181 (FDR = 0.02). This CpG was associated with three meQTLs, most significantly rs823114. This SNP, in the NUCKS1 gene, was strongly associated with PD in the meta-GWAS (meta p-value = 1.78E-13 (Lill et al., 2012;Nalls et al., 2014)). Furthermore, this meQTL is also an eQTL associated with expression levels of PM20D1, RAB29, SLC41A1, and NUCKS1 in blood, RAB29 in the basal ganglia, frontal cortex, and substantia nigra, and PM20D1 in the frontal cortex, i.e., brain areas of interest related to Table 5 History of clinical depression EWAS CpGs and Functional Genetics. QTL analysis, methylation QTLs (cis-meQTL) SNPs → CpGs (Dutch Biobank) and SNPs →expression (GTEx). Expression based on transcripts in samples from three brain regions (putamen basal ganglia, frontal cortex, and substantia nigra) from GTEx.  Purlyte et al., 2018). This association supports the notion that PD genetic risk loci may also influence symptom profiles among patients, which we have previously reported (Paul et al., 2018b) and also observe here with the PD GWAS SNP rs823114 influencing the methylation levels of cg15199181. Importantly, our findings suggest that part of the genetic risk may be conferred through an influence on mechanisms involving methylation and expression levels.
A history of clinical depression in PD was also related to a second DMR in chromosome 6 in the MHC class III region, near cg21769117, the CpG in the CLIC1 gene. The CpG was associated with four meQTLs that are also eQTLs associated with expression levels of 18 different transcripts in blood and 9 different transcripts in the brain (basal ganglia, frontal cortex, or substantia nigra). These transcripts include CLIC1, many other immune related proteins (HLAs, complement proteins), and notably CYP21A1P expression in the substantia nigra, which we determined based on our eQTL analysis with GTEx. cg00412337 in the nearby CLIC5 gene was also positively associated with depression. While these results are interesting, the complexity of the MHC class III region with hundreds of genes, requires additional evidence before drawing conclusions about its effects on depression in PD.
When we considered the complete gene set from all implicated CpGs, QTLs, and blood-based transcripts, all enriched biologic processes were related to immune function (e.g., interferon-gamma-mediated signaling pathway, antigen processing and presentation pathways, and immune response-activating pathways). Many of the 21 transcripts expressed in the brain regions were also related to lymphocyte mediated immunity and adaptive immune response, along with biosynthesis of mineralocorticoids and glucocorticoids via CYP21A2 and its pseudogene CYP21A1P. Glucocorticoids are immunoregulatory hormones generally synthesized in the adrenal cortex, where most CYP21A2/CYP21A1P expression occurs (Ahmed et al., 2019). However, glucocorticoid synthesis has been observed in the brain (Ahmed et al., 2019), seemingly allowing local regulation of immunologic response, and the GTEx analysis indicates both CYP21A2 and CYP21A1P are expressed in the basal ganglia, frontal cortex, and substantia nigra, with some variation in expression determined by eQTLs of interest. Furthermore, CYP21A1P is a transcribed pseudogene, or ancestral copy of the protein-coding gene that are thought to have lost the functional product of their parental gene due to accumulation of mutations (Milligan and Lipovich, 2015). While not all pseudogenes are transcribed, for those that are, they can provide a key mechanism for regulating the parental gene's expression (Milligan and Lipovich, 2015;Li et al., 2013). Overall, this differentially methylated region associated with clinical depression in PD seems to have genetic contributors (meQTLs/eQTLs) related to expression of important immune genes both peripherally and centrally.
Ultimately, our analysis demonstrates the dynamic and interdependent nature of biologic systems. Our EWAS linked distinct CpGs to depression among PD patients, functional genomics and QTL mapping broadened the scope of our investigation and strongly implicated the immune system in neurodegeneration and depression in PD. In the future, multi-omics measurements from the same study participants may allow us to link these QTLs and methylation levels to measured transcripts specifically in PD patients. This analysis also suggests that blood is an important tissue that provides clues into mechanisms that contribute to depression and neurodegenerative disease. Our analysis strongly implicated immune system function as being disrupted among PD patients with a history of clinical depression as suggested by differential leukocyte methylation.
Overall, the current study maps methylation signals associated with depression among Parkinson's disease. The findings provide evidence of immune system involvement in depression among Parkinson's patients, both for those experiencing current depressive symptoms and those with a history of clinical depression. The first expressed in leukocyte composition and the second in methylation level differences in specific immune system related chromosomal areas. These may be both a consequence of disease pathogenesis and a contributor to its progression. Future research should investigate whether such signals are specific to PD or also related to depression in other neurodegenerative disorders, such as Alzheimer's. PD is a highly heterogenous disorder with patients experiencing different symptom profiles, often including depression. Comprehensively mapping biologic pathways and perturbations associated with distinct symptoms may shed light onto PD's complex clinical heterogeneity and etiology.

Study population
This study was based on participants of the Parkinson's disease, Environment, and Genes (PEG) study (Narayan et al., 2013), a population-based PD study from three counties in California's Central Valley (Kern, Fresno, and Tulare). PEG was designed as a case-control study to investigate PD etiology (2001-2007 & 2010-2016; n = 849 PD patients early in disease; n = 1021 population-based controls). Further information on the patient population has been published (Duarte Folle et al., 2019). Informed consent was obtained from all subjects and the study was approved by the UCLA institutional review board. For this analysis, data were restricted to 465 PD patients and 227 controls of European ancestry with methylation data. All patients in the study were examined by movement disorder specialists (lead by J.B.) at least once and confirmed as having probable idiopathic PD based on published criteria (Hughes et al., 1992a(Hughes et al., , 1992b. Demographic characteristics for PD patients with and without methylation data were similar, including average age 70.4 (SD = 11.7) vs 70.5 (SD = 9.8) and 62% vs 65% male. We assessed two indicators of depression, first, having a history of clinical depression. This indicator was based on self-reported physician diagnosis during interviews conducted by trained study staff. Second, we assessed current depressive symptoms at the time of blood draw, measured by a self-administered 15-item geriatric depression scale (GDS).

Methylation and QTL analysis
For methylation, DNA was extracted from peripheral whole blood at the baseline visit for all study participants. We profiled and processed DNA samples using the Illumina Infinium 450k platform (486k CpGs) with standard settings. We used k nearest neighbors for imputation from the Impute R package, the background normalization method from the Genome Studio software to process DNA methylation β values, and corrected for type I/type II probe bias with BMIQ using the champ. norm function in the ChAMP R package (Morris et al., 2014). For analysis we removed x-reactive CpGs (29,233 CpGs), CpGs that had a SNP in the probe (probe_rs), CpG interrogation site (CpG_rs), or the single nucleotide extension (SBE_rs) (104,206 CpGs), and CpGs that were in the X or Y chromosome (11,648 CpGs). Analysis was therefore based on 352,325 CpGs. More detail has been published (Chuang et al., 2017;Paul et al., 2021) and the data are available on Gene Expression Omnibus (GEO), accession numbers GSE72774 and GSE72776.
We estimated whole blood cell composition from the DNAm using the Houseman estimation method to estimate the proportion of CD8 + T cells, CD4 + T cells, natural killer, B cells, monocytes, neutrophils, and eosinophils (Houseman et al., 2012(Houseman et al., , 2014. The neutrophil to lymphocyte ratio was calculated by taking the ratio of neutrophils to lymphocytes (CD8 + T, CD4 + T, natural killer, and B cells). We also estimated DNAm epigenetic immune age acceleration using extrinsic epigenetic age acceleration (EEAA). EEAA is a measure of biologic aging in immune cells that is based on the Hannum clock (Hannum et al., 2013), but dependent on white blood cell concentrations (Horvath and Levine, 2015).
To establish meQTLs associated with the EWAS CpGs, we used two sources, first, the BIOS public database of meQTLs, which is based on whole blood methylation from 3,841 samples from five Dutch biobanks (Bonder et al., 2017). Second, we used meQTLs established in our PEG PD study population. SNP data in PEG were derived from genome-wide association data generated using the Global Screening Array (Illumina, Inc.) and subjected to standard quality control and genotype imputation by Minimac3 (64). For details on the genotyping and post-genotyping procedures applied to the PEG dataset see ref (Hong et al., 2020). To establish the meQTLs, we limited to the 17 EWAS CpGs (Table 2) and ran pair-wise linear regression models (i.e., every CpG ~ SNP pair) among those of European ancestry, controlling for age, sex, and AIMs fractional ancestry. CpG ~ SNP pairs with an FDR<0.05 were considered meQTLs. Two additional meQTLs were determined from this analysis. The majority of meQTLs were detected in both study populations (Supplemental Table 3).
To establish if the meQTL was also an eQTL, we used GTEx v8, which provides characterization of genetic associations and gene expression and splicing in 838 individuals, 52 tissues, two cell lines (whole-genome sequence and RNA-sequence from approximately 960 deceased adult donors; 85% European Ancestry, 66% male, mean age 53.4 (21-70), mean BMI 27.3 (both men and women)) (Aguet et al., 2017;Lonsdale et al., 2013). As expression is tissue-dependent, we determined eQTLs based on whole-blood and three brain regions of importance in PD: substantia nigra, putamen basal ganglia, and frontal cortex.

Statistical analysis
Mean differences in cell composition markers between groups (PD patients with depression versus PD patients without depression and controls without depression) were compared with a Wilcoxon test. Using logistic regression models controlling for covariates listed below, we confirmed the mean comparisons, i.e., associations remained after confounder adjustment. For the epigenome-wide association analysis (differentially methylation positions, DMP), we used the meffil R package to test for association between the binary depression indicators and each CpG site using linear regression models fit with limma (Ritchie et al., 2015). To control for potential confounding, we included the following covariates: age at blood draw, sex, smoking, number of years with PD at blood draw, AIMs-based Caucasian fractional ancestry, an indicator for PEG study wave, and blood cell composition (CD4T and neutrophil proportions). We also included an indicator for pesticide exposure in sensitivity analysis. With the EWAS DMP associations, we next performed gene set enrichment analysis using the methylGSA R package, which accounts for the number of CpGs per gene included the Illumina 450k array. We assessed enrichment of KEGG pathways, Reactome pathways, and Panther pathways. Multiple testing for both the DMP EWAS and GSEA was adjusted for with a false discovery rate (FDR). To assess differentially methylated regions (DMR), we used the ChAMP R package, applying bumphunter and 1000 bootstraps to estimate regions for which the methylation genomic profile deviates between groups. The DMR function detects differentially methylation regions between two populations (i.e., patients with and without depression), and returns the DMRs and estimated p-values. Multiple testing was adjusted for with a family-wise error rate (FWER). Finally, for gene set enrichment based on all genes implicated by the EWAS and QTL analysis, we used WebGestalt (WEB-based Gene Set Analysis Toolkit), implemented via R. We assessed over-representation of the gene set for gene ontology terms and network topology-based analysis.

Ethics approval and consent to participate
The PEG study was approved by the UCLA Institutional Review Board (IRB# 11-001530) and informed consent was obtained from all individuals. Our research conformed to the Declaration of Helsinki.