Aberrant epigenome in iPSC‐derived dopaminergic neurons from Parkinson's disease patients

Abstract The epigenomic landscape of Parkinson's disease (PD) remains unknown. We performed a genomewide DNA methylation and a transcriptome studies in induced pluripotent stem cell (iPSC)‐derived dopaminergic neurons (DAn) generated by cell reprogramming of somatic skin cells from patients with monogenic LRRK2‐associated PD (L2PD) or sporadic PD (sPD), and healthy subjects. We observed extensive DNA methylation changes in PD DAn, and of RNA expression, which were common in L2PD and sPD. No significant methylation differences were present in parental skin cells, undifferentiated iPSCs nor iPSC‐derived neural cultures not‐enriched‐in‐DAn. These findings suggest the presence of molecular defects in PD somatic cells which manifest only upon differentiation into the DAn cells targeted in PD. The methylation profile from PD DAn, but not from controls, resembled that of neural cultures not‐enriched‐in‐DAn indicating a failure to fully acquire the epigenetic identity own to healthy DAn in PD. The PD‐associated hypermethylation was prominent in gene regulatory regions such as enhancers and was related to the RNA and/or protein downregulation of a network of transcription factors relevant to PD (FOXA1, NR3C1, HNF4A, and FOSL2). Using a patient‐specific iPSC‐based DAn model, our study provides the first evidence that epigenetic deregulation is associated with monogenic and sporadic PD.


Introduction
Parkinson's disease (PD) is a neurodegenerative disorder associated with the progressive loss of dopaminergic neurons (DAn) in the substantia nigra pars compacta (SNpc) (Lang & Lozano, 1998a,b). Yet PD is recognized as a systemic disease affecting other tissues apart from the nervous system (Hoepken et al, 2008;Beach et al, 2010;Shannon et al, 2012). Although most cases are sporadic, around 5-10% encompass monogenic forms caused by pathogenic mutations in PD-associated genes (Farrer, 2006). Among these, missense mutations in the leucine-rich repeat kinase 2 (LRRK2) gene are the most frequent cause of familial PD (Paisan-Ruiz et al, 2004;Zimprich et al, 2004) and also of the common sporadic form. The LRRK2 G2019S mutation alone explains up to 6% familial and 3% sporadic PD cases in Europeans (Di Fonzo et al, 2005;Gilks et al, 2005) and up to 20% of total cases among Ashkenazy Jews (Ozelius et al, 2006) or 40% in North African Berbers (Lesage et al, 2006). In addition, LRRK2-associated PD (L2PD) is clinical and neuropathologically similar to sporadic PD (sPD) lacking LRRK2 mutations (Healy et al, 2008), thus representing a valuable system to investigate the most common form of disease. Moreover, the reduced penetrance of G2019S in L2PD suggests the involvement, akin to sPD, of yet unknown disease-modifying factors (Healy et al, 2008).
Genetic and epigenetic alterations contribute to the physiopathology of diseases (Bergman & Cedar, 2013). In PD, beyond largely studied genetic defects, the epigenomic landscape of disease remains unknown (Urdinguio et al, 2009;van Heesbeen et al, 2013). Epigenetic modifications are inheritable changes of gene expression without alterations in the DNA sequence which can virtually capture the influence of environmental factors (Feil & Fraga, 2011). Thus, epigenetic alterations can reflect the relationship among factors which have been postulated to play a role in complex neurodegenerative disorders such as the individual genetic background, the environment, and the aging process. To date, epigenetic studies in central nervous system disorders have been hampered by the inaccessibility to disease targeted cells from patients, and especially of DAn from the SNpc in PD. The few published reports in PD were performed in blood cells using single gene candidate approaches (Kontopoulos et al, 2006;Pieper et al, 2008;de Boni et al, 2011;Jin et al, 2014) or human postmortem cortex and cerebellum representing end-points of disease (Desplats et al, 2011;IPDGC, 2011;Masliah et al, 2013). Results from these studies have been sometimes conflicting and did not yield robust associations. In addition, previous studies on iPSC-derived DAn did not explore in the role of the epigenome in PD Nguyen et al, 2011;Seibler et al, 2011;Cooper et al, 2012;Jiang et al, 2012;Sanchez-Danes et al, 2012b;Rakovic et al, 2013;Reinhardt et al, 2013;Ryan et al, 2013;Schondorf et al, 2014).
In the present study, we investigated the epigenome of PD by performing a genomewide DNA methylation study of CpG dinucleotides and a transcriptome study using an in vitro PD model of patient-specific disease-relevant cells (DAn). This cell system consisted in induced pluripotent stem cell (iPSC)-derived DAn generated upon cell reprogramming of parental skin cells from L2PD patients carrying the G2019S mutation (n = 4), sPD patients without LRRK2 mutations (n = 6), and gender-and age-matched healthy subjects (n = 4) (Sanchez-Danes et al, 2012a,b). Studied cell lines were similar in PD and controls as regards their properties and their maturation state and included 30-day morphologically and functionally mature ventromedial (vm)-DAn which were mostly of the A9 subtype (Sanchez-Danes et al, 2012b). The goal of our study was to explore for the first time the epigenomic landscape of PD using iPSC-derived DAn, and to compare the methylome of the common sPD form with the uniquely resembling L2PD monogenic form.

Results
Studied cell lines from PD patients and healthy controls were generated and characterized in parallel blind to researcher in a previous study (Sanchez-Danes et al, 2012b) using a 30-days differentiation protocol (Sanchez-Danes et al, 2012a) (Table 1, Fig 1 and Materials and Methods). Resulting iPSC-derived DAn had similar morphological and functional properties as well as similar full DAn maturation state in PD and controls (Fig 1) (Sanchez-Danes et al, 2012b). Yet consistently with the late onset of disease, the DAn cells from PD patients developed specific neurodegenerative phenotypes upon long-term culture (75-days) including impaired axonal outgrowth, deficient autophagic vacuole clearance, and accumulation of a-synuclein (SNCA) (Sanchez-Danes et al, 2012b;Orenstein et al, 2013).

Epigenetic changes are associated with monogenic and sporadic PD
We performed a comprehensive genomewide DNA methylation analysis of 30-days iPSC-derived DAn using the Illumina 450k methylation platform (Bibikova et al, 2011). Unsupervised hierarchical clustering of CpGs methylation values showed different DNA methylation profiles between both forms of PD (L2PD and sPD) and controls indicating robust differences between PD and controls (Figs 2A and EV1A). We further detected 1,261 differentially methylated CpG sites (DMCpGs) in L2PD and 2,512 in sPD with respect to controls under an absolute mean methylation difference above 0.25 (Bibikova et al, 2011) and an adjusted P below 0.05 (Figs 2B and EV1B,and Table EV1). Most DMCpGs in L2PD were common to sPD (78%) and no significant methylation differences were found when comparing L2PD and sPD using the same criteria mentioned above, indicating that L2PD and sPD share similar methylation profiles. Accordingly, both groups were merged for further analysis. In all PD subjects, we identified 2,087 DMCpGs as compared to controls including hypermethylation in 1,046 regions and hypomethylation in 1,041. DMCpGs mostly affected gene bodies and promoters but were also enriched at intergenic regions. Hypermethylated DMCpGs were more often located outside CpG islands, shores, or shelves (73% vs. 31% in background, P < 0.4 × 10 À14 ) ( Fig 2C). Notably, genes associated with DMCpGs were largely involved in neural functions and transcription factor (TF) activity (Table EV2). These data indicate that DAn from PD patients show epigenetic abnormalities. They also indicate that monogenic L2PD and the sporadic form of PD share similar DNA methylation changes.

Epigenetic changes manifest only in DAn from PD patients
We further investigated whether those DNA methylation changes observed in DAn were already present in fibroblasts or in undifferentiated iPSCs from the same subjects by analyzing a subset of representative individuals (two L2PD, two sPD, and three controls). Isogenic fibroblasts and iPSCs from these subjects showed no methylation differences between PD and controls neither for the 2,087 DMCpGs detected in iPSC-derived DAn ( Fig 2D) nor for any other given CpG ( Fig EV1C) under an absolute mean methylation difference above 0.25 (Bibikova et al, 2011) and an adjusted P below 0.05. Yet fibroblasts, undifferentiated iPSCs and iPSC-derived DAn showed distinct DNA methylomes as expected for each specific cell type (Doi et al, 2009;Hochedlinger & Plath, 2009) (Fig EV1D). Interestingly, a total of 75% of all the 2,087 DMCpGs did not change in the differentiation from iPSCs to DAn in PD patients, whereas only 40% DMCpGs remained unchanged in controls, suggesting an incomplete epigenomic remodeling in PD in spite of the successful reprogramming and differentiation into mature DAn (Sanchez-Danes et al, 2012b) (Table EV3 and Fig EV2). These findings suggest that latent molecular defects in skin cells from PD patients become uncovered upon differentiation into DAn. This is compliant with the idea of PD as a systemic disease affecting other tissues apart from the nervous system (Hoepken et al, 2008;Beach et al, 2010;Shannon et al, 2012) including molecular defects which were previously reported in fibroblasts from L2PD or sPD patients (Papkovskaia et al, 2012;Ambrosi et al, 2014;Yakhine-Diop et al, 2014). In addition, to address the extent of the methylation changes attributed to DAn alone in the context of cell heterogeneity inherent to iPSC-derived DAn systems, we differentiated iPSCs from the same subjects into neural cultures not-enriched-in-DAn by omitting the lentiviral-mediated forced expression of the ventral midbrain determinant LMX1A and the supplementation with DAn patterning factors. This protocol results in a > 4-fold depletion in the number of DAn (See Materials and Methods and (Sanchez-Danes et al, 2012a)). In neural cultures not-enriched-in-DAn, we did not find methylation differences between PD and controls under an absolute mean methylation difference above 0.25 (Bibikova et al, 2011) and an adjusted P below 0.05 ( Fig EV1C). Moreover, the methylation profile from PD DAn was closer to that from neural cultures notenriched-in-DAn as compared to control DAn (Fig 2D). For any given comparison and using the same restrictive cutoffs mentioned above, we found that the overall methylation variability referred to all samples was attributable, in decreasing order, to (i) the different cell types as expected, (ii) the condition health/disease only in iPSC-derived DAn, and (iii) inter-individual differences in a relative lesser extent. Altogether, these results indicate that the identified PD epigenetic changes are specific for DAn cells and consist in the failure of PD DAn to fully acquire the mature epigenetic identity own to healthy DAn.

Gene and protein expression changes occur along with methylation changes
We next studied the transcriptome of iPSC-derived DAn by using gene expression microarrays. Compliant to the methylome analysis, no differentially expressed gene (DEG) was found between L2PD and sPD, being most of DEGS in L2PD (93%) shared in sPD. As compared to controls, we identified 437 DEGs in the PD group as a whole under an adjusted P below 0.05 (Fig 3A and B, and Table EV4). These findings are in line with two previous studies reporting expression changes associated with PD in DAn, at least with L2PD (sPD not studied) Reinhardt et al, 2013). Upregulated DEGs (n = 254) were largely involved in neural functions and transcription factor regulatory activity, whereas downregulated DEGs (n = 183) affected basic homeostasis (Table EV5). One upregulated DEG was SNCA (> 2.5-fold), a gene involved in familial PD and sPD whose encoded protein, a-synuclein, aggregates in Lewy body inclusions which represent a hallmark of PD (Lang & Lozano, 1998a,b). Another upregulated DEG was SYT11 (> 5-fold) which has been top-linked associated to PD across genomewide association studies (Nalls et al, 2014). We subsequently selected seven genes involved in neural functions from the top-30 list of upregulated DEGs (OTX2, PAX6, ZIC1, SYT11, DCT, DCC, and NEFL)  A-G After iPSC-based reprogramming of fibroblasts, we generated DAn from PD patients (n = 10) and gender-and age-matched healthy controls (n = 4) using a 30-day differentiation protocol. Blind to researcher, resulting iPSC-derived DAn showed similar properties and maturation state in PD and controls. Specifically, studied iPSC-derived DAn encompassed 30-day ventromedial (vm)-DAn-enriched cultures of morphologically mature DAn showing typical bipolar morphology, lacking PD neural phenotypes, mostly of the A9 subtype, and showing properties of tyrosine hydroxylase (TH) + mature neurons. Representative immunocytochemical analyses of vmDAn from PD patients and controls showed that (A) TUJ1-positive cells expressing neural tubulin III (TUJ1) co-expressed TH (Scale bar, 25 lm); (B) coexpressed the A9 subtype marker of inward rectifier K + channel GIRK2 (Scale bar, 25 lm); (C) co-expressed the A9 subtype marker of transcription factor Forkhead Box A2 (FOXA2) (Scale bar, 12.5 lm); and (F) co-expressed the synaptic vesicle protein synaptophysin (SYP) which was located along neurites (Scale bar, 5 lm).
(D) The 1-week neural progenitor cells (NPCs) from which iPSC-derived DAn were generated strongly expressed DAn progenitor markers such as the nuclear-related receptor 1 NURR1 (Scale bar, 12.5 lm). (E) We observed similar and comparable cell counts of total neurons in PD and controls as quantified by TUJ-1/DAPI and of DAn as quantified by TH/TUJ1. Selected DAn cells were generated upon 2-6 iPSC lines per subject (See methods). Error bars indicate SEM. C, control; TUJ1, neuronspecific class III b-tubulin. (G) Calcium imaging assay revealed strong spontaneous neuronal activity. %DF/F indicates the relative change in fluorescence of the monitored neurons. Collectively, these results indicate that the iPSC-derived DAn used our study were sufficiently mature for the studied parameters to be functional and to spontaneously form active neural networks in a similar and comparable manner for PD and controls.    . This analysis shows different DNA methylation profiles between PD patients and controls, and no differences between L2PD and sPD. PC, principal component. B Bar diagram showing the number of differentially methylated CpGs (DMCpGs) among iPSC-derived DAn lines (n = 14) as detected in multiple comparisons. This analysis indicates that L2PD and sPD share similar DNA methylation changes with respect to controls. Hypermeth, hypermethylation; hypometh, hypomethylation. C Relative distribution of DMCpGs in PD iPSC-derived DAn lines (n = 10) across inter-and intragenic regions (left), across different gene-related regions (middle), and within CpG islands (CGI), CGI shores, CGI shelves, or outside CGI (right). Bkg, background (methylation platform). D DNA methylation analysis across the cell types involved in the cell reprogramming and differentiation processes and also iPSC-derived neural cultures not-enrichedin-DAn. Scheme of the cell reprogramming and differentiation protocols (upper part of the panel). Heatmaps and density color code of the 2,087 DMCpGs identified in iPSC-derived DAn from PD patients with respect to controls (n = 14), as well as their methylation status in parental fibroblasts (n = 9), undifferentiated iPSCs (n = 9), and iPSC-derived neurons (not-DAn) (n = 9) (lower part of the panel) (Wilcoxon rank test for independent samples with delta-beta above |0.25| and FDR-adjusted P < 0.05). A Heatmap showing the 437 differentially expressed genes (DEGs) detected in PD patients with either L2PD (L2) or sPD (s) as compared to controls (C) (n = 14), and density color code for gene expression levels. This analysis indicates that L2PD and sPD share similar gene expression changes with respect to controls. A total of 254 DEGs were upregulated in PD whereas 183 showed significant downregulation (Linear model with empirical bayes moderation of the variance similar to ANOVA with FDR-adjusted P < 0.05). B Venn diagram representing the number of DEGs in iPSC-derived DAn from PD patients with respect to healthy subjects (n = 14) showing that most DMCpGs in L2PD were common to sPD. C Technical validation of 10 top DEGs in PD iPSC-derived DAn as compared to controls (n = 14) identified by genomewide gene expression study (upper heat-map) and validated by real-time qPCR (lower heat-map). Relative gene expression and statistical significances were calculated using the DDCt method and GAPDH, ACTB, and PPIA as endogenous controls. Pearson correlation analysis between microarray and real-time qPCR data showed a high degree of correlation (r 2 = 0.99, P = 7.7 × 10 À7 ) (bottom graph). L2, LRRK2-associated PD; s, sporadic PD. and SNCA to validate the array data by real-time qPCR ( Fig 3C) and to study their protein expression levels by immunoblot. We detected a > 2-fold protein upregulation of all genes except NEFL ( Fig EV3A). Moreover, protein expression of some DEGs co-localized at the single-cell level with the DAn marker tyrosine hydroxylase ( Fig EV3C). These findings point toward the presence of gene and also protein expression changes in DAn from PD patients which occur simultaneously along with DNA methylation changes.
PD DNA mehtylation changes are associated with gene expression We then analyzed the relationship between gene expression and DNA methylation levels in iPSC-derived DAn from PD patients. We found a significant correlation in 17% of the 2,087 DMCpGs (n = 353) involving 239 different genes (Fig 4A and B). The percentage of correlation between methylation and expression was similar to that reported in other studies (Kulis et al, 2012). Of the 353 unique correlating DMCpGs, a total of 73% showed an inverse association of DNA methylation and expression, whereas the remaining 27% DMCpGs were positively associated (Table EV6). DMCpGs showing inverse association were situated in introns and 5 0 untranslated regions, whereas those with positive association were enriched in introns ( Fig 4C). Globally, DMCpGs showing association with gene expression were more frequently located in introns (49%) than in 5 0 regions (24%). These findings suggest that, as reported in cell differentiation and cancer (Jones, 2012;Kulis et al, 2012Kulis et al, , 2013, PD gene-body DNA methylation changes play a role in regulating gene expression, both inversely and positively.

Differential DNA methylation in PD is enriched in enhancer elements
We further studied the potential molecular mechanisms underlying the association between DNA methylation and gene expression. We analyzed our data in the context of recently available functional chromatin states (Ernst et al, 2011). Hypermethylated DMCpGs were highly enriched in enhancer elements (35% vs. 12% in the background, P = 0.3 × 10 À14 ), while hypomethylated DMCpGs were enriched in Polycomb-repressed regions (31% vs. 12% in the background, P = 0.2 × 10 À14 ) ( Fig 4D). When considering only the DMCpGs correlating with gene expression, DMCpGs with inverse correlation were located in promoters and enhancers, whereas DMCpGs with positive correlation were associated with repressed regions (Fig 4E). These results suggest that PD-associated DNA methylation changes target functionally active sequences.

PD enhancer hypermethylation is associated with the downregulation of a TFs network
We next inquired whether the DMCpGs identified in PD patients were enriched for transcription factor (TF)-binding sites (TFBSs). To this end, we overlapped our data with TFBS clusters generated by ChIP-seq in the Encyclopedia of DNA Elements project (Dunham et al, 2012;Gerstein et al, 2012;Lee et al, 2012). We found enrichment for binding sites of 23 TFs in PD-hypermethylated DMCpGs and of only two TFs in PD-hypomethylated DMCpGs (Fig 5A and  Table EV7). Since the PD-associated hypermethylation affected most prominently to enhancer elements (35%) (Fig 4D), we focused on this chromatin state and observed that 65% of all PD hypermethylated enhancers became demethylated from iPSCs to DAn in controls, whereas only 8% significantly lost methylation in PD patients which overall retained higher methylation levels (Table EV3 and Fig 2D). Furthermore, among the 23 TFs showing TFBS enrichment, we found reduced gene expression of four TFs, namely FOXA1, NR3C1, HNF4A, and FOSL2, whose downregulation was also significantly associated with increased methylation levels at enhancer DMCpGs in PD (Fig 5B-D). Moreover, variable expression of these TFs together with the expression of other TFs was coordinated in PD DAn (Fig 6), suggesting that the expression level of a TF network, rather than that of individual TFs, was associated with the hypermethylation of enhancers in PD. In addition, HNF4A and FOSL2 showed a significant downregulation of protein levels as detected by immunoblot, whereas NR3C1 showed a downregulation trend which did not reach significance (Figs 7A and EV3B, and Source data for Fig 7). These data complement recent work linking TF binding to enhancers and tissue-specific hypomethylation (Stadler et al, 2011;Hon et al, 2013;Xie et al, 2013;Ziller et al, 2013), and also downregulation of TF networks to enhancer hypermethylation (Agirre et al, 2015). Our findings suggest that the incomplete epigenetic remodeling observed for the 2,087 DMCpGs identified in iPSC-derived DAn from PD patients might be mediated by the aberrant downregulation of a network of key TFs whose deficiency could prevent their target sites to become demethylated during the differentiation from iPSCs to DAn.

Discussion
We report the first genomewide DNA methylation study of iPSCderived DAn in PD. We found that DAn from PD patients are epigenetically altered with respect to healthy controls and that DAn from monogenic L2PD and sPD share similar DNA methylation abnormalities. These methylation changes cannot be explained by technical biases since all iPSC-derived DAn lines were equally generated and characterized in parallel, blind to the researcher, and differentiated DAn had similar morphological and functional properties (Sanchez-Danes et al, 2012b). We also found that the PD-associated methylation changes were not present in parental skin cells or undifferentiated iPSCs and become uncovered only upon differentiation into the DAn cells targeted in PD. In addition, the methylation profile of DAn from PD patients, with either monogenic or sporadic PD, resembled that of neural cultures not-enriched-in-DAn indicating a failure to fully acquire the epigenetic identity own to healthy DAn in PD. We also detected that gene and protein expression changes occur simultaneously along with methylation alterations in PD and that these alterations are similar in monogenic L2PD and sPD. Finally, we found that the DNA methylation changes present in PD DAn are partially associated with gene expression, target gene regulatory sequences such as enhancers, and correlate with the RNA and protein downregulation of a network of TFs. The lack of DNA methylation differences between DAn from patients with L2PD and sPD provides the proof-of-concept of common or at least converging epigenomic changes in these disease forms (Zhu et al, 2011). However, one limitation of our work is that we studied only one monogenic form of PD (L2PD) and therefore    FOSL1  ESR1  RXRA  NR3C1  JUN  GATA3  MEF2A  FOSL2  PBX3  BATF  FOXA2  CEBPB  BCL11A  EP300  TEAD4  MEF2C  JUND  ATF3  FOXA1  RAD21  HNF4A  NFIC  TCF12 OCT4 NANOG * FDR P < 0.05 * * * * * * * * * * * * * * * * * * * * * * *  (Healy et al, 2008). Yet the interpretation of our results requires further considerations. First, the epigenome is expected to reflect the influence from the environment and thus one possibility is that unknown environmental factors could trigger epigenomic changes in PD patients. This environmental contribution would be expectable for sPD (Feil & Fraga, 2011), but the agedependent reduced penetrance of the LRRK2 G2019S mutation together with its identification not only in monogenic but also in sPD cases (Healy et al, 2008) could also support an environmental involvement in L2PD (Farrer, 2006;Urdinguio et al, 2009). Second, an alternative possibility is that the genetic background from PD patients could drive DNA methylation changes not only in monogenic L2PD as expectable but also in sPD. In this regard, it has been described that sPD could be caused by an accumulation of common polygenic alleles with relatively low effect sizes (Escott-Price et al, 2015). Thus, a possible integrative interpretation of our findings is that cumulative genetic risk factors such as single nucleotide polymorphisms (SNPs) or copy number variants (CNVs) in sPD, or LRRK2 mutations in L2PD, alone or in combination with still unknown environmental factors could be underlying the common epigenetic changes detected in both disease forms. These environmental and/or genetic factors arising from different but converging pathways could ultimately cause common end-point alterations in L2PD and sPD. The PD epigenetic changes were present only in iPSC-derived DAn but not in undifferentiated iPSC, in fibroblasts, nor in iPSC-derived neural cultures not-enriched-in-DAn. These results indicate that the PD-associated DNA methylation changes are specific of DAn but not of other types since neural cultures not-enrichedin-DAn did not reveal methylation differences. They also suggest that molecular defects, of an environmental and/or genetic origin as discussed above, should be carried from the PD patient, latent in fibroblasts, and manifested only upon differentiation into DAn. Here it should be mentioned that although in PD when and where the degenerative process starts is unclear, and the related underlying   Sanchez-Danes et al, 2012b). This is in agreement with previous studies describing normal dopaminergic phenotypes in PD iPSCderived DAn Nguyen et al, 2011;Cooper et al, 2012;Rakovic et al, 2013;Reinhardt et al, 2013) with the exception of one Parkin study reporting altered dopamine release and uptake (Jiang et al, 2012). However, we found that the epigenetic changes detected in our early 30-days DAn cultures antedated late (75-days), spontaneous (not-drug induced), PD-associated phenotypes which were previously described in our model and included impaired axonal outgrowth, deficient autophagic vacuole clearance, and accumulation of a-synuclein (Sanchez-Danes et al, 2012b;Orenstein et al, 2013). Consistent with our methylation findings, these PD phenotypes were also common in DAn from L2PD and sPD patients (Sanchez-Danes et al, 2012b). Yet although the early DNA methylation changes reported here antedated these long-term PD phenotypes potential causality needs to be explored in future works. As reported in cell differentiation and cancer (Jones, 2012;Kulis et al, 2012Kulis et al, , 2013, we found that the identified epigenetic changes correlate partially with expression, indicating that PD gene-body DNA methylation changes play a role in the regulation of gene expression. We also found that the PD DNA hypermethylation changes are prominent at enhancer regulatory regions. In line with other studies (Agirre et al, 2015), we also found that in PD, the hypermethylation of enhancers appears to be related to the downregulation of a network of key TFs, rather than of individual TFs. These data complement recent work linking TF binding to enhancers and tissue-specific hypomethylation (Stadler et al, 2011;Hon et al, 2013;Xie et al, 2013;Ziller et al, 2013) and suggest a theoretic model in which the incomplete epigenetic remodeling in PD DAn might be related to the downregulation of a network of key TFs whose deficiency could prevent their binding sites to become demethylated during the differentiation from iPSCs to DAn. Key TFs from this network have been previously associated with the specification of the substantia nigra (FOXA1, NR3C1, and HNF4) and of fetal brain (FOSL2) (Ziller et al, 2013). Of these, FOXA1 is involved in maintaining the dopaminergic phenotype in adult mesodienchepalic DAn (Stott et al, 2013;Domanskyi et al, 2014), whereas the glucocorticoid receptor NR3C1 regulates DAn neurodegeneration in PD (Ros-Bernal et al, 2011). Downregulation of HNF4A has been reported in blood from PD patients as a PD biomarker correlating with motor symptoms severity (Potashkin et al, 2012;Santiago & Potashkin, 2015), whereas FOSL2 has been linked to dyskinesia, a major side effect in the DA substitutive treatment with L-DOPA (Cao et al, 2010). Since these TFs seem to be relevant to the development and maintenance of midbrain DAn in PD (Lahti et al, 2011(Lahti et al, , 2012, their deficiency might be related to impairments in the maintenance of a differentiated DAn epigenetic cellular identity in PD (Holmberg & Perlmann, 2012). In this theoretic model, the upregulation of gene and protein expression of other TFs (OTX2, PAX6, or ZIC1) and genes (SNCA, DCC, or DCT) (Figs 7B and EV3A, and Source data for Fig 7) might transiently circumvent deficits of the network of key TFs, at least in early culture stages.
One hypothesis proposes the neurodevelopmental origin of neurodegenerative diseases in the sense that molecular mechanisms occurring during development are abnormally recapitulated in neurodegenerative disease processes (Goedert et al, 1993;Schafer & Stevens, 2010). For example, a recent study showed that neural degeneration in Alzheimer's disease involves the abnormal re-activation of cellular self-destruction mechanisms which take place during neural development (Nikolaev et al, 2009). Since we found that the epigenetic pattern from PD iPSC-derived DAn showed certain similarities to neural cultures not-enriched-in-DAn, our data could be interpreted in light of this theory pointing towards the presence of possible developmental epigenetic defects associated with PD. In addition, developmental deficits of TFs which are key in the differentiation of DAn have been recently linked to PD (Laguna et al, 2015).
We showed that iPSC-derived DAn from PD patients exhibit epigenomic and transcriptomic alterations, and therefore, our findings may have implications for future cell replacement therapies. The use of pluripotent stem cells has proved as a valuable in vitro system to investigate disease, as well as a promising tool for brain transplantation and dopaminergic restoration in PD. Yet the generation of DAn for cell replacement is still matter of research (Lindvall & Bjorklund, 2011). Recently, the autologous transplantation of iPSC-derived neurons into the striatum of healthy monkeys has been suggested to be advantageous as compared to allogenic grafts, at least in terms of reduced immunogenicity (Morizane et al, 2013). Our study suggests that future cell therapeutic strategies should pay attention to the correct epigenomic status of the reprogrammed dopaminergic cells, especially when using patient own cells. In summary, using a patient-specific iPSC-based DAn system, our study provides the first evidence that epigenetic deregulation is associated with both monogenic and sporadic PD.

PD patients and generation of iPSC-derived DAn
We used mature iPSC-derived DAn lines of 30-days of differentiation generated and characterized in parallel blind to researcher previously (Sanchez-Danes et al, 2012b) using a published protocol (Sanchez-Danes et al, 2012a). Expanded subject information, cell characterization, and technical details are extensively described in these precedent studies where they should be consulted. Of these, a summary presented in Table 1, Fig 1, and here as it follows. Arm surface skin biopsies of 3 mm of diameter were obtained from LRRK2 PD patients carrying the G2019S mutation (L2PD, n = 4); sporadic PD patients lacking family history of PD and mutations in known PD genes (sPD, n = 6); and gender-and age-matched healthy individuals without neurological disease history (controls, n = 4). Primary cultures of keratinocytes or fibroblasts were reprogrammed using retroviral delivery of OCT4, KLF4, and SOX2 to generate 2-6 iPSC lines per individual (n = 50 lines). Of these, we selected the two more homogenous lines per subject which were thoroughly characterized and shown to be fully reprogrammed to pluripotency. For the directed differentiation of iPSC to ventral midbrain dopaminergic neurons (vmDAn), we used a 30-days protocol based on the lentiviral-mediated forced expression of the vmDAn determinant LMX1A together with DAn patterning factors and co-culture with mouse PA6 feeding cells to provide trophic factor support. DAn , mostly of the A9 subtype which showed properties of TH + mature neurons: (i) electrophysiological analysis showed action potentials, (ii) the majority of DAn presented the typical bipolar morphology, (iii) co-expressed the A9 subtype marker GIRK2, and (iv) co-expressed the DA transporter (DAT). Collectively, studied iPSC-derived DAn were sufficiently mature, in a similar and comparable manner for both PD and control cultures, to be functional and to spontaneously form active neural networks.

Calcium fluorescence imaging of iPSC-derived DAn
We used the cell-permeant fluorescence dye Fluo-4-AM (Life Technologies) in combination with an imaging device to monitor spontaneous activity in PD and control cells. This fluorescence probe becomes active upon binding with calcium, therefore signaling the occurrence of action potentials in the neurons. Cultures of either PD or control cells were cultured in 30-mm-diameter petri dishes for 30-days and their activity monitored as follows. Prior to imaging, the culture dish was first washed with 4 ml PBS at room temperature to remove the original culture medium. The culture was next incubated for 30 min in a solution that contained 1 ml of recording medium (RM, consisting of 128 mM NaCl, 1 mM CaCl 2 , 1 mM MgCl 2 , 45 mM sucrose, 10 mM glucose, and 0.01 M Hepes; treated to pH 7.4) and 4 lg/ml of Fluo-4. The culture was then washed with fresh RM to remove residual free Fluo-4, and finally, a volume of 4 ml of RM was left in the dish for actual recordings. The culture dish was mounted on a Zeiss inverted microscope equipped with a CMOS camera (Hamamatsu Orca Flash 2.8) and an arc lamp for fluorescence. Gray-scale images of neuronal activity were acquired at intervals of 50 ms, and in a field of view that contained on the order of 100 cells. Cultures were visualized for 30 min at room temperature. The acquired data was next analyzed to extract the fluorescence traces for each neuron. In order to compare the fluorescence traces among neurons and cultures, the raw fluorescence signal F(t) of a given neuron was normalized as F*(t) = %DF/F = 100 * (F (t) À F0)/F0, where F0 is the fluorescence signal of the neuron at rest. Spontaneous neuronal activations were revealed in the fluorescence signal as a sharp increase of typically a 10% respect to the resting state.

Generation of iPSC-derived neural cultures not-enriched-in-DAn
iPSC-derived neural cultures not-enriched-in-DAn were generated as technical controls from a subset of PD patients (n = 6) and healthy subjects (n = 3) using the previously described protocol (Sanchez-Danes et al, 2012a) but omitting the lentiviral expression of the vmDAn determinant LMX1A and DAn patterning factors. The lentiviral delivery of LMX1A results in a > 4-fold enrichment of the total final number of DAn to~30% as compared to onlỹ 7% without LMX1A (Sanchez-Danes et al, 2012a). Embryoid body (EB) formation was induced for forced aggregation and maintained in suspension in the presence of MEF condition medium for 3 days. EBs were then maintained for 10 days in ultralow attachment plates in N2B27 medium in the presence of FGF2. After that, EBs were transferred to matrigel-coated plastic chamber slides and cultured in the differentiation medium in the absence of FGF2 for 30-days. The medium for each condition was changed every other day.
DNA, RNA, and protein isolation DNA, RNA, and proteins were isolated from one million cells using the Allprep DNA-RNA-Protein kit (QIAGEN). Concentration and quality of DNA and RNA were, respectively, determined in a Nanodrop 1000 Spectrophotometer and in an Agilent 2100 Bioanalyzer. Total protein concentration was assessed with the BCA assay (Thermo Scientific).

DNA methylation analysis
The EZ DNA Methylation Kit (Zymo Research) was used for bisulfite conversion of 1,000 ng genomic DNA. Bisulfite-converted DNA was hybridized onto the Infinium Human Methylation 450K BeadChip Kit (Illumina) which covers 99% of the RefSeq genes and 96% of CpG islands at a single-base resolution (refer to Product Datasheet: http://www.illumina.com/products/methyla tion_450_beadchip_kits.ilmn). We have previously demonstrated a correlation coefficient of findings above 95% between the Illumina 450k platform and whole-genome bisulfite sequencing (Kulis et al, 2012). The Infinium methylation assay was carried out following criteria requested for this platform (Bibikova et al, 2011). All the array experiments were performed simultaneously at the same day, and PD and control samples were randomized at each BeadChip, ruling out any potential batch effect. Array data were analyzed by minfi package available through Bioconductor (Kasper Daniel Hansen and Martin Aryee. Minfi: Analyze Illumina's 450k methylation arrays. R package version 3.0.1). To exclude technical biases, we used an optimized pipeline with several filters developed at the Unidad de Hematopatología at IDIBAPS (Kulis et al, 2012). From the initial dataset of 485,512 sites (excluding probes detecting SNPs), we removed those with poor detection P-values (P > 0.01, n = 670) and those with sexspecific DNA methylation (n = 6,614). The remaining 478,228 sites were used for downstream analyses. Single CpG quantitative methylation values (beta values) in each sample were calculated as the ratio of the methylated signal intensity to the sum of methylated and unmethylated signals. Resulting quantitative methylation values ranged from 0, fully unmethylated, to 1, fully methylated CpGs. DNA from pure mouse PA6 cell cultures was included as a control sample that did not hybridize at readable levels in the human Illumina 450k methylation array, therefore not influencing methylation findings.

Genomic annotation of CpGs
The 450k Human Methylation Array data was annotated using the current hg19 version of the UCSC Genome Browser as reference sequence. We assigned the CpGs into discrete categories according to their location with respect to gene-related regions: 10-kb 5 0 region (10,000 to 1,501 bp upstream of the transcriptional start site (TSS)), 1.5-kb 5  body (from the first intron to the last exon), 3 0 UTR, 1.5-kb 3 0 region (1 to 1,500 bp downstream of the 3 0 UTR), 10-kb 3 0 region (1,501 to 10,000 bp downstream of the 3 0 UTR), and intergenic regions. Owing to the presence of alternative transcription start sites and regions containing more than one gene, some CpGs were assigned multiple gene annotations. For the DMCpGs location relative to a CpG island (CGI), we used the categories: within CGI, in CGI shore (up to 2 kb from the CGI edge), in CGI shelf (from 2 to 4 kb from the CGI edge), and outside CGI.

Statistics: identification of differentially methylated CpGs
For each CpG, we computed the difference between the DNA methylation level in the two groups under comparison. Subsequently, CpGs were considered differentially methylated (DMCpGs) when showing a delta-beta methylation difference above |0.25| and a false discovery rate (FDR)-adjusted P Wilcoxon rank test for independent samples below 0.05. A delta-beta above |0.20| is detectable with a 99% confidence in our methylation platform (Bibikova et al, 2011). Consistently, the same cutoff of delta-beta above |0.25| and FDRadjusted P < 0.05 was used in our study to identify statistically significant methylation differences (DMCpGs) between cases and controls in iPSC-derived DAn, fibroblasts, undifferentiated iPSC, and iPSC-derived neural cultures not-enriched-in-DAn for any given comparison.

Statistics: Gene expression analysis and identification of differentially expressed genes
We hybridized 100 ng of RNA onto the Genechip Human Exon 1.0 ST Array (Affymetrix) which covers > 96% of the human transcriptome. Using a robust multi-array analysis (RMA) algorithm (Irizarry et al, 2003), pre-processing of microarray data resulted in single gene log 2 -transformed values from 36,079 transcripts (22,014 Entrez RefSeq. genes). We filtered out (i) genes with group mean signals < 50 percentile of all signals, and (ii) genes with standard deviation (SD) < 50 percentile of all SDs. The remaining 4,686 Entrez RefSeq genes were used for downstream analysis. RNA from pure mouse PA6 cell cultures was included as a control sample that did not hybridize at readable levels in the Affymetrix human array, therefore not influencing expression findings. We used the Bioconductor Limma package to detect differentially expressed genes (DEGs) under a linear model with empirical bayes moderation of the variance similar to ANOVA developed to adjust for small sample size in microarray studies (Smyth, 2004). A FDR-adjusted P < 0.05 was used to identify DEGs.

Correlation of DNA methylation and expression data
Based on common gene annotations, DNA methylation data from 2,087 PD-associated DMCpGs and gene expression data from all 36,079 transcripts were overlapped. We detected 2,430 annotating pairs. Subsequently, Spearman correlation analysis was done under a FDR-adjusted P < 0.05.

Transcriptional and epigenomic characterization of differentially methylated sites
We analyzed our data in the context of functional chromatin states using the NHEK keratinocyte as background cell line (Ernst et al, 2011) which is similar to the cells used to generate iPSCs in our study. We considered regions with states 1 and 2 (designated as "active promoter" and "weak promoter") as "promoter regions"; state 3 as "poised promoters"; states 4, 5, 6, and 7 (designated as "strong enhancer" and "weak/poised enhancer") as "enhancer regions"; state 8 as "insulator"; states 9 and 10 (designated as "transcriptional transition" and "transcriptional elongation") as "strong transcription regions"; state 11 as "weak transcription"; state 12 as "Polycomb repressed"; and states 13, 14, and 15 (designated as "heterochromatin (nuclear lamina)" and "repetitive heterochromatin") as "heterochromatin regions".

Transcription factor analysis
We overlapped the methylation data from 2,087 PD DMCpGs with transcription factor (TF)-binding site (TFBS) clusters generated by ChIP-seq in the ENCODE project (Dunham et al, 2012;Gerstein et al, 2012;Lee et al, 2012) and available at the UCSC Genome Browser (http://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&g= wgEncodeHaibTfbs). For each TF, we calculated the log 2 ratio of the proportion of TFBSs in the lists of hyper-and hypomethylated DMCpGs as compared to that in the whole methylation array (

Study approval
This study has been conducted conforming the principles of the Declaration of Helsinki and the Belmont Report. The Commission on Guarantees for Donation and Use of Human Tissues and Cells of the Instituto de Salud Carlos III (ISCIII) and the local ethics committee at the Hospital Clínic de Barcelona approved the study. All subjects gave written informed consent prior to their participation in the study.

Accession number
Genomewide DNA methylation and gene expression datasets generated in this study have been deposited in the Gene Expression Omnibus (GEO) under accession GSE16453. Summaries of analyses on these crude data may be found at Tables EV1, EV3, EV4, and EV6.
Expanded View for this article is available online.

Results
Upon cell reprogramming of somatic skin cells (fibroblasts) into induced pluripotent stem cells (iPSC), we generated DAn from patients with sporadic PD (sPD) (90-95% of cases) and patients with a less frequent familial form caused by mutations in the gene LRRK2 (L2PD).
Using genomewide approaches, we found large epigenomic (DNA methylation) and gene expression changes in DAn from PD patients as compared to healthy subjects. Interestingly, these changes were largely similar in sPD and familial L2PD. In addition, the PD epigenetic changes were specific for DAn since fibroblasts, iPSCs, or other neural types (not-DAn) derived from the same PD patients did not show epigenetic abnormalities. Moreover, the PD epigenetic changes affected prominently to regulatory regions (hypermethylation of enhancers) and were related to the gene or protein downregulation of a network of transcription factors which is relevant to PD.

Impact
Using a patient-specific iPSC-derived DAn model, our findings indicate that: (i) epigenetic deregulation occurs in PD, (ii) PD-associated changes are similar in the sporadic and the L2PD familial form suggesting common etiologic processes and potential applicability of common therapeutic treatments, and (iii) future iPSC-based cell therapy strategies should pay attention to the correct epigenomic status of the reprogrammed DAn if using patient own cells.