Abstract
Background:
Stochastic epigenetic mutations (SEM) reflect a deviation from normal site-specific methylation patterns. Epigenetic mutation load (EML) captures the accumulation of SEMs across an individual’s genome and may reflect dysfunction of the epigenetic maintenance system in response to epigenetic challenges.
Objective:
We investigate whether EML is associated with PD risk and time to events (i.e., death and motor symptom decline).
Methods:
We employed logistic regression and Cox proportional hazards regression to assess the association between EML and several outcomes. Our analyses are based on 568 PD patients and 238 controls from the Parkinson’s disease, Environment and Genes (PEG) study, for whom blood-based methylation data was available.
Results:
We found an association for PD onset and EML in all genes (OR = 1.90; 95%CI 1.52-2.37) and PD-related genes (OR = 1.87; 95%CI 1.50-2.32). EML was also associated with time to a minimum score of 35 points on the motor UPDRS exam (OR = 1.28; 95%CI 1.06-1.56) and time to death (OR = 1.29, 95%CI 1.11-1.49). An analysis of PD related genes only revealed five intragenic hotspots of high SEM density associated with PD risk.
Conclusion:
Our findings suggest an enrichment of methylation dysregulation in PD patients in general and specifically in five PD related genes. EML may also be associated with time to death and motor symptom progression in PD patients.
Keywords
INTRODUCTION
Parkinson’s disease (PD), the second most common neurodegenerative disorder after Alzheimer’s disease, affects approximately 2%of individuals over 60 years of age [1]. PD is clinically characterized by motor symptoms such as rigidity, resting tremor, bradykinesia, and postural instability, but little is known about factors that influence symptom worsening [2]. Neuropathology in PD is marked by the loss of dopamine producing neurons in the substantia nigra pars compacta. While evidence for genetic contributions to PD risk has been accumulating over decades [3–5], research on epigenetic processes and their role in PD etiology is still incomplete [6–8]. There is also mounting evidence that the environment and its interplay with genetic variation and susceptibility plays a role in increasing PD risk and some of this may be through epigenetic mechanisms that affect protein expression levels. For example, epigenetic processes involving methylation of cytosine nucleotides can suppress transcription of genes to proteins. Epigenetic mutations have previously been shown to be useful biomarkers of biological aging which in turn potentially contributes to both risk and progression of neurodegenerative disease [9]. Based on the hypothesis that methylation patterns across the genome play a role in biologic mechanisms that contribute to aging and disease risk, investigators have conducted epigenome-wide association studies (EWAS) to map regions to traits of interest [10]. The EWAS approach is predicated on the notion that expression levels for a limited number of genes are primary drivers for disease risk, possibly mediated by environmental factors. Our group as well as others who have investigated links between PD and methylation patterns have employed EWAS to map loci that confer risk for PD [7, 11]. For example, Chuang et al. identified 82 CpGs and 5 CpGs in blood and saliva respectively that were either hyper or hypomethylated in PD cases when compared to controls [7].
Complementary work has suggested that overall dysregulation of methylation levels, rather than changes in methylation levels at a specific set of loci, are important drivers for outcomes such as aging [12], pre-term birth [13], or cancer [14]. Recently, a metric known as a SEM (stochastic epigenetic mutation) has been suggested as one approach for quantifying dysregulation at a particular genomic site [15, 16]. SEM reflects the accumulation of epigenetic maintenance machinery failures (stochastic mutations) and the number of SEMs increases exponentially with chronological age [12]. As an example of the applicability of this metric in earlier work, we aggregated SEM status across all methylation sites into a genome-wide score we called EML (epigenetic mutation load) to investigate the relationship between methylation levels and biological clocks [17]. Here, we employ this measure (denoted in this paper as EMLgenome) to improve our understanding of PD etiology and progression. Additionally, we devise a new measure called EMLPD that is purely composed of methylation sites that map to a set of curated genes previously linked to PD risk. Our primary objective is to assess whether (log-transformed) EMLgenome and EMLPD scores are associated with PD risk, motor symptom progression, and/or time-to-death.
MATERIALS AND METHODS
Study population
We draw on blood samples from participants of the Parkinson’s disease, Environment and Genes (PEG) study, including 511 individuals of European ancestry (292 PD patients and 219 controls) and 61 individuals of Hispanic ancestry (42 PD patients and 19 controls). PEG was designed to investigate links between environmental factors and PD in agricultural regions of central California [18, 19]. Eligible patients recruited via community outreach, local neurologists, and clinics were enrolled into the PEG study in two waves, denoted in this paper as PEG1 and PEG2. Analyses related to PD progression were restricted to either patients (PEG1 and PEG2 combined) or PEG1 controls (for comparison purposes). Standardized interviews were used for collecting information such as demographic and lifestyle factors, comorbidities, and medication use. PD patients enrolled in PEG1 were recently diagnosed (within 3 years) and in PEG2 on average 4 years before they were examined by UCLA movement disorder specialists (JB, AK) at least once (i.e., at baseline). Patients and controls were ascertained from Kern, Tulare, and Fresno counties (2001–2007 for PEG1 and 2011-2015 for PEG2). PEG1 patients were identified by neurologists, large medical groups, or public service announcements. PEG2 patients were identified through the California PD Registry operating in the same counties.
Of 1,167 PD patients initially invited in PEG1, 563 patients were eligible based on the threshold of diagnosis of 3 years or less. 363 PEG1 patients were confirmed as having probable idiopathic PD through examination by movement disorder specialists. At first re-contact for follow-up examination, 83 (22.9%) of the 363 patients were deceased, 25 (6.9%) withdrew, 9 (2.5%) could not be found, and 4 (1.1%) were too ill to participate. For follow-up, 242 (66.6%) were re-examined with 6 patients participating by mail/phone only (no in-person follow-up exam and therefore not contributing to this analysis) and 3 patients were re-classified by us during follow-up as not having idiopathic PD. Of 233 patients re-examined and confirmed to have idiopathic PD, 178 (76.4%) were seen twice and 55 were seen only once during follow-up (21 were deceased prior to and 14 pending a second exams, 8 were too ill to participate, 4 could not be located, 3 withdrew, and 5 participated by mail/phone only).
From 2011 to 2015, 734 PD patients were eligible at first screening for entry into the PEG2 cohort. Of these, 85 potential PEG2 patients refused to participate, 47 could not be examined by UCLA movement disorder specialists at a local clinic, and 119 were excluded after an exam as they did not have idiopathic PD, leaving a baseline of 483 PD patients who met our enrollment criteria Of the 227 patients with methylation data available, 134 (59%) provided some follow-up data (124 of whom have motor symptom data), while 59 (26%) were deceased, 10 (4%) refused to participate, 17 (8%) could not be found, 5 (2%) are pending to have a second exam, and 2 (1%) were too ill to participate.
Population controls from PEG1 from the same communities were randomly sampled from Medicare lists (before HIPAA enactment) and through residential tax assessor’s records as well as through community outreach. The controls for PEG2 only contributed saliva samples for DNA extraction and thus are not part of this investigation. PD risk analysis was restricted to n = 806 participants (n = 238 controls, n = 341 PD patients from PEG1, n = 227 PD patients from PEG2) for whom we currently have quality control validated blood-based DNA methylation data available.
Methylation data
DNA methylation data were obtained from peripheral blood samples. Methylation status was measured in bisulfite-treated genomic DNA using the Illumina Infinium HumanMethylation-450 assay. For the main analyses assessing the relationship between EML and all PD outcomes, we included all participants from PEG1 and PEG2, controlling for study wave status. To assess robustness of our inferences, we also conducted a sensitivity analysis, limiting to only PEG1 cases for the following reasons. We computed the distribution of our predictor of interest (i.e., the SEM score as defined in the following sections) based on our reference group, i.e., the community controls sampled along cases in PEG1 only, because only PEG1 required control subjects to provide a blood sample while controls for PEG2 only provided a saliva sample. Furthermore, we observed differences in the overall methylation levels at certain sites between PEG1 and PEG2 patients, possibly indicating either technical artifacts (e.g., a batch effect) or biological differences (e.g., differences in age, environmental factors between the two birth cohorts), which could potentially skew distributions of methylation events between PEG1 and PEG2 cases. Therefore, we present results for the whole study population and also for PEG1 only.
Quality control
Raw signal intensities were retrieved from Illumina iDat format files using the function read.metharray.exp of the R package minfi from the open-source library Bioconductor [20], followed by linear dye bias correction, noob background correction, and BMIQ normalization using the same R package [21–23]. The β-values, determined by calculating the ratio of intensities between methylated and unmethylated sites, were used for all analyses. Probes on the array were excluded based on the following criteria: 845 probes that were statistically likely to be detected above the background signal level regardless of methylation status (based on alpha-level p < 0.05), 645 probes with bead counts less that 3 in at least 5%of the samples, 11,334 probes mapping to sex chromosomes, 7,306 probes that included a known SNP within the CpG interrogation site, and 27,332 probes that were cross-reactive. In total, 438,050 remaining probes were included for downstream analyses.
Outcomes
In this paper, we examined the relationship be-tween methylation pattern and four outcomes: first, overall risk of PD, then disease progression and time to death in PD patients. PD diagnosis was determined by UCLA movement disorder specialists as described above and we compared cases to our community controls. Patients were examined clinically while not receiving PD medications, i.e., being functionally in an ‘off’ state (82%at the baseline examination and 80%at follow-up). For patients whom we could only examine on medication, we estimated an off-score by adding the difference of the whole study population’s mean off- and mean on- scores at the time of exam to the patient’s on-score, as described in previous work [24]. Second, we modelled time to death as a continuous variable indicating the number of years from the baseline exam to death. The last two outcomes were based on two complementary measures, capturing PD symptom progression that have been widely used in PD research: the Unified Parkinson’s Disease Rating Scale part 3 (motor –mUPDRS) and Hoehn & Yahr staging [25], rated by UCLA movement disorder specialists who conducted physical examinations at each examination to assess PD progression. We used time to a score of 35 or higher on the mUPDRS to indicate progression to more severe PD motor symptoms, as we did in previous analyses [26]. A mUPDRS score≥35 (higher score represents worse motor function) was chosen as a meaningful threshold for motor progression because it represents, on average, motor progression to a stage where patients start presenting some dependency for functional activities, equivalent to a HY stage 3 and to 60%on the Schwab and England scale [27]. In our analysis, we considered time to conversion to mUPDRS > = 35 and time to conversion to Hoehn & Yahr stage 3 or greater (HY3).
PD-specific gene set
The Parkinson’s UK Annotation Initiative has provided the research community a free resource in the form of a list of 811 curated genes involved in over 4,500 GO annotations that relate to processes relevant to PD [28]. Based on this data, we hypothesized that methylation patterns at these genes may be informative and would be associated with PD disease risk. We intersected this list of 811 curated genes to probes annotated with the same gene symbols based on a methylation probe descriptor file provided by Illumina, resulting in 747 candidate regions.
Metrics for assessing methylation dysregulation: SEM and EML
After exclusion criteria were applied on methylation probes that failed quality control, at each subject-probe datapoint across 438,050 probes and 806 study subjects, we computed its SEM status as a binary score. Specifically, at each locus, we first determined the interquartile range (IQR) among all 238 controls. Relying on a previously published definition, a SEM is observed for a given person at a specific CpG site if an individual’s methylation level is more than three times the interquartile range (IQR) lower than the 25th percentile (Q1 – 3×IQR), or more than three times the IQR higher than the 75th percentile (Q3 + 3×IQR) [12].
We assessed whether downstream statistical inferences would remain stable when using different threshold criteria via sensitivity analyses. Namely, we also considered two generalized scores, SEM
α and SEM
β, defined as:
where we re-fit all our logistic regression models for predicting PD risk after varying values of α and β.
Based on the SEM scores defined at each locus, for each subject, we constructed three EML scores to capture the overall burden of SEM counts for selected sets of valid methylation probes for each subject [17]. We denote these EML measures as the EMLgenome, EMLPD, and EML <window > scores. The EMLgenome score is equivalent to the definition introduced in Yan et al. (2020), i.e., the EML taken over all available probes spanning the genome. A subset of probes that map to the 747 genes from the PD geneset described earlier are used to construct the EMLPD score. We hypothesized that the EMLPD score may also predict PD disease risk and compared results of analyses using both EML scores. Finally, we applied a sliding window approach in order to consider the possibility that there are certain hotspots, defined as regions mapped to a subgroup of probes within the PD geneset that are most strongly associated with PD risk. We define notation for this measure as EML <window >, where the subscript <window> (e.g., EMLChr5 :100502 - 120000) represents the physical base pair coordinates of the start and end boundaries of the sliding window.
Figure 1 suggests a strong right skewing in the distribution of the EMLgenome score. This was confirmed based on a Kolmogorov-Smirnov (KS) test, comparing to a normal distribution, for both EMLgenome and EMLPD (deviance of 1.0 (p < 1e-16)). Hence, for regression models, a log transformation was applied to all EML definitions to ameliorate violations to normality in regression model residuals. To aid in interpretation of the unit size of regression coefficients, we further applied a second transformation, mapping all log-transformed EML scores to a Z-score.
Distribution of EMLgenome comparing cases (PD status = 1) and controls (PD status = 0) for PEG study wave 1 and 2.
Statistical analyses
To study the relationship between our EML measures and PD status or disease severity (i.e., time-to-event as described above), we constructed logistic regression models to estimate odds ratios in cases and controls, and Cox proportional hazards models for estimating hazard ratios among PD cases only. We also conducted additional analyses, stratified by gender. All regression models were adjusted for potential confounding variables: age (continuous variable with 31–99), sex (binary), ethnicity (white/non-white), and five quantities that estimated blood cell composition since methylation is tissue cell type composition dependent. For the time to event (Cox) models, to minimize confounding on methylation due to disease duration, we also included a variable measuring the number of years between diagnosis and the baseline interview date. The cell composition variables distinguish between naïve CD8 cells, CD8 + CD28-CD45RA- T cells, plasma blasts, CD4 T cells, and granulocytes. These variables were estimated from an online tool (https://horvath.genetics.ucla.edu/html/dnamage/), based on methods descri-bed by Horvath [29].
To assess robustness of the results, as a sensitivity analysis we added three additional variables based on knowledge gleaned from previous work investigating potential factors for methylation changes: educational attainment, smoking status, and a biomarker measuring epigenetic age acceleration. An EWAS of 4152 participants demonstrated an education-related gradient of methylome changes that reflect cigarette smoke exposure signatures [30]. For this study, we encoded years of educational attainment and smoking status as categorical variables: (educational attainment: less than 12, 12, or greater than 12 years of education and smoking status: ever smoked, never smoked, or quit smoking). Previous studies have shown that a biomarker known as extrinsic epigenetic age acceleration (EEAA) was associated with increased mortality [31], PD risk [32], and EML [17]. Thus, given its association with EML, we considered EEAA as a potential confounder. We compared regression results from the standard adjusted model (i.e., with all control variables described above) to a model with the added continuous variable measuring EEAA, and two categorical measures of educational attainment and smoking status in the sensitivity analysis.
In models of disease severity (i.e., time to death and PD symptom progression), we also included total levodopa equivalent dosage (LED) reported at the baseline exam as levodopa medication may potentially affect methylation levels, and patients with more severe symptoms of PD may be prescribed a more aggressive levodopa regiment. Previous analyses suggested that methylation levels confer a risk of developing PD, and it has been shown that methylation is affected by genes, i.e., that the average heritability of DNA methylation sites is 18.7%[33]. Thus, we also included as a covariate in our models a polygenic risk score (PRS) to control for any heritable variation. The PRS was constructed as a weighted sum of allele counts, where the weights are effect estimates from known PD risk SNPs [34].
For the purpose of testing associations between PD and EML <window > (i.e., window-specific EML scores), we applied the R package wgscan to test a series of candidate sliding windows, comprising contiguous methylation probes, for association with either PD disease status or disease severity. The wgscan program, which was originally applied for mapping a contiguous set of rare SNPs to a given outcome [35], can be applied to other measures (e.g., CNV, methylation counts) under the assumption that a cluster of variation at neighboring sites is relevant to the trait. Methods for rare variant analysis in general apply either a burden test (i.e., scores are summed) or a test of overdispersion in the variance of the effect sizes on a selected set of contiguous variants, where the latter test can potentially be more powerful than the former if variants within a window do not share the same directionality of effects [36]. In the context of our analysis, we treat SEM calls at each methylation probe as the variant of interest. The program applies a sliding window approach, systematically considering various window sizes. Candidate window sizes can be configured by the user, and significance levels are adjusted for the increased type 1 error due to the large number of windows (hypotheses) being tested. For each window tested against the trait of interest (PD), we included in the model the standard set of control variables described above. At each window in the provided region, the program reports p-values for the burden and overdispersion tests. A global region-specific significance level, adjusted for the number of windows tested is also reported to the user. Finally, we computed a genome-wide significance level αglobal to correct for the number of regions tested. For k regions (genes), this threshold is defined as αglobal = 1/(1/αgene - 1 +1/αgene - 2 +...1/αgene - k) where the region-specific significance levels (e.g., αgene - 1) is computed to account for the multiple comparisons inherent in testing multiple sliding windows.
RESULTS
Table 1 summarizes the characteristics of the 806 subjects included in this study whose quality checked methylation data were available for analyses.
Characteristics of PEG participants
1Time span in years between first and last exams for motor progression. 2Time span in years from diagnosis until baseline exam.
Distribution of SEM frequency across loci
Inspection of the distribution of SEM frequency across all methylation probes showed that the probability of a SEM event at any particular locus was rare. A histogram of SEM counts at each methylation probe is depicted in Supplementary Figure 1. The median frequency at a given site was 0.74%(6/806) with a range of 0-211.
Estimated effects of EML genome on PD risk
A visual comparison of the distributions of the EMLgenome score between cases and controls, as depicted in the violin plot of Fig. 1, suggests that EML levels are higher in cases than controls. Furthermore, the large range in scores for the full dataset appears to be driven more strongly by female cases. To test our hypothesis that overall methylation dysregulation is enriched in PD cases, we fitted regression models with PD status as the outcome of interest.
Table 2 summarizes results from logistic regression models for all PEG1/PEG2 subjects and men and women separately. For comparison, we also present results after excluding PEG2 subjects as part of a sensitivity analysis. As expected, for effects estimated with different case inclusion or exclusion criteria (i.e., all PEG waves versus PEG1 subjects only), the directionality of effects was consistent, with somewhat diminished effect sizes for the latter comparison. To quantify differences in the distribution of EML scores between the two datasets, we regressed the EML measures on PEG wave status, controlling for the same confounding variables in the logistic regression models. We did not observe any formally statistically significant heterogeneity of effects as indicated in Supplementary Table 1.
Odds ratios of PD onset for the variables EMLgenome and EMLPD
Models control for sex, age, ethnicity, and cell counts (Gran, NK, CD8pCD28nCD45RAn, CD4T).
Risk of PD was statistically significantly associated with both EML scoring criteria when comparing all cases to controls (EMLgenome: OR = 1.90; 95%CI (1.52–2.37); EMLPD: OR = 1.87; 95%CI (1.50–2.32). The EMLgenome variable exhibited the strongest estimated effect size on PD risk for female cases only (all subjects: OR = 2.28; 95%CI 1.59–3.25; PEG1 only: OR = 2.05; 95%CI 1.45–2.92). The association between EMLgenome and PD risk was weaker in males (all subjects: OR = 1.66; 95%CI 1.22–2.24; PEG1 only: OR = 1.46; 95%CI 1.08–1.97). The magnitude of the estimated effect size of PD risk was smaller with EMLPD than with EMLgenome, but the direction of the effect was consistent, with females always showing a stronger association than males. To assess the possibility of sex specific heterogeneity, we fitted models after adding an interaction term between sex and either EML variable. Based on this interaction term, we observed moderate evidence of heterogeneity by sex for EMLgenome (OR = 1.53; 95%CI 1.03–2.26; p = 0.0338) and weaker evidence for EMLPD (OR = 1.41; 95%CI 0.96–2.06; p = 0.0787), suggesting a stronger joint effect of female sex and EMLgenome than expected.
We conducted additional sensitivity analyses on the full dataset to assess the robustness of our findings. First, adding the potential confounders EEAA age acceleration, educational attainment, and smoking status did not change our results (Supplementary Table 2). Second, when we defined EML according to SEM α and SEM α odds ratio estimates remained very similar (Supplementary Table 3).
Mapping of EML <window > to PD risk
We split the genome-wide methylation dataset into 747 regions, each of which corresponds to a gene that has previously been associated with PD and received a GO annotation as listed on the Parkinson’s UK Annotation Initiative website. A total of 9400 windows ranging from size 5 kB to 50 kB were tested across this subset of the genome.
We did not find any formally statistically significant associations between the EML <window > scores and the continuous outcomes measuring disease severity. However, burden and dispersion tests indicated that five of these EML <window > scores were associated with PD risk, based on the genome-wide significance level αglobal of 3e-7, as shown in Table 3. Notably, these regions map to the genes encoding for proteins that are related to cellular stress response (e.g., STK39), endocytosis (e.g., NEDD4), ER protein trafficking (e.g., PI4KA), and neuronal /microglia function (e.g., P2Y12). Although we configured the program to construct candidate window sizes ranging from 5kB to 50 kB, none of the scores corresponding to windows greater than 5 kB in size were significantly associated with PD risk. The Manhattan scatter plot in Supplementary Figure 2 depicts the landscape of -log10 p-values of window-specific burden (shown in red) and dispersion (shown in blue) tests.
Table3
*Chromosome.GeneSymbol.
Effect of EML genome /EML PD on PD progression and mortality
Next, we assessed whether either EMLgenome or EMLPD were associated with measures of disease progression in PD cases only, applying a Cox proportional hazards ratio time to event model for each of the three outcome variables. We did not observe statistically significant associations between the EML measures and the HY3 time to event variable. For the time to a mUPDRS score of 35 or greater, we estimated increased hazard ratios and their confidence intervals were formally statistically significant for all cases (EMLgenome HR = 1.28, 95%CI 1.06–1.56 and EMLPD HR = 1.20, 95%CI 1.08–1.56) and the female subset of cases (EMLgenome HR = 1.43, 95%CI 1.10–1.86 and EMLPD HR = 1.41, 95%CI 1.10–1.86), as indicated in Table 4. Neither EMLgenome nor EMLPD were associated with time to a HY stage > = 3 (EMLgenome HR = 1.13, 95%CI 0.94–1.37 and EMLPD HR = 1.13, 95%CI 0.94–1.36).
Hazard ratios for EMLgenome/EMLPD measures and time to event - (having a score of 35+ on mUPDRS) and time to death in PD patients
Models control for sex, age, ethnicity, cell counts (Gran, NK, CD8pCD28nCD45RAn, CD4T), total levodopa equivalent (total_led), PRS (polygenic risk score), and PD duration in years.
Likewise, hazard ratios were increased for time to death for all cases (EMLgenome HR = 1.29, 95%CI 1.11–1.49 and EMLPD HR = 1.29, 95%CI 1.12–1.48) and the female subset (EMLgenome HR = 1.33, 95%CI 1.07–1.64 and EMLPD HR = 1.31, 95%CI 1.07–1.60), also presented in Table 4. Effect sizes did not change substantively in models including different sets of control variables (Supplementary Table 4). To assess whether this result was specific to PD cases, we also assessed time to death between EMLgenome and EMLPD in controls. While we did not observe associations in controls that were formally statistically significant, the effect estimates were above the null value of one (EMLgenome HR = 1.21, 95%CI 0.85–1.74 and EMLPD HR = 1.22, 95%CI 0.86–1.73).
DISCUSSION
We investigated whether dysregulation of blood cell methylation plays a role in PD etiology and contributes to mortality in patients with PD in a population-based study of patients and their community controls from California’s Central Valley. For these participants, we have been able to gather a rich set of outcome data with patients followed for progression, and all participants for mortality over decades. Using a high-density assay from Illumina, we queried over 450 k probes after extensive QC filtering to assess SEM, a surrogate for measuring methylation dysregulation. Our findings suggest an increase ranging from 66–128%in risk of PD for each standard deviation of EML, with a stronger estimated effect size in women than in men. Although our analyses were based on pooled cases from two different birth cohorts (PEG1 and PEG2) compared against only one reference group (PEG1 controls), our sensitivity analyses (as indicated in the Table 2 and Supplementary Table 1) suggest that our findings are robust against confounding due to a cohort effect in that effect sizes were diminished in the smaller case control set but generally findings and trends remained the same (i.e., EMLgenome had stronger associations than EMLPD, and females showed stronger associations than males).SEM events in blood cells are stochastic in nature, affect the genome overall, and are driven by some external process but seem to nevertheless not only affect PD risk but also motor decline especially among female patients as well as time to death in all patients. Similar to PD risk, for time to death and motor symptom progression, we observed positive associations with both EML measures (i.e., EMLPD and EMLgenome).
Using an agnostic sliding window analysis, where only contiguous probes within a window were used to construct a score called EML <window > from all PD risk genes previously identified, several hotspots of dysregulation emerged as genome-wide statistically significant according to both a burden and a dispersion test for the region. Stochastic epimutations in and around five PD risk genes were identified as associated with PD risk. The serine/threonine kinase (STK39) phosphorylates alpha-synuclein, and a hyperphosphorylated state has been implicated in PD [37]. NEDD4 protein, which attaches ubiquitin chains to alpha synuclein, is critical to marking potential synuclein aggregates for degradation by the lysosomal machinery [38]. Deficiency in expression of this gene may be a driving factor in the alpha-synuclein inclusions that are a hallmark of PD. Knockout experiments in mice of the phosphoinositide kinase PI4KA resulted in diminished myelination in Schwann cells [39]; notably the authors noted that these mice displayed a neurodegenerative-like phenotype. The PI4KA enzyme can also play an important role in the immune system, where it critically binds to the CD7 glycoprotein to participate in proper intracellular signaling that modulates T cell functionality [40]. Injured neurons, which release nucleotides (purines), normally trigger an inflammatory response where microglia migrate toward the source of these nucleotides, mediated via the purinergic receptor P2Y12 [41]. Since neurodegenerative diseases such as PD are characterized by pro-inflammatory responses [42], over-expression of P2Y12 may be a driver of disease risk and/or progression. The Dickkopf-related protein 3 (DKK3), which is involved in the Wnt receptor signaling pathway, has been suggested as a therapeutic target for PD [43]. Both the ‘classic’ (through β-catenin) and ‘non-classical’ pathways play roles in dopaminergic cell development and synaptic function. PD-associated proteins encoded by PARK2 (protein: parkin) and LRRK2 have been shown to modify classic Wnt signaling [44].While our study did not find evidence for epigenetic modifications in these PD genes, we identified other genes (APC and AXIN1) involved in classic Wnt signaling in enrichment analysis. Interestingly, decreased Wnt signaling has also been reported in AD [45]. CALM2 is one of the AD-related genes we associated in blood based methylation analysis with PD status [7] and it has a function in calcium ion binding which is important in the non canonical Wnt signaling pathway [46].
We note some caveats to our results and interpretations. Reported p-values for the various regression models were not corrected for the number of outcomes tested as these outcomes are highly correlated as measures of disease progression/risk and applying a Bonferroni correction may be too stringent. Although we estimated relatively strong effect sizes for PD risk, the heterogeneity in effect sizes between sexes and their wide confidence intervals need to be addressed in further studies with a larger sample size. For example, a lack of a statistically significant association between EML and time to death in controls may have been due to the small sample size relative to cases. To consider the possibility of heterogeneity between cases and controls, we fitted a time to event model for cases and controls that included an interaction term between case status and EML and the interaction term was not statistically significant (p > 0.5). We should note that not all PEG patients with methylation data had follow-up data for mUPDRS scores (time to death follow-up data was complete however). For this reason, we considered the possibility that hazard ratios estimated for time to an increase of 35+ in the UPDRS score may have been biased, as our outcome of interest (EMLgenome) was associated with probability of follow-up based on a logistic regression. After computing case weights based on the inverse probability censoring weights method, we refit our Cox proportional hazards model that incorporated these weights. The revised hazard ratios for this model, however, remained stable at 1.28 (95%CI: 1.0451–1.564) at a precision of 3 significant figures. Independent replication will be crucial in determining whether and how blood methylation dysregulation contributes to PD risk, motor symptom progression, and/or time to death. Notably, as these signatures of methylation dysregulation were measured in blood and not brain tissue, they may suggest the involvement of peripheral organs such as the immune system in progression and time to death in PD. Also, as methylation patterns were assessed at only one time point in our study after PD onset, we cannot determine directions of action, i.e., whether disease onset caused these blood cell methylation changes or whether the methylation dysregulation contributed to PD risk representing a potential biomarker of processes leading up to PD. It is plausible that SEMs could accumulate in higher density in certain regions (e.g., the five regions identified in the sliding window analysis) in blood and other tissues of patients if certain repair mechanisms fail. Indeed, it has been shown that epigenetic fingerprints in blood can reflect upstream epigenetic changes in other tissues [47, 48]. We have previously observed a significant increase in DNA methylation age in PD patients compared to controls, as assessed by measures of both intrinsic and extrinsic epigenetic age acceleration of blood [32]. In our recent study using the same methylation dataset, we found evidence for a correlation between EML and dysfunction of the epigenetic maintenance system as captured by various measures of methylation age [17].
Given the mounting evidence that methylation plays a major role in aging and aging-related diseases, findings from this study provide interesting implications for potential mechanisms contributing to PD onset, which represent an interplay between genes and the environment. Our results suggest that certain dysfunctions in epigenome repair mechanisms affecting PD gene methylation and pathways may not only be critical in keeping the blood methylome stable but may also influence motor symptom progression in women and time to death in all PD patients.
Footnotes
ACKNOWLEDGMENTS
This work was supported by the National Institute of Environmental Health Science (grant numbers R21ES024356, 2R01ES010544, U54ES012078).
CONFLICT OF INTEREST
The authors have no conflict of interest to report.
