The SH-SY5Y human neuroblastoma cell line serves as a model for neurodegenerative diseases since the cells can be differentiated to mature neurons by culturing with specific compounds. Although the neuroblastoma cells possess some basic qualities of human neurons, the active proliferation feature of the cells does not seem to be a characteristic of neurons. Since neurons are cells that are regarded as postmitotic and do not proliferate once they are terminally differentiated. Therefore, most recent studies are emphasizing the need to differentiate the neuroblastoma cells to develop homogenous neurons that exhibit morphology and biochemical characteristics of neurons. Consistent with this, we have reported that differentiated neuroblastoma cells exhibited dopaminergic neuron characteristics with increased in neurite projections, dopamine neurotransmitters, α-synuclein neuroprotein and up-regulation of dopamine receptor D2 expression when maintained in low serum retinoic acid culture medium. Our study also showed that differentiated neuroblastoma cells were more susceptible to 6-hydroxydopamine (6-OHDA) than undifferentiated neuroblastoma cells [1]. The neurotoxin 6-OHDA has significantly contributed to developing a new drug therapy for PD. It has the structural analogue of catecholamines, dopamine and noradrenaline and able to exert a potent toxic effect on catecholaminergic neurons [2]. The 6-OHDA has been frequently used in vitro experiments to induce apoptosis on neuroblastoma cells and mimic PD in a cellular model. The 6-OHDA and human neuroblastoma combo has been widely established and accepted in neuroscience research as an in vitro platform to understand a new drug’s toxicity and protective mechanism before embarking on clinical research. Neuroprotective effects of numerous plant extracts [3], flavonoids [4], vitamins [5] and phytonutrients [6] have been elucidated using the 6-OHDA and SH-SY5Y neuroblastoma cells.

Various mechanisms have been suggested in the 6-OHDA-induced toxicity of neuroblastoma cells. A recent study reported that 6-OHDA-induced toxicity resulted in degradation of Cu (copper) exporter ATP7A and Cu chaperone Atox1, leading to disruption of cellular Cu trafficking. Moreover, 6-OHDA has been found to attenuate the Cu proenzyme dopamine β-hydroxylase that converts dopamine to noradrenaline [7]. Using quantitative proteomic analysis, our research group has documented that the 6-OHDA induced apoptosis remarkably upregulated the expression of VDAC1 and 17 histone proteins via overexpression of nucleosomal degradation pathways in differentiated neuroblastoma cells [8]. Besides that, He et al., have shown that 6-OHDA impairs the autophagic flux by downregulating the transcription factor EB (TFEB), a critical component controlling lysosome biogenesis. These damaging effects of 6-OHDA significantly reduce the expression of lysosome-associated membrane protein 1 (Lamp1) and the lysosomal proteases activity [9].

However, evidence of the epigenetic changes induced by 6-OHDA in in vitro PD model remains void. Epigenetic changes can be defined as any alterations in the gene expression pattern following DNA methylation, acetylation, phosphorylation, ubiquitylations or sumolyation that do not change the sequence of DNA building blocks. The well-studied epigenetic modification is DNA methylation [10]. In the DNA methylation process, DNA methyltransferase covalently transfers a methyl group to the C-5 position of the cytosine ring in CpG dinucleotide region of DNA. The genome region with the largest number of CpG dinucleotide repeats is known as CpG islands (CGI), spanning 1000~1500 base pairs (bp) in length. CGI shores and shelves are defined as 2-kb and 4-kb long bp covering both sides of CGI. Studies have shown that the epigenetic signature in disease may spread in any CpG region, including islands, shore, shelf or open seas (outside CGI regions). Furthermore, DNA methylation may occur in various gene areas, including upstream/downstream regions of the transcription start site (TSS), 3’UTR (untranslated region), 5’UTR, exon, gene body and introns [11]. A recent study has demonstrated that hypomethylation (most cytosines do not have 5-mC compared to control) of the promoter, first exon, 5’UTR, sixth exon, N_shore and S_Shore resulted in promotion of key oncogenes, that lead to onset of cancer [12]. In regards to PD, a genome-wide study has revealed DNA methylation changes in brain tissues obtained from PD patients compared to healthy individuals. In this study, a total of 234 significant DMRs were identified in the dorsal motor nucleus of the vagus nerve, with 44 in the substantia nigra and 141 in the cingulate gyrus were discovered in PD patients [13]. Therefore, it is essential to understand if DNA methylations in the human brain also exist in the cellular model widely used in the investigation of PD. A new drug that could reverse the DNA methylation manifested in the PD model and protect neuronal cells may have a promising outcome as a potential drug for the treatment or management of PD.

In the present study, we first developed the differentiated neuroblastoma cells and subsequently exposed the cells to 6-OHDA to create a PD model. Then, the cells were subjected to Illumina Human DNA methylation EPIC BeadChip assay and data were preprocessed in Genome. Studio and analysed using R Bioconductor minfi package. We generated the CpG probes that displayed significant differences from the control group and analyzed these CpG probes using the Metascape bioinformatic tool. Our finding revealed the overrepresentation of significant PD-related CpGs that play crucial roles in neuron growth, metabolism, calcium signaling and survival mechanism.

Cell culture

The SH-SY5Y human neuroblastoma cells (Cat # CRL-2266) were cultured in a complete culture medium containing 88% of Dulbecco Modified Eagle Medium (DMEM) supplemented with 4.5 g/L glucose and L-glutamine without sodium pyruvate (Corning, Corning, NY, U.S.A), 10% heat-inactivated fetal bovine serum (FBS, Biosera, Nuaille, France), 1% Penicillin-streptomycin (P/S) (Gibco, Carlsbad, U.S.A) and 1% non-essential amino acid (NEAA, GIBCO, Carlsbad, U.S.A) 37°C in a humidified 5% CO₂ incubator. The culture medium was replaced every alternate day, and cells were sub-cultured at 70% confluence.

Differentiation of human neuroblastoma cells

The differentiated neuroblastoma cells were established by seeding 1×10⁵ cells/ml of neuroblastoma cells in a complete culture medium in a T75 flask. After 24 hours of cell seeding, the cells were exposed to a differentiation medium containing 95% DMEM, 3% heat-activated FBS, 1% P/S, 1% NEAA and 10 μM of all-trans-retinoic acid (RA, Sigma Aldrich, St. Louis, U.S.A) for 6-days. After 6-days of incubation in differentiating medium, the differentiated neuroblastoma cells expresses dopaminergic characteristics at morphological, biochemical and genetic levels [1].

Preparation and treatment protocol of 6-OHDA

On day 7 (after the 6-days of differentiation phase), the differentiated neuroblastoma cells were exposed to a culture medium containing 40 μM of 6-OHDA (Sigma Aldrich, St. Louis, U.S.A) for 24 hours. The concentration of 40 μM of 6-OHDA was selected based on the preliminary studies that demonstrated significant inhibition (IC50) on differentiated human neuroblastoma cells. The DNA extraction was performed after 24 hr treatment of differentiated neuroblastoma cells with 6-OHDA.

Extraction and purification of genomic DNA

DNA extraction was performed on untreated differentiated neuroblastoma cells (Neg Control) and differentiated neuroblastoma cells treated with 6-OHDA using the Lucigen MasterPure DNA purification kit (Thermofisher, MA, USA). Briefly, the untreated and 6-OHDA treated differentiated neuroblastoma cells were pelleted by centrifugation and the supernatant removed, leaving approximately 25 µl of liquid in each tube. The tubes were vortexed for 10 s to resuspend the cells and 300 µl of tissue and cell lysis solution containing 1 µl of Proteinase K (provided with the kit) was added to each tube. The tubes were mixed thoroughly and incubated at 65°C for 15 minutes, accompanied with a vortex mix every 5 minutes. Following this, the samples were cooled to 37°C and 1 µl of 5 µg/µl RNAase A (provided with the kit) was added to each tube. The tubes were mixed thoroughly and incubated at 37°C for 30 minutes. The tubes were placed on ice for 3~5 minutes before the DNA precipitation step.

To precipitate total DNA, 175 µl of MPC protein precipitation reagent (provided with the kit) was added to each tube containing the lysed samples. Then, the tubes were vortexed vigorously for 10 s and centrifuged (≥10,000 g) for 10 minutes at 4°C. The supernatant from each tube was transferred to a microcentrifuge tube and 500 µl of isopropanol (Sigma-Aldrich, St. Louis, USA) was added to the recovered supernatant. The tubes were inverted several times and centrifuged (≥10,000 g) for 10 minutes at 4°C to pellet total DNA in each sample. After centrifugation, the isopropanol was carefully discarded without dislodging the DNA pellet. The DNA pellet was rinsed with 70% ethanol and the residual ethanol was removed entirely from the tube. Finally, the DNA was resuspended in 35 µl of TE buffer (provided with the kit). The purity and quality of the extracted DNA were determined using a Nanoquant spectrophotometer (A260/A280 >1.8 & <2.1) and validated with agarose gel electrophoresis (no dispersion DNA fragment band). The extracted genomic DNA (gDNA) was stored at -80°C prior to use.

Genomic DNA bisulphite conversion

The purified gDNA samples were thawed to room temperature before the bisulphite conversion procedure using the EZ DNA Methylation Kit (Zymo Research, CA, USA). To initiate the bisulphite conversion step, 5 µl of M-dilution buffer (provided with the kit) was added to the 1 μg DNA sample (20 ng/μl). The samples were mixed by pipetting and incubated at 37°C for 15 minutes. Following this, 100 µl of the prepared CT conversion reagent (provided with the kit) was added to each tube and mixed thoroughly. The tubes were incubated in a thermocycler at 95°C for 30 sec, 50°C for 60 min for 16 cycles, then “hold” at 4°C. Next, 400 µl of M-binding buffer (provided with the kit) was added to Zymo-Spin IC columns (provided with the kit) and each column was placed into a collection tube. The columns were appropriately labelled. The samples were loaded into the respective Zymo-Spin IC column containing the M-Binding buffer and mixed by inverting the column. Then, the columns in the collection tubes were centrifuged (≥10,000 g for 30 seconds). The flow-through was discarded and 100 µl of M-Wash buffer (provided with the kit) was added to the columns and the columns were centrifuged (≥10,000 g for 30 seconds). Following this, 200 µl of M-Desulphonation buffer (provided with the kit) was added to each column. The columns were left at room temperature for 15~20 minutes and centrifuged (≥10,000 g for 30 seconds). The flow-through was discarded and 200 µl of M-Wash buffer was added to the column and centrifuged (≥10,000 g for 30 seconds). The wash step was repeated again. Finally, the column was placed into an appropriately labelled sterile 1.5 ml tube and 10 µl of M-Elution buffer (provided with the kit) was added to the column matrix. The column was centrifuged (≥10,000 g for 30 seconds) and the eluted DNA was stored at -20°C until further use.

Infinium human methylation Epic BeadChip methylation measurements

The methylation level of 865 918 CpGs in the human genome was assessed on 1 μg of bisulphite-converted gDNA using the Human DNA methylation EPIC BeadChip according to the Infinium Methylation Assay Protocol. The 6 samples were equally loaded in the EPIC Beadchip to prevent the batch effect. Assayed bead chips were scanned to generate the data in the form of .idat files for each biological sample via Illumina’s GenomeStudio software. Following this, the .idat files were further processed using the R (4.1.2)/Minfi Bioconductor package (The R Project for Statistical Computing, Auckland, New Zealand). The β-value between 0 (totally unmethylated) and 1 (totally methylated) were generated after the background subtraction of control probes, followed by normalization with a subset quantile normalization (SWAN) [14]. The quality control procedure eliminated the ensuing probes: i) probes with sex-specific CpG loci, ii) probes with a detection rate at p>0.01, iii) probes containing SNP. After these steps, 715 810 probes from a total of 865 918 CpG probes remained for methylome-wide analysis [15]. The DMPs between untreated control and 6-OHDA treated group were determined using dmpFinder algorithm from the minfi package in R Bioconductor. A pooled t-test, with a p value of ≤0.01 level of significance; bootstrapping, 300 (iterations), smoothing; between untreated and 6-OHDA treated group was applied as a cut-off. CpGs with average Δbeta (delta beta; difference in beta values between 6-OHDA treatment and untreated control group) of >0.1 was regarded as hypermethylation and Δbeta value <-0.1 being hypomethylation [16]. A chromosomal distribution and volcano plot depicting the difference in the methylation between 6-OHDA and the control groups (p<0.01) were generated using R Bioconductor Minfi package.

Bioinformatic analysis

The statistically significant DMRs with p<0.01 and Δbeta values cut-off of 0.1 were imputed into the Metascape tool (http://metascape.org/gp/index.html) for elucidation of gene ontology (GO) using GO biological Processes, KEGG, REactome Gene Sets, Canonical Pathways, PANTHER Pathway, CORUM and WikiPathways. The p values between the GO ontology clusters were calculated using hypergeometric distribution, while q-values were calculated using the Benjamini-Hochberg procedure to account for multiple tests. Then, the protein-protein interaction (PPI) enrichment analysis for the significant genes was performed using STRING, BioGrid, OmnoPath and InWeb_IM databases. The Molecular Complex Detection (MCODE) algorithm has been applied for a network that contains between 3 and 500 proteins. The best scoring terms by p values for each MCODE were identified as the functional description of each component [17].

Statistical analysis

The statistical difference between the study variables (untreated vs. 6-OHDA) was analysed using the R Bioconductor package minfi [18]. The CpG probes that demonstrated Δbeta <-0.1 (hypomethylated) and >0.1 (hypermethylated) with p<0.01 between 6-OHDA-treated and untreated differentiated neuroblastoma cells groups were regarded as statistically significant.

Genome-wide association analyses of 6-OHDA induced cytotoxicity in differentiated neuroblastoma cells

To determine the differences in the DNA methylome between untreated and 6-OHDA treated differentiated neuroblastoma cells, we performed genome-wide DNA methylation using the Infinium MethylationEPIC beadchip method. Among 865 918 CpG probes interrogated, 715 810 CpG probes were obtained after data filtration and normalization using the R Bioconductor/Minfi package. Then, these CpGs were further processed using the minfi package, and we discovered that 236 DMPs or 163 DMRs were significantly methylated by 6-OHDA with Δbeta cut-off 0.1 and p<0.01. The Manhattan plot of -log10 transformed p values of the 715 810 CpG probes at respected chromosomes is shown in Fig. 1A. This plot shows the distribution of DNA methylation changes induced by 6-OHDA in differentiated neuroblastoma cells across all chromosomes, except the sex chromosomes. The grey line indicates the cut-off of statistical significance of p<1×10^-3 and the red line indicates p<1×10^-5.

The volcano plot in Fig. 1B shows the bird’s view visualization of DNA methylation changes in 6-OHDA exposed differentiated neuroblastoma cells. The volcano plot shows the log2 fold change of DMPs against -log10 (p values). Among 236 significant DMPs, 62 CpG probes displayed log2 (fold change) of more than 0.2, and 75 CpG sites with log2 (fold change) of less than -0.2.

We also explored the genomic distribution of the differentially methylated CpGs in relation to UCSC (University of California, Santa Cruz) CpG island (Shelf, shore, island, open seas) in 6-OHDA treated differentiated human neuroblastoma cells compared to the reference group as shown in Fig. 2A. The differentially methylated CpGs were overrepresented in open seas (65%) and under-represented in N_Shelf (1%) and S_Shelf (2%). We also investigated the methylation patterns at the UCSC refgene groups. The ref gene groups explored were the gene body, ExonBand (Exon boundary), TSS200 (within 200 bases of transcription start sites), TSS1500 (within 1500 bases of transcription start site), 5’UTR (5’-untranslated region) and 3’UTR (3’-untranslated region). Our data displayed that the CpG methylation was overrepresented at the gene body (80 CpG sites) and underrepresented at the ExonBand (3) and 3’UTR (5) regions (Fig. 2B).

We identified the top 10 hyper- and hypomethylated CpG sites with p<0.01 in 6-OHDA-treated differentiated neuroblastoma cells, as shown in Table 1 and 2. In the top 10 hypermethylated CpGs, 3 CpGs are from chromosome 3, in which two CpGs are from the gene body and one from 5’UTR. The top 5 hypermethylated CpGs with p<0.01, Δbeta>0.1 are ABHD5 (0.3), ITPR1 (0.29), NRBP1 (0.26), ARID1B (0.26) and FADS3 (0.26) gene regions. Of note, ABHD5 (α/β-Hydrolase Domain Containing 5, Lysophosphatidic Acid Acyltransferase) gene is a key gene that encodes a protein that activates an enzyme called triglyceride lipase (ATGL). The ATGL is an important triacylglycerol (TG) hydrolase that catalyzes the generation of glycerol and free fatty acids. Mutation in this gene causes the pathologic accumulation of TG and leads to detrimental changes in neurons, such as a lack of axon regeneration and poor neuronal outcomes [19]. Hence, the hypermethylation of the ABHD5 gene body in the 6-OHDA-induced differentiated neuroblastoma cells may indicate lipid metabolism defects in neurons. In hypomethylated CpGs, nine out of ten CpGs are hypomethylated at the gene body. The top five hypomethylated CpG are RAD51B (-0.15), PMS1 (-0.15), C11orf49 (-0.13), STRN (-0.12) and SHISA7 (-0.11). Interestingly, the top hypomethylated CpG, RAD51B is a DNA repair gene that maintains genetic stability during the cell cycle [20]. Since DNA damage is one of the biochemical hallmarks in oxidative stress-induced apoptosis, the differential methylation of DNA repair genes may be an indication of differential gene regulation to counteract the cell death in the PD model.

Bioinformatic analysis of differentially methylated regions regulated by 6-OHDA

The functional GO enrichment analysis (GO, KEGG, WikiPathways, Reactome terms) was conducted on the total set of 163 statistically significant DMRs as well as the hyper- and hypomethylated DMRs using Metascape bioinformatic database (Fig. 3). The overall DMRs analysis shows 20 GO-biological terms; whereby the top three overexpressed gene clusters included the regulation of neuron projection development (AKT1, ALK, BDNF, GAK, GRID2, IGF1R, SERPINI1, ROBO1, SDC2, BRSK2, CUX2, SIPA1L1, SEMA5B, DDX56), long-term depression (GNAS, GRID2, IGF1R, ITPR1, PLCB4, MAPK1) and serotonin receptor 4/6/7 and NRC signaling (EGR1, GNAS, NR3C1, MAPK1). The enrichment of GO:0010975: regulation of neuron projection development reflects the significant methylation of genes modulating the rate or extend of neuron projection development regulated by 6-OHDA on differentiated neuroblastoma cells. Loss of dopaminergic innervation that disrupts striatal projection systems’ response to cortical and thalamic signals has been evidenced in PD [21]. The discovery of differentially DNA methylated genes regulating the neuron projection development in our experimental model is crucial as it is congruent with the pathological event in PD (Fig. 3A).

The functional GO enrichment analysis on 71 hypermethylated genes depicted 11 GO-biological terms, with the top three overexpressed gene clusters including the regulation of ARF protein signal transduction (MAP4K4, IQSEC1, IQSEC3), response to fluid shear stress (AKT1, ACE, SOCS5) and antigen processing and presentation of peptide antigen via MHC class 1 (ACE, MICB, SAR1B) (Fig. 3B). Our study showed a significantly increased DNA methylation in genes regulating ARF (ADP-ribosylation factor) signal transduction that functions in recruiting cargo-sorting coat proteins to regulate membrane traffic and organelle structure and modulate membrane lipid composition via guanine nucleotide exchange. This data is in line with the identification of protein aggregates and fragmentation of the Golgi apparatus in the brain of PD and AD that is associated with the disruption in ARF GTPases membrane trafficking functions leading to neurodegeneration [22]. Therefore, we can conclude that the 6-OHDA-induced apoptosis in our study is mediated by changes in crucial regulator ARF GTPase proteins expression that disrupts constant membrane trafficking pathways in neurons and eliminates dysfunctional proteins [23].

We also conducted a functional GO enrichment analysis on 92 hypomethylated genes induced by 6-OHDA on differentiated neuroblastoma cells. Our analysis revealed 17 GO-biological terms, where the top three enriched gene clusters are regulation of protein transport (PIK3C3, SREBF1, SUMO1, GPR68, NROB3, HYAL2, SERGEF, RBM22), mitotic spindle organization (EML1, KIF2A, NEK6, CLASP1) and regulation of cell division (CXCR5, FGF1, PIK3C3, TP63, KIF13A) (Fig. 3B). The DNA methylation alteration of genes regulating protein transport within cells or between cells may cause disruption in protein translation, chaperone-assisted protein folding and protein degradation pathways [24].

We performed the PPI analysis on the 6-OHDA induced DNA methylation changes on 163 DMRs and obtained a rather messy hairball network (Fig. 4A). To further interpret the PPI, we focused on the MCODE algorithm with significant p values that was extracted by the database according to functional enrichment analysis. The input of 163 DMRs in the Metascape interface generated four MCODE modules that comprising 12 genes, from which 2 seed nodes were identified; ITPR1 and TDRD3. The seed proteins are proteins that are responsible in forming the clusters. The top three pathway clusters overexpressed in the MCODE 1 are three KEGG pathways that include Cholinergic synapse (hsa04725), dopaminergic synapse (hsa04728) and Apelin signaling pathway (hsa04371). The Cholinergic synapse pathway functions in formation, distribution and uptake of acetylcholine (ACh) neurotransmitter in neuronal synapse. In the human brain, ACh facilitates memory, learning and motor control [25]. Dopamine is generated in dopaminergic synapses and interacts with dopamine receptors in the central nervous system (CNS) to control locomotor, motivation, reward, sleep and endocrine regulation [26]. The Apelin signalling pathway is a neuronal survival response mechanism mediated by an endogenous native peptide ligand, Apelin [27]. The three genes that were differentially methylated in these three pathways are Inositol 1,4,5-triphosphate receptor type 1 (ITPR1), Guanine Nucleotide Binding Protein Subunit Gamma 7 (GNG7) and AKT Serine/Threonine Kinase 1 (AKT1) (Fig. 4B). The 3 GO terms/pathways displayed p values of less than 0.05, with cholinergic synapse expressing log10 p value=-7.3, dopaminergic synapse showing log10 p value=-7.1) and Apelin signaling pathway with log10 p value=-7.1 as curated by Metascape customized analysis (Fig. 4C). The increased DNA methylation levels on these genes in differentiated neuroblastoma cells exposed to 6-OHDA indicate the changes in gene regulation in these 3 pathways. Multitude studies have shown that progressive loss of synapses within striatum leads to neuronal loss in neurodegenerative diseases, such as PD [28]. Therefore, our findings have further confirmed that the loss of synapses in a neurodegeneration event may have association with changes in DNA methylations in genes involved in maintaining the number of axons and synapses for optimal function of neurons in central nervous system (CNS).

The other two pathways overexpressed by 6-OHDA are chromatin remodeling (GO:0006338) and chromatin organization (GO:0006325) in MCODE 2 module and protein phosphorylation (GO:0006468) processes in MCODE 3 modules.

The comparison of beta values of the ITPR1, AKT1 and GNG7 genes between 6-OHDA and neg control are illustrated in Fig. 5. All DMRs displayed an increased DNA methylation in 6-OHDA treated differentiated neuroblastoma cells compared to the negative control group. The Δbeta ITPR1 is 0.29 (p=0.002), Δbeta AKT1 is 0.11 (p=0.0000064) and Δbeta GNG7 is 0.1 (p=0.000198).

We further verified the overexpressed pathways by the 163 significantly methylated DMRs in our 6-OHDA induced cytotoxicity on differentiated neuroblastoma cells using KEGG online database. Upon uploading the DMRs in the KEGG mapper interface, we found that the metabolic pathway was remarkably enriched with significant methylation of 12 DMRs, followed by neurodegeneration-multiple diseases (5 DMRs) and Ras signaling pathway (5 DMRs). Interestingly, we found that most of the DMRs were found in more than one pathway, such as GNG7 which was identified in the Ras signaling pathway, apelin signaling pathway and serotonergic synapse; AKT1 was enriched in Ras signaling pathway, Apelin signaling pathway, autophagy and cAMP signaling pathway; ITPR1 gene in pathways of neurodegeneration, apelin signaling pathway, autophagy, serotonergic synapse and calcium signaling pathway (Table 3). This data signifies the multiple roles these genes play in the biological pathways activated in the 6-OHDA-induced damage in differentiated neuroblastoma cells.

Epigenetic modification, particularly DNA methylation, is an indispensable molecular mechanism that occurs in response to genetics, environment, aging, diet and environmental toxins. Recently, studies have been directed to understanding the DNA methylation on PD-associated genes that alters the gene expression leading to neurodegeneration. An emerging study has confirmed that hypomethylation of cytochrome P450 2E1 (CYP2E1) gene was associated with increased expression of CYP2E1 protein in PD brains but not in blood cells and skin fibroblasts [29]. Considering the significance of DNA methylation in PD, we embarked on a study to investigate the DNA methylation alterations that arise in a cellular model of PD. However, the short-term treatment (24 h) of 6-OHDA on differentiated human neuroblastoma cells was expected to yield a mild to moderate level of methylation changes. Therefore, we report all CpG probes with Δbeta cut-off of 0.1 with p<0.01. Our finding showed that the 6-OHDA induced neurotoxicity model yielded significant DNA methylation changes on 236 CpG probes or 163 DMRs, p<0.01. Among these CpG probes, 110 CpGs were increased, while 126 decreased in DNA methylation level in 6-OHDA treated cells. The DNA methylation was primarily discovered in open sea region and gene body according to UCSC classification of CGI regions (Fig. 2B) and refgene groups (Fig. 2C), respectively.

The 6-OHDA implicated DNA methylation changes on these 236 CpG probes was further scrutinized using Metascape, a web-based portal designed to provide functional gene annotation and unique pathways [30]. Our analysis revealed three significant DMRs (ITPR1, AKT1, GNG7) that are involved in three neuronal pathways as curated in KEGG pathway database, namely Cholinergic synapse (hsa04725), Dopaminergic synapse (hsa04728) and Apelin signaling pathway (hsa04371). These three genes are significantly hypermethylated by 6-OHDA than in negative control. The changes in DNA methylation levels in genes regulating the growth and maintenance of neuron synapses may alter the expression of these gene and eventually affects neuronal functions. Moreover, our study demonstrated that the neurodegeneration process in 6-OHDA-induced PD model is mediate by DNA methylation of critical genes that plays vital roles in cholinergic and dopaminergic synapse formation and functions.

ITPR1 is predominantly detected in the CNS and plays an important role in regulating critical neuronal functions in synaptic plasticity and dendritic outgrowth [31]. Inositol 1,4,5-triphosphate receptors consist of a family of Ca²⁺ channels that regulate intracellular Ca²⁺ levels for physiological events to occur. We revealed that the 6-OHDA induced toxicity on differentiated neuroblastoma cells implicated significant DNA methylation changes (Δbeta=0.29, p=0.0019) on ITPR1 gene body-open seas. A previous study reported that misregulated DNA methylation cycle on ITPR1, caused downregulation of ITPR1 expression that resulted in cerebellar ataxia in murine model [32]. Furthermore, a compelling study has revealed that ITPR1-knock-out (KO) mice presented with severe cerebellar ataxia, opisthotonos, repetitive rigid posture and tonic contractions of the neck and trunk [33]. The role of ITPR1 in experimental PD has been well clarified by Fedorenko et al. According to their study, the 6-OHDA-induced state of experimental hemiparkinsonism in rats damaged the dopaminergic neurons of the substantia nigra and resulted in higher expression of ITPR1 gene. The overexpression of ITPR1 gene was associated with abnormal release of Ca²⁺ from intracellular calcium storage during the neurodegenerative process [34]. Aberrant Ca²⁺ levels in neurons affected by dysregulated ITPR1 has been identified as one of the possible mechanisms linking polyglutamine expansion of Huntingtin (Htt) protein and degeneration of GABAergic medium spiny striatal neurons in Huntington’s disease [35].

The AKT1, a serine-threonine protein kinase, has been known for its role in cellular growth (division, differentiation), survival and apoptosis in neuronal cells. We showed that the 6-OHDA caused a significant increase in DNA methylation (Δbeta=0.11, p=0.0000064) at TSS1500-S shore compared to the control group. The effect of DNA shore methylation was reported in cancer research, whereby increased shore methylation alters tumor suppressor genes in breast tumors more than in normal tissue [36]. The effect of DNA methylation in AKT1 on gene expression was reported in a methamphetamine-induced psychosis study in human subjects. In this study, the DNA methylation on promoter regions of AKT1 was associated with increased expression of this gene in psychosis patients [37]. Nevertheless, compelling data from previous studies have suggested that AKT1 is involved in various signaling pathways that protect against PD [38]. Moreover, the previous study has shown that AKT activation is reduced in the PD brain, and its expression is mainly regulated by the PD-causing gene, PINK1. The lack of PINK1 activity causes the suppression of phospholipid PI (3,4,5)P₃ levels that give rise to fragmentation of Golgi apparatus, the pathological changes observed in PD. PINK1 gene protects Golgi apparatus from fragmentation through activation of AKT via kinase-dependent activity regulations of the primary activator, phospholipid PI(3,4,5)P₃ [39]. Besides that, Grider et al. have reported that increased activity of AKT in primary sensory neurons promotes dendritic branching and elongation via activation of localized Rac1 and mTOR signaling cascades [40]. Disruption in AKT activation or activity results in apoptosis in cultured neurons. At the same time, transfection of an active form of AKT promotes survival signaling pathway, neuron size and innervation of the striatum in postnatal development of dopamine neurons in mice [41].

The GNG7 encodes for γ7 protein that binds with the G-protein to form a heterotrimer (Gα_olfβ₂γ₇) complex. The motor activity at the striatum is regulated via hierarchical assembly of the Gα_olfβ₂γ₇ with dopamine receptor D₁ (DRD1) or adenosine A_2A receptor (AA2AR) that induce the activation of adenyl cyclase [42]. In our study, we showed a significant increase in DNA methylation at the 5’UTR-open seas region of GNG7 (Δbeta=0.1, p=0.000198). Although, no previous study reported on modulation of the GNG7 in neurological disorder, our bioinformatic analysis using Metascape tool has revealed a significant PPI in MCODE 1 interaction (Fig. 4). Nevertheless, a previous investigation by Schwindinger et al. on mice suggested that knock-down of GNG7 and GNG3 genes presented with a progressive seizure that declined the mice’s life span to only 75 days. These effects were absent in a single knock-out of GNG7 or GNG3 gene alone which showed normal mice survival. This report concluded that the disparity in gene knock-down effects was attributed to the different functions elicited by G-protein constituted γ3 and γ7 subunits [43]. A compelling report from recent studies suggested that the selective loss of GNG7 (γ₇) gene expression disrupts the heterotrimer assembly of Gα_olfβ₂γ₇, which affect the coupling with DRD1 or AA2AR receptors that interferes with the cAMP signaling and eventually compromises the downstream motor behaviors. It is important to note that dysfunction in the cAMP signaling in the human striatum has been associated with neurological or neurodegenerative diseases such as PD, dystonia, schizophrenia and drug addiction [44].

Although we are considered the pioneers in reporting the GWAS on DNA methylation in 6-OHDA induced PD cell-based model, this study was not without its limitations. Since this study was aimed to understand the different patterns of DNA methylation in PD-associated DMRs in the cellular model of PD, hence our data were not validated using other method. Secondly, the treatment duration with 6-OHDA on differentiated human neuroblastoma cells was only 24 hours. This protocol was well established for PD related studies [1]. As epigenetic changes require a sufficient amount of time to occur, we found the CpG methylation levels in our study were in the range of mild to moderate. Therefore, we set the Δbeta threshold as 0.1, with p<0.01 significance level. This finding can serve as preliminary data for studies investigating the epigenetic changes implicated by therapeutic targets using a cell-based model of PD.

In conclusion, our study provides the first evidence for DNA methylation changes in the 6-OHDA induced cellular model of PD. Furthermore, our findings suggest that there was almost equal level of hypermethylation and hypomethylation in the identified DMPs. We also demonstrated that this experimental model significantly hypermethylated PD-related CpGs: AKT1, ITPR1 and GNG7.

This research was supported by JCSMHS internal grant scheme.

Top 10 significantly hypermethylated DMRs (Δbeta>0.1) in 6-OHDA treated differentiated neuroblastoma cells with p<0.01

Top 10 significantly hypomethylated DMRs (Δbeta>-0.1) in 6-OHDA treated differentiated neuroblastoma cells with p<0.01

Top overexpressed biological pathways obtained from KEGG Mapper