Alpha-synuclein overexpression induces epigenomic dysregulation of glutamate signaling and locomotor pathways

Abstract Parkinson’s disease (PD) is a neurological disorder with complex interindividual etiology that is becoming increasingly prevalent worldwide. Elevated alpha-synuclein levels can increase risk of PD and may influence epigenetic regulation of PD pathways. Here, we report genome-wide DNA methylation and hydroxymethylation alterations associated with overexpression of two PD-linked alpha-synuclein variants (wild-type and A30P) in LUHMES cells differentiated to dopaminergic neurons. Alpha-synuclein altered DNA methylation at thousands of CpGs and DNA hydroxymethylation at hundreds of CpGs in both genotypes, primarily in locomotor behavior and glutamate signaling pathway genes. In some cases, epigenetic changes were associated with transcription. SMITE network analysis incorporating H3K4me1 ChIP-seq to score DNA methylation and hydroxymethylation changes across promoters, enhancers, and gene bodies confirmed epigenetic and transcriptional deregulation of glutamate signaling modules in both genotypes. Our results identify distinct and shared impacts of alpha-synuclein variants on the epigenome, and associate alpha-synuclein with the epigenetic etiology of PD.


Introduction
Alpha-synuclein (aSyn) plays crucial roles in neurodevelopment, neuronal health, and synaptic transmission (1). Mutations, multiplications, and single nucleotide polymorphisms (SNPs) in SNCA, the gene encoding aSyn, are associated with Parkinson's disease (PD), an ageassociated, and therefore increasingly prevalent, neurological disorder (2,3). SNCA variants have also been linked to several other neurodegenerative conditions, including Lewy body dementia (LBD), multiple system atrophy (MSA), and Alzheimer's disease (4,5). The full range of functions of aSyn are still unclear, as a number of different cellular roles have been reported, including SNARE complex assembly, regulation of neuronal differentiation, glucose level, and dopamine biosynthesis, as well as modulation of calmodulin activity/G protein-coupled receptor kinase activity (1). aSyn point mutations can impair one or several of these processes, and thus disrupt neuronal health. Interestingly, we recently reported that nuclear localization and phosphorylation modulate the pathological effects of aSyn (6). However, the molecular mechanisms leading to PD may differ according to the variant. For example, the A30P mutation occurs in the membrane binding domain of aSyn, which may affect its ability to act as a presynaptic chaperone and interact with membrane-bound receptors (7). Although A30P aSyn appears less likely to aggregate than wild-type (WT) aSyn, this mutation can still result in PD, possibly due to a loss of aSyn function (7). Conversely, multiplications of the SNCA gene can lead to excess protein production, promoting aggregations and fibrils that impair synaptic function and can lead to neuronal death (8). Duplications and triplications of SNCA have been reported in patients with familial PD, LBD, and MSA (2,4).
Despite progress in the understanding the relations between aSyn dysfunction and PD, much remains to be elucidated. For example, although aSyn-and PD-associated transcriptional deregulation have been demonstrated in cell culture, animal models, and human patients, the mechanisms underlying these transcriptional changes are not yet fully understood (9)(10)(11)(12)(13)(14)(15)(16). As different aSyn mutations can lead to PD through different mechanisms, it is also important to understand which aspects of these variants are unique and which are shared. Studies by our group and others have shown that expression of WT or A30P mutant aSyn in human neurons is associated with DNA damage and transcriptional dysregulation of genes involved in cell survival and DNA repair, which may be mediated by altered histone H3 acetylation (9,10). In addition, altered DNA methylation (DNAm) patterns have been observed in both the SNCA gene itself and genome-wide in blood, saliva, and brain of PD patients (17)(18)(19)(20)(21). Taken together, these observations suggest that genome-wide transcriptional and epigenetic dysregulation may be consequences of aSyn overexpression and could potentially be involved in PD susceptibility and pathogenesis.
Epigenetic studies of PD provide the opportunity to assess genetic and environmental inf luences on disease risk concurrently (22)(23)(24)(25)(26)(27)(28). DNA methylation (DNAm), which refers to the attachment of a methyl group to DNA on the 5 carbon of cytosine, frequently in the context of cytosine-guanine (CpG) dinucleotides, is a well-studied epigenetic mark that changes during development and is inf luenced by genes and environment (22)(23)(24)(25)29,30). Although DNAm patterns can be unrelated to mRNA expression patterns or laid down as a consequence of gene expression, in some cases DNAm can also impact transcription by altering the ability of transcription factors to bind to gene regulatory regions, changing chromosomal interactions, and inf luencing splicing (22,(29)(30)(31). This can alter the regulation of pathways involved in disease susceptibility (32)(33)(34). In addition, DNAm patterns can be useful as biomarkers of aging, disease states, or exposures regardless of association with transcription (35)(36)(37)(38)(39)(40)(41).
Although DNAm is inf luenced by a myriad of factors that differ between individuals, the existing literature on DNAm in PD has primarily adopted a casecontrol approach (17)(18)(19)(20)(21). In addition, the role of DNA hydroxymethylation (DNAhm), which accounts for up to 40% of all modified cytosines in brain tissue, remains unclear (42). DNAhm is an oxidized form of DNAm, and an intermediate of the DNA demethylation process; however, DNAhm has also been suggested to be a stable, independent epigenetic mark in the brain (43). The stability of DNAhm and its ability to influence transcription by eliminating DNAm-protein interactions, introducing DNAhm-protein interactions, altering chromatin state, and altering splicing suggest that it could impact the etiology of neurological diseases (44,45). It is important to distinguish DNAhm from DNAm, as conventional bisulfite conversion techniques measure both DNAhm and DNAm in one compound signal, which may bias interpretations (46). DNAhm is stable in postmortem formalin-fixed tissue, and some initial studies characterized DNAhm patterns in the brains of deceased PD patients (47). Bulk DNAhm levels are unchanged in the cortex, substantia nigra, and brainstem of PD patients, while DNAhm levels and expression of TET2, the enzyme that forms this mark, are elevated in neurons of the prefrontal cortex of PD patients (48,49). Importantly, DNAm and DNAhm patterns can be influenced by genetics, and genetic alterations at loci, such as SNCA, influence PD susceptibility, suggesting that it will be important to incorporate genetic considerations into future epigenetic studies of PD (50)(51)(52). To date, only a few studies have assessed the impact of genetic background on DNAm patterns in PD (53- 55).
Accounting for the range of genetic, environmental, and lifestyle factors that can influence DNAm and DNAhm are major challenges in human epigenetic studies of PD (22,23,29,50). Furthermore, the use of postmortem tissue to study PD also typically limits the analysis to cases of advanced disease. Model systems, such as rodents or cell culture, can help to address these issues, providing contexts in which both genetics and environment can be tightly controlled, and allowing the examination of early impacts of neurological disease in the brain.
In this study, we characterized the influence of two molecularly distinct aSyn variants, WT and A30P aSyn, on the DNA methylome and hydroxymethylome of human dopaminergic neurons. We profiled genome-wide DNAm and DNAhm patterns in three groups of Lund human mesencephalic (LUHMES) cells differentiated into dopaminergic neurons: control LUHMES cells, LUHMES cells overexpressing WT aSyn at levels similar to those seen in SNCA multiplication carriers (WT aSyn cells), and LUHMES cells overexpressing A30P aSyn (A30P aSyn cells). In addition, we integrated our DNAm and DNAhm data with transcriptomic data from the same cells, incorporating H3K4me1 ChIP-seq to score DNA(h)m changes across gene regulatory features (9). WT and A30P aSyn expression were associated with thousands of DNA(h)m changes, particularly affecting genes that regulate locomotor behavior and glutamate signaling pathways. These observations suggested that familial PD-associated aSyn mutations have widespread epigenomic effects, and may contribute to molecular and transcriptional dysregulation associated with the etiology and heterogeneity of PD.

aSyn overexpression altered genome-wide DNAm patterns in dopaminergic neurons
We used differentiated LUHMES cells as a model system in which to evaluate the epigenome-wide impacts of  aSyn expression in dopaminergic neurons. Genome-wide DNAm and DNAhm were analyzed in the same cells as used in a previously published transcriptomic analysis by our group (GSE89115, GSE181126) (9). Previously, we verified dopaminergic neuronal differentiation in these cells by immunostaining, which demonstrated the presence of tyrosine hydroxylase (TH)-positive axonal and dendritic networks, as well as endogenous aSyn expression (Supplementary Material, Fig. S1). Proliferating LUHMES cells were infected with lentiviral constructs containing IRES-GFP, WT aSyn-IRES-GFP, or A30P aSyn-IRES-GFP, and positive clones were selected by f luorescence activated cell sorting (FACS). As expected, the levels of SNCA mRNA and aSyn protein expression were higher in WT aSyn cells and A30P aSyn cells than in control cells, as determined by RNA-seq (approximately 4-fold in both genotypes) (Fig. 1A, B, Supplementary Material, Table S1) and immunoblotting (approximately 6-fold in both genotypes, see ref. 9). Elevated SNCA mRNA expression was primarily driven by three transcripts in WT and A30P aSyn cells (Fig. 1B, Supplementary Material , Table S1).
To investigate the impact of WT and A30P aSyn expression on genome-wide DNAm patterns, we assessed DNAm differences between control, WT aSyn, and A30P aSyn cells at 813589 EPIC probes. An effect size threshold of |delta beta| ≥ 0.05 was selected to be well above technical noise (maximum 2.2% RMSE between technical replicates), and probes were considered to be statistically significant at a Benjamini-Hochberg adjusted p (padj) ≤ 0.05, corrected for multiple comparisons. First, we examined whether overexpression of aSyn was associated with reduced DNAm at the first intron of SNCA, as reported previously in PD patients (20,21). We expected DNAm to be lower at this region in aSynexpressing LUHMES cells as intron 1 has been established as a putative promoter for SNCA containing GATA binding sites, and enzymatic removal of DNAm in this region is associated with increased levels of SNCA mRNA and protein expression (56,57). In addition, to specifically examine DNAm levels at the endogenous gene without confounding by CpGs also located on the transgene, we removed five CpGs mapping to the cDNA sequence used in the aSyn constructs (NM_000345) and the remaining 35 SNCA CpGs are shown in Fig. 1C, D. As expected, both aSyn genotypes had significantly lower levels of DNAm at several CpGs in the region corresponding to intron 1 for the majority of SNCA transcripts, and upstream of the transcription start site (TSS) for several other transcripts (Fig. 1C, D). This region covered binding sites for GATA1, EZH2, TAL1, and EGR1 (Fig. 1D).
We next assessed genome-wide DNAm changes induced by overexpression of each aSyn variant. First, we used a site-specific approach, with linear modeling applied to all CpGs passing quality control (QC) filters (Fig. 2, Table 1). In comparison to controls, WT aSyn cells had 18 521 sites with decreased DNAm and 10 812 sites with increased DNAm (padj ≤ 0.05) ( Fig. 2A Table S5). To confirm an example from our findings by pyrosequencing, we selected a region of the TUBA8 gene that included the top two probes with the greatest changes in DNAm across both genotypes; differential DNAm at these two positions, one additional position on the array, and four positions not on the array were all confirmed (  Table S5).
The proportion of differentially methylated probes shared across all comparisons and any two of three comparisons was greater than expected by chance, and in most cases, the differences from controls were more pronounced for WT aSyn cells than for A30P aSyn cells (padj < 0.001, 10 000 permutations) (Fig. 2B, Table 1). Interestingly, WT and A30P aSyn overexpression primarily affected CpGs outside of genes and enhancers as defined by H3K4me1 ChIP-seq (Fig. 2C). However, differentially methylated probes were enriched in TSS1500 regions for both genotypes (statistically significant in WT only), and trends for differentially methylated probe feature enrichment differed between genotypes for 3 UTRs and CpG North Shores (Fig. 2C, D).
Following the investigation of differential DNAm patterns at individual CpG sites, we additionally examined genotype-dependent DNAm differences at comethylated regions (CMRs) (Supplementary Material, Fig. S3). CMRs are defined as regions with correlated DNAm across individuals, and such regional patterns of DNAm are more likely to be biologically meaningful than DNAm at single CpGs (58). We defined 57 941 custom DNAm CMRs in LUHMES neurons and applied  Fig. S3A, B). Approximately 40% of the total differentially methylated probes discovered with the site-specific approach mapped to CMRs in any comparison, and overlaps between differentially methylated CMRs and individual CpGs were more likely than expected by chance

aSyn overexpression altered genome-wide DNAhm patterns to a lesser extent than DNAm patterns
After assessing genome-wide DNAm alterations, we examined differential DNAhm at 233440 sites that passed our detection threshold of 3.6% (calculated according to the 95% quantile of negative hydroxymethylation (hmC) values after subtraction; see Methods). Fewer changes were observed in DNAhm for each comparison compared to DNAm, most of which were increases; these increases in DNAhm were significant when tested for skewness (p < 2.2e −16 ) (Fig. 3, Table 2,  Supplementary Material, Tables S6, S7). Only one probe was differentially hydroxymethylated in comparison between WT aSyn and A30P aSyn cells (cg11833293 in the SHANK2 3 UTR, padj = 0.045, delta beta = 0.089) (Fig. 3A, B). The overlaps of 16 CpG sites between control vs. WT aSyn and control vs. A30P aSyn were not greater than expected by chance (Fig. 3B, padj > 0.05, 10 000 permutations). Both genotypes also had fewer differen- tially hydroxymethylated probes in the first exon than expected by chance (Fig. 3C). When defining CMRs in the DNAhm data and applying linear modeling, we also found greater increases in CMR DNAhm in control vs. WT aSyn and control vs. A30P aSyn comparisons (Supplementary Material, Fig. S4). Although fewer CMRs were found in DNAhm data overall, differentially hydroxymethylated CMRs were also more likely to overlap with differentially hydroxymethylated CpGs from site-specific analysis than expected by chance (Supplementary Material, Fig. S4A).

WT and A30P aSyn altered DNAm at locomotor and glutamate signaling pathway genes
To identify possible functional consequences associated with the DNAm and DNAhm changes in each genotype, we performed Gene Ontology (GO) enrichment analysis on all differentially methylated or hydroxymethylated genes within each comparison using over-representation analysis in ermineR (59). Thirty-four GO biological processes were enriched in differentially methylated sites shared between the control vs. WT aSyn and control vs. A30P aSyn comparisons (Fig. 4A, padj < 0.05). A total of 442 CpG sites were annotated to the top GO term, 'locomotory behavior,' and differentially methylated in at least one genotype (420 sites in WT aSyn cells and 89 sites in A30P aSyn cells) (Fig. 4B-E). A total of 129 differentially methylated sites were annotated to the second highest ranking GO term, 'glutamate receptor signaling pathway' (127 sites in WT aSyn cells and 18 sites in A30P aSyn cells) (Supplementary Material, Fig. S5). No significant multifunctionality-corrected enrichments were observed for differentially methylated or hydroxymethylated sites in any of the other groups.
We additionally employed SMITE, which captures modules of functionally related genes with changes to at least one of DNAm, DNAhm, or gene expression in each comparison (61). We included H3K4me1 ChIP-seq data from the same LUHMES cells used for DNA(h)m and mRNA analyses in our SMITE workf low, allowing us to consider DNAm and DNAhm changes at enhancers in addition to the default promoter and gene body regions, and weighted the importance of each modification (expression: 0.4; DNAm: 0.35; DNAhm: 0.25). Eighteen modules were identified in WT aSyn cells, encompassing functions including DNA damage repair, cell cycle control, neuronal differentiation, chaperone activity, and cell signaling (insulin, glutamate, and PPAR) (Fig. 5A,  Supplementary Material, Fig. S7). Twenty-four modules were identified in A30P aSyn cells, which had similar functions to the WT modules as well as some additional modules related to the urea cycle, sumoylation, PDGF signaling, and K + channel transport (Fig. 5B, Supplementary Material, Fig. S8). Modules from both comparisons were primarily driven by gene expression, with some contributions from DNAm and/or DNAhm (Fig. 5). This effect was consistent even when DNAm or DNAhm was weighted higher in the SMITE model (Supplementary Material, Tables S11, S12).
Glutamate receptor signaling modules discovered using SMITE had more concurrent changes to DNAm, DNAhm, and transcription in both genotypes. Two modules centered around GRIK2 were each significantly altered in WT and A30P cells, consistent with the direct comparison of differentially (hydroxy)methylated versus differentially expressed genes. GRIK2 promoter DNAm and mRNA levels were negatively correlated in both genotypes (Fig. 6). Several differentially methylated GRIK2 promoter CpGs in both genotypes also belonged to a CMR (chr6: 101846707-101 846 916), which contained five individually significant CpGs from the control vs. WT aSyn comparison and one individually significant CpG from the control vs. A30P aSyn comparison.

Discussion
It is important to elucidate the mechanisms underlying the development and progression of PD in different individuals to develop preventative strategies and treatments. Here, we explored whether DNAm and DNAhm contributed to interindividual differences in PD etiology among carriers of aSyn variants. We assessed the impacts of overexpressing WT or A30P aSyn on genomewide DNAm and DNAhm patterns in dopaminergic neurons, the primary cell type affected in PD. Our results showed that overexpressing either aSyn variant induced thousands of DNAm changes and hundreds of DNAhm changes in pathways related to PD and neurodegeneration, and that WT aSyn particularly impacted glutamate receptor signaling genes at both epigenetic and transcriptional levels. The distinct characteristics of each aSyn protein may explain why both similar and unique effects on DNA(h)m were observed. This study enhances our understanding of the wide-ranging genomic impacts of different aSyn forms and provides further insights into the possible molecular mechanisms of PD.
The association of aSyn expression with differential DNAm of genes associated with locomotor behavior observed here was consistent with previous reports that aSyn plays a role in dopamine biosynthesis and that dopaminergic neuron activity has impacts on locomotor phenotypes (62,63). The altered DNAm, DNAhm, and expression of neuronal differentiation genes in both genotypes also agreed with previous work in SNCA triplication PD patient-derived induced pluripotent stem cells (iPSCs), which indicated reduced expression of genes involved in neuronal differentiation and impaired differentiation capacity associated with aSyn overexpression (64). Several genes identified in our SMITE network analysis have altered epigenetic regulation in the brains of PD patients, including DLG4 (differentially methylated in PD cingulate gyrus) and MAGI2 (differentially methylated in the frontal cortex and blood of PD patients, and shows increased H3K27ac in the prefrontal cortex of PD patients) (13,15,17). Glutamate signaling pathway genes have also been implicated in PD, primarily through excessive glutamatergic transmission associated with excitotoxicity in dopaminergic neurons (65). Consistent with these observations, one of the differentially methylated genes shared between WT aSyn-expressing LUHMES neurons and mouse hippocampus was a metabotropic glutamate receptor, and differentially methylated and hydroxymethylated genes in WT aSyn-expressing mice had similar functions to genes with differential methylation, hydroxymethylation, and/or expression in our SMITE analysis. These observations suggested that although DNAm and DNAhm patterns are highly cell type-specific and can change with age and environment, some of the major signaling, transcriptional, and metabolic pathways altered by aSyn expression in this model are more likely to have functional relevance across different physiological contexts, possibly including the human brain. In addition, many of the CpGs identified in site-specific differential (hydroxy)methylation analyses belonged to CMRs, including SNCA and GRIK2, suggesting that these CpGs are likely to impact regulation of glutamate receptor genes and other pathways affected by aSyn as a unit (58).
Ionotropic and metabotropic glutamate receptor expression are at least in part epigenetically controlled, as transcription of both receptor types is modulated by H3K4 methylation across the human lifespan (66). Signaling cascades activated as a result of glutamate receptor stimulation can also result in gene expression changes, possibly including expression of glutamate receptor genes themselves (67). Therefore, it is possible that aSyn-induced changes in DNAm of the promoters of glutamate receptor genes may be one mechanism underlying the regulation of their expression and subsequent glutamate release, thus facilitating excitotoxicity. It is also possible that existing high levels of glutamate release in the context of aSyn expression resulted in activation of signaling cascades in this study, which altered glutamate receptor expression, and DNAm changes occurred after this expression change. Although it is not possible to determine which of these scenarios is true based on our data, this study provided reasonable support for a role of aSyn-mediated epigenetic deregulation of glutamate receptor expression, with potential downstream impacts on synaptic activity and excitotoxicity. Indeed, preclinical studies have highlighted glutamate receptors as therapeutic targets for PD, and regulation of glutamate signaling genes may be influenced by lifestyle factors, such as dietary exposure to neurotoxins or tea polyphenols (26,(68)(69)(70). Consistent with our findings, aSyn overexpression was also recently shown to induce glutamate release from mouse neurons and astrocytes in vitro and in vivo and to activate glutamate receptors (71). While we did not directly examine glutamatergic transmission in this study, our results provide evidence from an epigenetic perspective to expand existing research on the role of glutamate signaling in the etiology and prevention of PD, and support a potential role for aSyn in this pathway.
In addition to these core similarities in the pathways affected by each aSyn variant, WT and A30P aSyn also had unique effects on the DNA methylome and hydroxymethylome. Broadly, many more DNAm and DNAhm changes were seen in WT aSyn cells than in A30P cells in comparison with controls. This was somewhat surprising in light of our previous work, which showed greater effects at the transcriptional level in A30P cells (9). These differences were unlikely to be related to aSyn protein levels, as both WT and A30P LUHMES neurons expressed similar levels of aSyn (9). It was also unexpected that there were almost no changes to non-CpG methylation with aSyn overexpression (one differentially methylated CH probe in WT vs. A30P aSyn cells) (Supplementary Material, Table S3) despite known enrichment of non-CpG methylation in neurons (72,73).
Several mechanisms may explain the greater number of DNAm changes observed in WT aSyn cells than in A30P aSyn cells. First, DNA damage occurs to a greater extent in WT aSyn LUHMES neurons than in A30P aSyn LUHMES neurons, and it is possible that some DNAm changes in WT aSyn cells are a ref lection of this damage and/or subsequent repair attempts (9,74). Following DNA damage, protective mechanisms in areas of active transcription maintain DNA repair and restore epigenetic modifications, which may also explain why we observed more DNAm changes in intergenic regions but fewer in intragenic regions than expected by chance in both genotypes (74,75). Second, differences in characteristics of each protein could also explain the greater number of DNAm alterations in WT cells. WT aSyn could influence DNAm through various mechanisms, including binding to membrane-bound G protein-coupled receptors and initiation of signaling cascades; facilitation of SNARE complex assembly at presynaptic membranes, which could allow neurotransmitters to activate signaling cascades at postsynaptic neurons; and binding to DNA (6,8,76). The A30P mutation prevents aSyn from binding to membranes, thus removing some of these avenues of influence (7). Finally, WT aSyn can sequester Dnmt1 from the nucleus in mice, resulting in global loss of DNAm in the brain (77). Differences in DNMT1 sequestration between aSyn variants could reduce the capacity for maintenance of DNAm in a nonspecific manner. aSyn aggregation, which is more likely with the WT protein, has also been shown to increase this sequestration (77). Although A30P aSyn is more likely to be localized to the nucleus than WT aSyn, the ability of A30P aSyn to sequester Dnmt1 has not been assessed (76)(77)(78). In addition to inf luencing DNAm loss in WT aSyn cells and gain in A30P aSyn cells, this potential difference in nuclear aSyn may also affect the amount of DNAm at transcription-associated regions due to its ability to bind DNA (6,(76)(77)(78)(79).
In contrast to DNAm, the patterns of differential hydroxymethylation were similar between aSyn genotypes, and fewer DNAhm changes were seen at our statistical thresholds. This was expected due to the fetal origin of the LUHMES cell line (43,80). Some of these DNAhm patterns in differentiated LUHMES cells may still be informative for later life stages as certain genes, particularly cell type-specific genes, retain similar brain DNAhm signatures through early development into adulthood (72). The overall increase in DNAhm with aSyn overexpression observed here was consistent with previous reports showing that DNAhm levels are increased in neurons and cerebellar white matter from PD patients (48,49). The opposite direction of change for DNAm and DNAhm at the same loci was also expected, as DNAm must be oxidized for the formation of DNAhm (60). From the results of this study, it is not possible to determine whether this DNAhm was a transient intermediate or represented a stable epigenetic mark; however, correlating DNAhm patterns with gene expression levels has been shown to provide some insight regarding loci where it is associated with transcription (43).
It is important to understand the interplay between DNAm, DNAhm, and gene expression to gain insight into the molecular etiology of PD, and cell culture models, such as LUHMES cells, provide appropriate platforms for concurrent multi-omics profiling. The relatively weak correlations between DNA(h)m changes and expression changes in this study (< 0.1% of genes in either genotype) were expected and were consistent with the literature, where correlations with expression at only 0.6%-15% of CpG sites in blood and up to 0.3% of CpG sites in brain have been reported (81)(82)(83)(84)(85). Therefore, it was also unsurprising that there was no significant enrichment for the overlap between differentially methylated, differentially hydroxymethylated, and differentially expressed genes. However, the few loci where transcriptional changes were correlated with DNA(h)m changes may still have either functional relevance for PD or represent potential biomarkers of neurodegeneration-relevant pathways. For the subset of genes with changes to both DNAm and mRNA expression, DNAm changes may have occurred first, thus impacting expression of disease-relevant genes, or the reverse may be true, i.e. DNAm changes may reflect alterations in gene expression (81).
When changes in DNAm precede changes in transcription, DNAm may represent a target for modulation of gene expression to prevent, for example, upregulation of glutamate signaling genes in individuals with aSyn multiplications. Indeed, an epigenetic editing approach using site-directed DNAm via a dCas9-DNMT3A fusion construct was reported previously to successfully reduce aSyn levels in SNCA-triplication patient-derived iPSC dopaminergic neurons, demonstrating therapeutic possibilities for manipulation of DNAm in modulating aberrant transcription (86). Conversely, when DNAm changes follow transcriptional changes, DNAm is more likely to represent a biomarker of expression.
Finally, DNAm alterations that did not affect gene expression at this cross-sectional time point may have an impact on expression at later times or upon exposure to particular stimuli (87). Further studies involving exposure of aSyn-expressing neurons to substances and stimuli relevant to PD risk and/or neuroprotection may be useful to examine whether the role of DNA(h)m in the toxicity of aSyn is related to an inf luence on the ability of the cell to respond appropriately to either risk or protective factors for PD.
The assessment of whether changes in DNA(h)m associated with aSyn overexpression are correlated with transcriptional changes will provide insights regarding whether the roles of cytosine modifications in the toxicity of aSyn are mediated by effects on gene expression. However, it should be noted that this study was not designed to elucidate the mechanisms underlying differential methylation in certain loci in our aSyn-expressing LUHMES neurons. In the general context of the complex and multi-layered regulation of chromatin biology, and given the specific limitations of existing epigenome-wide association study approaches, it is only possible to speculate about the pathways and signals targeting specific CpG loci across the human genome, including the observed alterations in DNAm levels at the endogenous SNCA gene itself (33).
This study provided novel insights into the epigenomic impact of aSyn overexpression and set the stage for future studies in samples derived from PD patients. Nevertheless, it should be taken into consideration that integration of the SNCA transgene occurred at random positions in the genome in our cells, which could affect DNAm levels surrounding the integration positions. To minimize this variability and represent a range of integration positions, we used eight replicates of each LUHMES cell genotype and seven replicates of control LUHMES cells for EPIC array profiling, with each replicate consisting of a pool of thousands of cells. Cell culture has also been shown to inf luence DNAm (88)(89)(90). We reduced the effects of such variability by using replicates obtained from a range of passages and culture dishes, and regressed out the effects of passage before analysis of differential DNA(h)m. In addition, as expected, determination of DNAhm using a combination of oxidative bisulfite-and bisulfite-converted samples introduced technical noise, which prompted us to reduce the number of sites for DNAhm analysis (91). Finally, the LUHMES cells used in this study are an approximation of dopaminergic neurons in the brains of PD patients, but are subject to the inherent limitations of any in vitro cellular model. However, LUHMES cells and other tissue culture models of aSyn overexpression (primarily iPSCs and neuroblastoma lines) have been used successfully in previous studies to uncover biochemical mechanisms of PD, including mitochondrial dysfunction, increased H3 acetylation, and cell-cell transfer of aSyn protein (9,92,93). LUHMES cells provide a useful platform in which to study the biology of aSyn and possibly gain insight into preclinical stages of PD, prior to aSyn aggregation and neurotoxicity on a large scale. As aSyn overexpression has been shown previously to affect cellular physiology even at the neural precursor stage, such studies are still valuable (94). In addition, our LUHMES cells have an advantage over iPSC studies in that the cells are of neural origin, reducing the potential for epigenetic changes associated with iPSC reprogramming (89,90).
Overall, this study demonstrated that WT and A30P aSyn overexpression each had significant impacts on the DNA methylome of human dopaminergic neurons, and that only a small proportion of these DNAm changes occur concomitantly with DNAhm and changes in expression of neurodegeneration-related genes. This study contributes to our understanding of the molecular genetic etiology of PD and lays groundwork for future studies to investigate whether epigenetic and transcriptional changes associated with increased expression or mutations in aSyn may be reversible by modification of environmental and/or lifestyle factors.

Generation of aSyn-expressing LUHMES cells
Expression of aSyn in LUHMES cells was achieved as described previously (9,95); the methods below are reproduced from (9) with permission. Full-length human aSyn c-DNA (SNCA, NM_000345) and A30P aSyn were subcloned into the second-generation bicistronic lentiviral vector, pWPI (Tronolab, Lausanne, Switzerland), under the control of the chicken β-actin (CBA) promoter. All cloned sequences were verified by direct sequencing (StarSeq, Mainz, Germany). The pWPI lentiviral vector containing only the internal ribosomal entry site (IRES)green f luorescent protein (GFP) cassette was used as a positive control for infection in all experiments. Second-generation lentiviral particles were generated as described previously (96). Brief ly, 293 T cells were transiently co-transfected with the modified transfer vector plasmids and second-generation packing system (Tronolab). Supernatant was collected 48-h posttransfection, concentrated by PEG-it Virus Precipitation Solution (System Biosciences, Palo Alto, CA, USA), and resuspended in Panserin 402 (Pan-Biotech, Aidenbach, Germany). Lentiviral vector titers were measured quantitatively by qPCR using primer sequences specific for the woodchuck hepatitis virus post-transcriptional regulatory element (WPRE), as described previously (97). To generate cells expressing aSyn, proliferating LUHMES cells were infected with equimolar amounts of IRES-GFP and WT aSyn-IRES-GFP or A30P aSyn-IRES-GFP lentivirus. Cells positive for green f luorescence were selected by FACS (FACSCalibur; BD Biosciences, Franklin Lakes, NJ, USA).

Immunocytochemistry
Immunocytochemistry was performed as described previously (9,95); the methods below are reproduced from (9) with permission. Cells were grown on coverslips and fixed at different time points of differentiation with 4% paraformaldehyde for 40 min at room temperature (RT). After washing with 1× phosphate buffered saline (PBS), the cells were permeabilized with 0.5% Triton/PBS for 15 min, and finally blocked with 1.5% normal goat serum (NGS)/PBS (blocking buffer) for 1 h at RT. Cells were incubated with anti-aSyn antibody (mouse, 1:1000; BD Transduction Laboratories, Franklin Lakes, NJ, USA) and anti-TH antibody (rabbit, 1:1000; Millipore, Billerica, MA, USA) for 3 h. The cells were then washed with 1× PBS and probed for 2 h at RT with secondary antibodies (rabbit or mouse, 1:1000 Alexa 488 and 1:1000 Alexa 555; Life Technologies, Carlsbad, CA, USA) diluted in blocking buffer. Finally, the nuclei were stained with DAPI and mounted with Mowiol (Sigma-Aldrich). The images were acquired on a Leica 6000B microscope (Leica Microsystems, Wetzlar, Germany).

DNA extraction and bisulfite conversion
Cell pellets were thawed on dry ice and homogenized with a 20G needle in Qiagen Buffer RLT Plus (Qiagen, Valencia, CA, USA) with β-mercaptoethanol. DNA extractions were performed with a Qiagen AllPrep DNA/RNA Mini Kit in accordance with the manufacturer's instructions. DNA quantity and purity were assessed by spectrophotometry.
DNA samples of 1.5 μg were split into two equal aliquots of 750 ng for sodium bisulfite (BS) conversion or oxidative bisulfite (oxBS) conversion using the NuGEN TrueMethyl oxBS Module (NuGEN Technologies, San Carlos, CA, USA). One aliquot per sample was treated with the oxidation protocol, while the second aliquot was subjected to mock oxidation. All samples were then sodium bisulfite converted according to the manufacturer's instructions.

DNA methylation and DNA hydroxymethylation analysis
Aliquots of 750 ng per sample of BS-or oxBS-converted DNA were run on Infinium HumanMethylationEPIC (EPIC) BeadChips (Illumina, San Diego, CA, USA) according to the manufacturer's instructions, producing data for 853 307 CpG and 2880 CNG sites (where N represents any other nucleotide) (data available at https://www. ncbi.nlm.nih.gov/geo/ with accession No. GSE181126). Beta values were generated from raw intensity signals using GenomeStudio software (Illumina) and exported into R ver. 3.6 (R Foundation for Statistical Computing, Vienna, Austria) for data analysis. Filtering was performed to remove 59 internal SNP control probes and 43 254 cross-hybridizing probes (98). The pfilter function in the wateRmelon package was additionally used to remove sites with a detection p > 0.05 in 1% of samples (7685) and sites with a bead count < 3 in 5% of samples (1626) (98). This left 813 589 EPIC probes in the final dataset. The dasen function in wateRmelon was used to normalize the oxBS and BS data separately (99). Batch effects from chip, position, and passage were removed using the ComBat function in the SVA package (100).
The hmC values were calculated by subtracting the normalized oxBS beta values from the normalized BS beta values. A threshold for hmC detectability (3.6%) was calculated using the 95% quantile of negative hmC values generated after subtraction; all probes with a mean hmC level below this threshold were discarded (91).
Differential methylation between control (n = 7), WT aSyn (n = 8), and A30P aSyn (n = 8) cells was calculated by linear regression using limma with the model DNA(h)m ∼ Genotype (101). Individual sites were considered significant at FDR ≤ 0.05 and |delta beta| ≥ 0.05 (calculated by subtracting the mean beta value across all replicates of the group in question from the reference group, e.g. WT aSyn beta − Control beta). The skewness.norm.test function from the normtest R package was used to assess delta beta skewness, with 1000 Monte Carlo simulations (https://cran.r-project.org/web/ packages/normtest/normtest.pdf).
CMRs were defined separately in pre-processed DNAm and DNAhm data from LUHMES neurons using the cmr function from the CoMeBack R package with a maximum distance of 1 kb, minimum Pearson's correlation coefficient of 0.4, and minimum of two CpGs per CMR (58). Composite betas were calculated with the cmr_comp function, using the default method (weighted average of the first principal component of CMR probe DNA(h)m levels). As in site-specific analysis, linear modeling with limma was applied to the DNAm and DNAhm composite betas, with one composite beta value and one CpG probe ID representing each CMR. The overlap in differentially (hydroxy)methylated probes from the CMR and site-specific analyses was permuted 10 000 times by randomly sampling a set of CpG probes from the preprocessed EPIC DNAm or DNAhm data the same size as the number of significant site-specific probes, randomly sampling a set of representative CMR probes from the DNAm or DNAhm CMRs the same size as the number of significant CMRs, and calculating how many iterations had greater or smaller overlap than the real value.
All code is publicly available at github.com/samschaf/ LUHMES. Pyrosequencing 1 μL of BS or oxBS-converted DNA per 30-μL reaction was used to perform PCR with HotStarTaq DNA Polymerase (Qiagen) for 45 cycles with an annealing temperature of 58 • C. PCR samples were visualized on a 1% agarose gel post-amplification to ensure integrity.
Pyrosequencing assays were designed with PyroMark ® Assay Design 2.0 software (Qiagen) (Supplementary Material, Table S5). 10 μL of PCR reaction mixture, 70 μL of binding buffer, and 12 μL of sequencing reaction mix per sample were prepared with PyroMark ® Gold Q96 Reagents (Qiagen) according to the manufacturer's instructions. Sequencing was performed with the PyroMark ® Q96 ID (Qiagen). Differential DNAm/DNAhm between groups was assessed by Welch's two-sample t test (two-tailed), using the t.test() function in R.

Generation of aSyn-expressing mice and genome-wide DNAm/DNAhm profiling
Generation of WT aSyn-expressing mice and DNA extraction were performed as described previously (14). Hippocampal DNA from four WT and four transgenic mice raised under standard housing conditions to the age of 12 months was profiled by RRBS for genome-wide DNAm and DNAhm. RRBS libraries were prepared following the protocol of Gu et al. (102), using 400 ng of genomic DNA per sample. 0.4 ng of unmethylated phage lambda DNA per sample was added during the MspI digestion step as a control for bisulfite conversion efficiency, and MspI reactions were performed with incubation at 37 • C for 16 h. Digested DNA was end-repaired, and then purified using AMPure XP magnetic beads (Beckman Coulter Life Sciences, Fullerton, CA, USA) at 2× concentration, and 0.1 μM methylated adapters (IDT, Coralville, IA, USA) were used in ligation reactions. Ligated DNA was purified with 1.5× magnetic beads and quantified using a Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). Eight samples per lane were pooled in equimolar fashion, and each pool was split into two aliquots for oxBS and BS conversion using the NuGEN TrueMethyl oxBS Module (NuGEN Technologies). One aliquot was subjected to oxidation, while the second aliquot was subjected to mock oxidation; the bisulfite conversion protocol was then applied to both. Libraries were then amplified by PCR with Pfu Turbo Cx HotStart DNA polymerase for 14 cycles (Agilent Technologies, Santa Clara, CA, USA). Size selection (200-550 bp) was performed using a 0.55× AMPure beads wash (Agencourt Bioscience, Beverly, MA, USA) followed by two 1.2× washes for cleanup. The final libraries were eluted in 25 μL of Qiagen EB Buffer (Qiagen), and quality was assessed by PCR, spectrophotometry, fluorometry, and DNA high-sensitivity chip (Agilent Technologies). 75 bp paired-end sequencing was conducted with Illumina HiSeq 2500 by the Michael Smith Genome Sciences Centre (Vancouver, bc, Canada).
BSMAP was used to align reads to the mm10 and phage lambda genomes and to calculate methylation ratios (103). Sequence data was assessed for quality before and after alignment using FastQC software (https:// www.bioinformatics.babraham.ac.uk/projects/fastqc/). Methylation ratio data was then read into R ver. 3.6.3 and filtered for: a minimum read depth of 10×; a maximum read depth of the 99.9% quantile of coverage; and sites covered in all samples. oxBS beta values were used to represent DNAm, while BS − oxBS values were used to represent DNAhm, with the 95% quantile of negative hmC values generated after subtraction used as a detectability threshold (91). We removed cytosines with ≤ 10% variability between groups to reduce the multiple testing burden for differential DNA(h)m analysis.
Differential DNA(h)m analysis was conducted with the betaRegression function from the BiSeq R package, using beta-binomial regression with the Wald test for significance (104). The p-values were adjusted for multiple comparisons using the Benjamini-Hochberg method. Cytosines were considered differentially methylated at padj ≤ 0.05 and |delta beta| ≥ 0.1.
To determine whether more genes differentially methylated in WT aSyn-expressing LUHMES neurons were also differentially methylated in mice than expected by chance, we first subset the RRBS data to genes differentially methylated in LUHMES cells. We considered the 'hit list' of differentially methylated cytosines in mice as the number of cytosines within this subset passing padj ≤ 0.05 and |delta beta| ≥ 0.1. Next, we randomly selected a set of cytosines the same size as this 'hit list' from the subsetted RRBS background, and determined the genes to which they mapped. The number of differentially methylated genes in the random selection was compared to the number of differentially methylated genes in the 'hit list.' This process was repeated 10 000 times in each of the DNAm and DNAhm datasets. Enrichment and depletion p-values were calculated by dividing the number of iterations with more or fewer differentially methylated genes than the real number by the number of permutations.

Gene ontology enrichment
Genes were assigned to CpGs if they belonged to the longest UCSC RefGene transcript annotated to each site, resulting in 30 325 unique coding and non-coding gene annotations for the DNAm dataset and 23 190 unique gene annotations for the DNAhm dataset. Gene Ontology enrichment analysis was performed with the over-representation (ORA) method in ermineR 1.0.1.9, using differentially methylated genes identified with the effect size and significance thresholds described above as input (59). Genes that were differentially methylated or hydroxymethylated in both WT aSyn and A30P aSyn cells in either direction were removed from the input list and analyzed separately. All code is publicly available at github.com/samschaf/LUHMES.

RNA sequencing and differential expression analysis
RNA-seq of LUHMES cells used for multi-omic integration analysis in this study was described previously (9,95) (data available at https://www.ncbi.nlm.nih.gov/ geo/ with accession No. GSE89115); the methods below are reproduced from (9) with permission. Total RNA from differentiated LUHMES cells was extracted and purified using an RNeasy mini kit (Qiagen) in accordance with the manufacturer's instructions. Three biological replicates were used for RNA-seq experiments. Sequencing and RNA quality analyses were performed as described previously (105). Brief ly, RNA quality was assessed using RNA 6000 Nano chips run on a 2100 Bioanalyzer (Agilent Technologies). Libraries were prepared using a TruSeq RNA Sample Preparation v2 kit (Illumina). The library quality was checked using High-Sensitivity DNA chips on a Bioanalyzer (Agilent Technologies). The sample concentration was measured with a Qubit dsDNA HS Assay Kit and adjusted to 2 nM before sequencing (50 bp) on a HiSeq 2000 sequencing system (Illumina) in accordance with the manufacturer's instructions.
Differential gene expression analysis was performed as described previously (105). Brief ly, RNA-seq data were aligned to the genome using STAR with the default options, which generated mapping files (BAM format). Differential expression read counts were generated using featureCounts, and samples were compared for differential expression using DESeq2 (106,107). Genes with p ≤ 0.05 and mean read count ≥10 were considered to be differentially expressed.

Chromatin immunoprecipitation sequencing (ChIP-seq)
Chromatin and DNA were cross-linked by adding formaldehyde to the culture dish to a final concentration of 1% and incubating for 10 min at 25 • C. The cross-linking reaction was quenched with glycine at a final concentration of 0.125 M at 25 • C. Cells were then collected, washed with cold PBS, and pelleted by centrifugation at 900 rcf for 5 min at 4 • C. Cell pellets were resuspended in RIPA SDS buffer (140 mM NaCl, 1 mM EDTA at pH 8.0, 1% Triton X-100, 0.1% sodium deoxycholate, 10 mM Tris-Cl at pH 8.0, 1% SDS) supplemented with Roche complete protease inhibitors, incubated at 4 • C for 10 min, and then used for shearing in Covaris ® instrument. The sheared chromatin was cleared by centrifugation at 16000 rcf for 5 min and an aliquot was used to check DNA size on an Agilent 2100 Bioanalyzer (Agilent Technologies). Immunoprecipitation was performed when DNA was sheared to fragments of 150-300 bp.

Multi-omic integration
The overlap between differentially methylated, hydroxymethylated, and expressed genes in each comparison was permuted 10 000 times to test whether it was greater or smaller than expected by chance. For this purpose, a set of EPIC probes from the preprocessed DNAm data the same size as the number of differentially methylated probes was randomly selected, a set of EPIC probes from the preprocessed DNAhm data the same size as the number of differentially hydroxymethylated probes was randomly selected, and a set of genes from the RNAseq data the same size as the number of differentially expressed genes was randomly selected. The gene-level overlap between randomly sampled DNAm, DNAhm, and RNA-seq data was then calculated, and this number was compared to the real number to calculate enrichment and depletion p-values.
DNAm, DNAhm, and RNA-seq data were also integrated using SMITE in R ver. 3.6, with adjusted p-values as significance input, delta betas as effect size input for DNAm and DNAhm, and log2FC values as effect size input for RNA-seq (62). Genomic regions were split into enhancers (± 5000 bp from H3K4me1 ChIPseq peaks), promoters (± 1000 bp from TSS), and gene bodies (remaining regions). DNAm and DNAhm p-values were weighted across genomic regions according to their significance using Stouffer's method, creating a combined p-value for each region. The p-values from all modifications were then used to create weighted genelevel scores, with application of the following weights: expression, 0.4; enhancer DNAm, 0.125; promoter DNAm, 0.125; body DNAm, 0.1; enhancer DNAhm, 0.125; promoter DNAhm, 0.125; body DNAhm, 0.1. Weights were chosen to return modules likely to contain differentially expressed genes, with equal probability of contribution from differential DNAm and/or differential DNAhm (total weight of 0.25 for each methylation type). Enhancers and promoters were weighted slightly higher than gene bodies due to the higher likelihood that DNAm at these regions would inf luence gene expression. Alternative SMITE models with equal weights for each modification, DNAm weighted highest, DNAhm weighted highest, or DNAm and expression weighted higher than DNAhm were also tested, and produced overall similar results (Supplementary Material, Tables S11, S12). Scores were annotated to a REACTOME proteinprotein interaction network, and a spin-glass algorithm was used to identify modules that had altered DNAm, DNAhm, and expression for each comparison. All code is publicly available at github.com/samschaf/LUHMES.

Supplementary Material
Supplementary Material is available at HMG online.