Dana M. Freeman
Zhibin Wang- Laboratory of Environmental Epigenomes, Department of Environmental Health and Engineering, Bloomberg School of Public Health, Johns Hopkins University, Baltimore, MD, United States
CCCTC-binding factor (CTCF) is a regulatory protein that binds DNA to control spatial organization and transcription. The sequence-specific binding of CTCF is variable and is impacted by nearby epigenetic patterns. It has been demonstrated that non-coding genetic variants cluster with CTCF sites in topological associating domains and thus can affect CTCF activity on gene expression. Therefore, environmental factors that alter epigenetic patterns at CTCF binding sites may dictate the interaction of non-coding genetic variants with regulatory proteins. To test this mechanism, we treated human cell line HEK293 with rotenone for 24 h and characterized its effect on global epigenetic patterns specifically at regulatory regions of Parkinson’s disease (PD) risk loci. We used RNA sequencing to examine changes in global transcription and identified over 2000 differentially expressed genes (DEGs, >1.5-fold change, FDR < 0.05). Among these DEGs, 13 were identified as PD-associated genes according to Genome-wide association studies meta-data. We focused on eight genes that have non-coding risk variants and a prominent CTCF binding site. We analyzed methylation of a total of 165 CGs surrounding CTCF binding sites and detected differential methylation (|>1%|, q < 0.05) in 45 CGs at 7 PD-associated genes. Of these 45 CGs, 47% were hypomethylated and 53% were hypermethylated. Interestingly, 5 out of the 7 genes had correlated gene upregulation with CG hypermethylation at CTCF and gene downregulation with CG hypomethylation at CTCF. We also investigated active H3K27ac surrounding the same CTCF binding sites within these seven genes. We observed a significant increase in H3K27ac in four genes (FDR < 0.05). Three genes (PARK2, GPRIN3, FER) showed increased CTCF binding in response to rotenone. Our data indicate that rotenone alters regulatory regions of PD-associated genes through changes in epigenetic patterns, and these changes impact high-order chromatin organization to increase the influence of non-coding variants on genome integrity and cellular survival.
Introduction
Parkinson’s disease (PD) is the second most common neurodegenerative disorder in the United States (de Lau and Breteler, 2006). More than 800 genetic association studies have been conducted to interpret genetic contribution to PD etiology (Lill et al., 2012; Coetzee et al., 2016). Genome-wide association studies (GWAS) evaluate the association of common genetic variants to a phenotype or disease outcome. Since 2005, thousands of variants have been identified to have a significant association with a disease and more than 1600 single-nucleotide polymorphisms (SNPs) have been identified as genetic risk variants for PD (Lill et al., 2012). However, unlike rare monogenetic associations, the functional consequences of most of these variants have yet to be determined. Over 90% of all indexed SNPs including those associated with PD occur in non-coding regions of the genome (Maurano et al., 2012; Verstraeten et al., 2015). This discovery led to the hypothesis that SNPs in the human genome interact with regulatory elements to control gene expression (Wang et al., 2019). This is supported by expression of quantitative trait loci (eQTLs) defined as genetic regions that are enriched at positive GWAS sites and explain variability in the expressivity of a gene (Nica et al., 2010). Despite these advances, it remains a challenge to determine which genetic variants in a broad region of variants is the driver of gene expression changes particularly when regulatory element interactions are long-range (Do et al., 2016). SNPs cluster within enhancers and can modify PD risk. These observations of PD-associated SNPs have been described in multiple cell types (Coetzee et al., 2016). With this new evidence, studies are now focusing on interactions of regulatory elements to understand how genetic associations trigger disease biology within the brain.
CCCTC-binding factor can play a long-range cis-regulatory role that insulates genes from their surrounding signaling environment by directing chromatin looping (Phillips and Corces, 2009). Functional CTCF binding sites are required for the formation of distinct structural domains within a three-dimensional chromosomal organization (Ong and Corces, 2014; Tang et al., 2015). CTCF binding is dependent upon DNA sequence (CCGCGNGGNGGCAG) and allelic hypomethylation (Wang et al., 2012). Thus, genetic variants and epigenetic patterns within binding sites can contribute to dysfunctional CTCF allele-specific binding (Tang et al., 2015; Wang et al., 2019).
Approximately 85% of PD cases cannot be explained by genetic predisposition alone (Franco et al., 2010; Verstraeten et al., 2015; Labbé et al., 2016). Therefore, it is likely that most cases are caused by the interplay of common SNPs with environmental factors. Environmental factors can modulate the association of a genetic variant with a disease (Lee et al., 2011). For instance, exposures that impact allele-specific methylated regions in the genome can influence CTCF binding and thus influence non-coding variants’ effect on genetic expression (Wang et al., 2019). GWAS association signals are complex in that they can cover a broad region of DNA with several polymorphisms, so we focused on environmentally induced epigenetic changes in CTCF binding regions nearby risk-associated genes to explore mechanisms of gene–environment interactions in PD.
In our pesticide-induced cellular model, we used rotenone, a naturally occurring insecticide and potent inhibitor of complex I in the mitochondrial electron transport chain. The primary use of rotenone today is as a piscicide to terminate invasive or noxious species of fish. Permissible application concentrations up to 250 ppb can be applied to public and recreational waters (United States Environmental Protection Agency [USEPA], 2007). Rotenone is a widely accepted PD toxicant and can robustly replicate pathology via depletion of ATP, generation of reactive oxygen species, damage of nigrostriatal tissues, and death of dopamine producing cells in the midbrain (Dawson et al., 2002; Cicchetti et al., 2009). It has also been shown to cause these types of cellular pathology in HEK293 (Orth et al., 2003; Teixeira et al., 2018). While cell line HEK293 is a human immortalized cell line derived originally from primary embryonic kidney cells, it has been found to have a genetic signature similar to neurons (Stepanenko and Dmitrenko, 2015). We chose this cell line given their well-characterized genome and ENCODE regulatory elements (ENCODE Project Consortium, 2012; Lin et al., 2014).
DNA methylation and histone acetylation are epigenetic modifications implicated in rotenone-induced neurotoxicity (Huang et al., 2019). DNA hypomethylation has been reported in response to pesticide exposure (Hou et al., 2012), and we discovered that rotenone reduces DNA methylation at DNMT1-dependent regions in the human genome (Freeman et al., 2020). Histone acetylation patterns have been more extensively studied in rotenone-induced PD due to its high correlation with gene expression and enhancer activation (Wang et al., 2008). Most studies agree that rotenone-induced neurodegeneration is associated with pathological hyperacetylation as a result of impaired homeostatic activity of HATs and HDACs (Feng et al., 2015; Park et al., 2016; Harrison et al., 2018; Wang et al., 2018; Huang et al., 2019).
In this study, we examined rotenone-induced changes in DNA methylation and histone acetylation patterns at CTCF binding sites adjacent to PD-associated genes. Eight selected genes had identified disease-risk SNPs in a non-coding region and were indexed by a meta-data analysis of over seven million human polymorphisms (Lill et al., 2012). We hypothesize that rotenone exposure modifies epigenetic patterns at CTCF binding motifs and affects its allele-specific transcription factor binding. We postulate that this mechanism could mediate the interchange between genetic variants and regulatory elements controlling transcription and genomic stability.
Materials and Methods
Cell Culture and Treatment of Human Cell Line HEK293
All media reagents and chemicals in cell culture were purchased from Sigma (St. Louis, MO, United States). Human cell line HEK293 was grown in Dulbecco’s Modified Eagle Medium with high glucose, L-glutamine, and sodium pyruvate. Media were supplemented with 10% (v/v) heat-inactivated fetal bovine serum and 1% (v/v) penicillin–streptomycin. HEK293 cells were confirmed by ATCC. Cells were treated at approximately 70% confluency with 200 nM rotenone or DMSO vehicle control (<0.001%) for 24 h. Cell viability was measured with trypan blue (0.4%) staining, and cells were counted manually with a hemocytometer. Cells used for experiments were at least 85% viable relative to the vehicle control.
Qualitative Analysis of Total 5mC Levels
Genomic DNA was extracted from two replicates of DMSO or rotenone-treated HEK293 using 1:1:1 phenol: chloroform: isoamyl alcohol (Sigma, St. Louis, MO, United States). Global DNA methylation was first measured by dot blot analysis. Bisulfite-treated DNA (30–60 ng/μL) was denatured at 95°C for 5 min and then cooled at 4°C for 5 min in a conventional thermocycler (MyCycler; Bio-Rad; Hercules, CA, United States). DNA was spotted onto 0.45-micron nitrocellulose paper as 1- or 2-μL drops and dried for 30 min at room temperature. The membrane was UV cross-linked at 3000 Hz and incubated in anti-5methylcytosine (5mC) primary antibody overnight at 4°C (Epigentek 33D3; Farmingdale, NY, United States). The membrane was washed with TBST and incubated with a secondary antibody conjugated to HRP for 1 h at room temperature (Santa-Cruz Biotechnology anti-mouse IgG sc-2005; Dallas, TX, United States). The membrane was washed again with TBST after secondary incubation and visualized with chemiluminescence (ProSignal Femto; Prometheus; Raleigh, NC, United States). We used a serial dilution of 100% 5mC standard (Zymo Research D5012; Irvine, CA, United States) to create a standard ladder (Supplementary Figure 1).
HEK293 Western Blot for Global H3K27ac
HEK293 cells were collected after a 24-h treatment, and histones were extracted using the Abcam Histone Extraction kit according to kit instructions (Cambridge, United Kingdom). Histone protein concentration was measured by Qubit Protein Assay from Thermo Fisher Scientific (Waltham, MA, United States). Protein (5 μg) was loaded onto 4–15% Bio-Rad Page Gels and transferred to 0.45 μm nitrocellulose (Bio-Rad, Hercules, CA, United States). The blots were incubated with H3K27ac primary antibody (1:1000; Abcam ab4729) overnight at 4°C and anti-rabbit IgG conjugated secondary antibody (1:5000, Santa Cruz Biotechnology Sc-2357) for 1 h at room temperature. The histone protein was normalized to total histone 3 (H3; Abcam ab1791) and quantified with ImageJ software. We tested significance by comparing the ratio of H3K27ac/H3 with a two-tailed Student’s paired t-test (p < 0.05).
RNA Extraction and RNA Sequencing Library Construction
Total RNA was extracted from two replicates of DMSO or rotenone-treated HEK293 using the TRIzol method (Invitrogen, Carlsbad, CA, United States). A total of 2 μg per sample was used for library construction using the TruSeq Sample Preparation kit from Illumina (San Diego, CA, United States). Poly-A-containing mRNA molecules were isolated from total RNA using oligo-dT attached magnetic beads. Isolated mRNA was then fragmented and synthesized into double-stranded cDNA according to kit instructions. Ligation of unique Illumina adapter indices was completed for each sample before bead purification. Libraries were loaded onto a 2% agarose gel, and library products between 200–800 bp were purified using the mini-Elute gel extraction kit from Qiagen (Hilden, Germany). Approximately 150 ng was sent for sequencing on a HiSeq 2000 platform with 100-bp paired-end reads.
RNA Sequencing Data Analysis
Adapter sequences were removed from the raw sequencing data, and individual libraries were converted to the fastq format. Sequencing reads were aligned to the human genome (hg19) with TopHat2 (v2.0.9) (Kim et al., 2013). For mRNA analyses, the RefSeq database (Build 37.3) was chosen as the annotation references. Read counts of annotated genes were obtained by the Python software HTSeq count (Anders et al., 2015). DEGs were defined as those with a 1.5-fold change in expression using FDR < 0.05 from the edgeR package (Robinson et al., 2010). Gene Ontology annotation was done with Gorilla online platform and visualized with Revigo and Cytoscape (Eden et al., 2009; Supek et al., 2011; Otasek et al., 2019).
RNA Sequencing Validation With Quantitative Reverse Transcription-PCR
Total RNA was extracted from an additional replicate of HEK293 treated with DMSO or rotenone using the same procedure as stated above. A total of 500 ng RNA was converted to cDNA with the PrimeScript RT reagent kit with gDNA eraser from Takara (Kusatsu, Japan). We selected 10 genes for quantitative PCR (qPCR) analysis using primers listed in Supplementary Table 1. All qPCR reactions were performed on a 7500 Real-Time PCR system from Applied Biosystems (Foster City, CA, United States) using the iTaq Universal SYBR Green Supermix from Bio-Rad (Hercules, CA, United States). The change in expression was normalized to the GAPDH housekeeping gene and expressed as fold change (2–ΔΔCT).
Identification and Selection of PD-Associated Genes
We identified PD-associated genes using the National Health Genomic Research Institute GWAS Catalog (Buniello et al., 2018). We searched for all associations both reported and mapped using the trait “Parkinson’s disease” (EFO_0002508) which included 39 publications investigating genomic signatures of both familial and environmentally driven PD as well as Lewy body pathology and Parkinsonism in frontotemporal lobe dementia (Supplementary File 2). We calculated the frequency for various region types (non-coding, regulatory, coding) within the 246 known genetic variants provided by GWAS Catalog. We compared 399 reported and mapped genes to our list of DEGs. We then cross-referenced these genes with the PD gene online resource which analyzed over 800 publications and seven million polymorphisms (Lill et al., 2012). We selected five genes that remained significant in the PD gene meta-analysis, were represented in at least two studies, and had their most significant variant in a non-coding region (Table 2). We also selected three additional genes from the PD gene database that were represented in our RNA sequencing data (Table 2). The first, UBOX5, was among the most significant polymorphisms identified by the meta-analysis (Lill et al., 2012; Nalls et al., 2014). The other two, PARK2 and CHCHD2, have significant polymorphisms according to the PD gene database but are also reported to have autosomal mutations that contribute to familial disease cases (Lill, 2016).
Region Selection for Bisulfite and ChIP Primer Design
CCCTC-binding factor transcription factor binding was observed using the Uniform Transcription Factor Binding data found in the ENCODE Regulation super track in UCSC Genome Browser. We selected all CTCF transcription factor binding sites detected with ChIP-seq experiments from the ENCODE consortium from 2007 to 2012 (ENCODE Project Consortium, 2012). We also predicted which cytosine would overlap the binding motif using the CTCF binding prediction tool database v2.0 (Ziebarth et al., 2012). Primer design was focused on CTCF binding sites for both bisulfite sequencing and ChIP-qPCR experiments (further described below).
Bisulfite-DNA Conversion and Bisulfite-Amplicon Sequencing Library Construction
Genomic DNA was extracted from two replicates of DMSO or rotenone-treated HEK293 using phenol: chloroform: isoamyl alcohol (Sigma, St. Louis, MO, United States). A total of 200 ng DNA was bisulfite-converted using the Sigma DNA Imprint Modification kit two-step protocol. Bisulfite-converted DNA (BS-DNA) was amplified with primers for selected regions designed with MethPrimer (Li and Dahiya, 2002) (Supplementary Table 2). Amplified BS-DNA products were run on a 2% EtBr agarose gel and purified using the mini-Elute gel extraction kit from Qiagen (Hilden, Germany). Purified products for each sample were pooled together, and 1 ng was used for library preparation using the Illumina Nextera XT DNA Library Preparation kit. Each sample was tagged with a unique Nextera XT adapter (San Diego, CA, United States). Sequencing libraries were quality checked via Bioanalyzer and run on an Illumina MiSeq platform to generate 150-bp paired-end reads.
Bisulfite-Amplicon Sequencing Analysis for Methylation Patterns at CTCF Binding Sites
The raw fastq files were imported into the Galaxy web platform (Afgan et al., 2016). Reads with quality score <30 were filtered out, and reads with quality score >30 were trimmed with Trim Galore (Krueger, 2015). Reads were mapped to the human genome (hg19) using bwa-meth (Pedersen et al., 2014). MethylDackel was used for methylation calling, and per-cytosine contexts were merged into per-CPG metrics1. Duplicates and singletons identified in alignment were ignored from the methylation call. Minimum and maximum per-base depths were 1000× and 100,000×, respectively. The output was selected for methylKit format. Coverage statistics and differentially methylated regions were calculated for CG sites with methylKit installed in R (v3.5) (Akalin et al., 2012). Differentially methylated cytosines were defined as being present in both biological replicates, having a minimum absolute difference of 1% using the coverage weighted mean and having a SLIM adjusted q-value < 0.01 using the methylKit logistic regression model (Ning et al., 2011). The change in mean percent methylation (Δme) for all CpG sites within a defined region was calculated by taking the mean number of methylated versus non-methylated CpG sites from the pooled control and treated samples and using Fisher’s exact test FDR < 0.05.
Chromatin Immunoprecipitation
All chemicals were purchased from Sigma unless otherwise noted (St. Louis, MO). HEK293 cells were harvested after a 24-h treatment and resuspended in fresh media at 10 × 106 cells/mL in a conical tube. Cells were fixed with 1% formaldehyde for 10 min at room temperature. Reaction was stopped with 0.2 M glycine and incubation at room temperature for 5 min. Fixed cells were centrifuged for 5 min at 300 × g and 4°C and washed with 1 mL cold PBS. Fixed cell pellet was stored at -80°C until chromatin immunoprecipitation (ChIP).
Cell pellets were resuspended at approximately 1 × 106 cells/0.1 mL with PBS + 0.5% Triton-X + 1% protease inhibitor cocktail and incubated on ice for 10 min prior to centrifugation for 5 min at 400 × g 4°C. The pellet was resuspended in TE buffer pH 8.0 with protease inhibitor and PMSF. Cells were sonicated at high intensity for 30 s on/60 s off until DNA fragments were within 200–800 bp as checked by 2% agarose gel. After sonication, samples were centrifuged for 15 min at 14,000 × g 4°C to pellet insoluble material. Sheared chromatin was transferred to RIPA buffer, and 10% of total chromatin was saved for input DNA extraction.
Chromatin immunoprecipitation was done with Dynabeads Protein A (Invitrogen, Carlsbad, CA, United States) and 4 μg of primary ChIP-grade antibody (H3K27ac Abcam ab4729; CTCF Millipore 07-729; Rabbit IgG Santa Cruz Biotechnology sc-2025). Beads were washed with lithium chloride (LiCl 0.25 M) buffer, and immunoprecipitated DNA was extracted from beads using the phenol: chloroform method. DNA was quantified using Qubit dsDNA high-sensitivity assay (Thermo Fisher Scientific, Waltham, MA, United States).
ChIP-qPCR Analysis
We selected eight genes for quantitative real-time PCR (qPCR) analysis using primers listed in Supplementary Table 3. Primers were designed with NCBI Primer Blast at H3K27ac peaks surrounding the predicted CTCF binding site (Ye et al., 2012). All qPCR reactions were performed on a 7500 Real-Time PCR system from Applied Biosystems (Foster City, CA, United States) using the iTaq Universal SYBR Green Supermix from Bio-Rad (Hercules, CA, United States). H3K27ac and CTCF enrichment was calculated from the Ct threshold value as a percent of the total input DNA. Rabbit IgG samples were used as a negative control (Figures 5, 6B).
Results
Rotenone-Induced Stress Alters Transcription Factor Intracellular Signaling
To characterize the changes in gene expression upon rotenone exposure, we treated HEK293 cells with rotenone 200 nM for 24 h. We then used RNA-seq analyses to identify over 2000 DEGs in response to rotenone (Supplementary File 3). To gain insights of impacted biological processes, we performed gene ontology enrichment analysis on 1853 of these DEGs with known HUGO (hgnc) symbol and cell description using a p-value threshold of p < 10–3. Among enriched biological processes, we observed a significant induction of the oxidative stress response, transcription factor activity, and chromatin organization (<10–3) (Figure 1C). We observed a significant enrichment of genes involved in nucleic acid binding (<10–7) and DNA binding (<10–5) with the nuclear cell component being most represented (<10–5) (Figures 1A,B and Supplementary Tables 5, 6). We analyzed pathway nine enrichment of the top 200 genes with the largest change in expression using reactome pathways (Fabregat et al., 2018; Supplementary Table 4). Three of the top five pathways enriched in our data were major transcription factor pathways including SMAD, NOTCH, and TP53 which have implications in PD reviewed in the discussion section.
Figure 1. Gene Ontology enrichment analysis of RNA-sequencing data. (A) Gene Ontology molecular functions enriched in our differentially expressed genes from rotenone-treated HEK293 cells. (B) Gene Ontology cellular component enriched. Color gradient indicates enrichment p-value: white, >10– 3; yellow, 10– 3–10– 5; light orange, 10– 5–10– 7; dark orange, <10– 7. Enrichment p-value and adjusted FDR for each term are listed in Supplementary Tables 5, 6. (C) Network analysis of significantly enriched biological processes. Blue boxes enrichment p-value < 10– 3 and orange boxes enrichment p-value < 10– 5. Rotenone treatment changes the expression of genes involved in transcription factor signaling in the nucleus.
Alteration of PD-Associated Genes Stands Out Upon Rotenone Exposure
The GWAS Catalog is a public database of approximately 72,000 variant–trait associations from over 3500 publications (Buniello et al., 2018). Out of 246 PD-associated variants with genetic sequence context information in the GWAS Catalog, 220 variants (89%) were in non-coding regions (intron, intergenic, regulatory, and exon) (Supplementary File 2). Intronic variants constituted most of the known polymorphisms. We searched our DEGs for PD-associated genes and identified 14 genes from the GWAS Catalog (Supplementary File 2). Of these genes, 13 were also considered significant PD-associated genes according to meta-analysis data in the PD gene (Lill et al., 2012). We validated the RNA sequencing results for 10 of these genes and were able to validate 8 of them with qPCR analysis (R2 = 0.96) (Figure 2). We selected five genes (ITGA8, GPRIN3, FER, CNKSR3, BMP4) and three additional genes from the PD gene meta-analysis (UBOX5, PARK2, CHCHD2) for further examination of epigenetic patterns (Table 1).
Figure 2. RNA-sequencing validation of PD-associated genes. Fold-change comparison of RNA sequencing results versus qPCR results and the linear calibration curve of RNA sequencing results with qPCR results expressed as fold change in expression.
Table 1. Selected non-coding genetic variants.
Table 2. Regulomedb results for SNPs within selected regions.
Selected Genes Contain Prominent CTCF Binding Sites in Their Regulatory Non-coding Regions
To examine any potential CTCF motifs within these selected genes, we visualized CTCF binding using experimental data from ENCODE and the CTCF binding prediction tool from the Cui Lab at the University of Tennessee (ENCODE Project Consortium, 2012; Ziebarth et al., 2012). Intriguingly, all selected genes had at least one prominent CTCF binding site in a regulatory non-coding region (Table 1). We designed both bisulfite primers and ChIP primers at these sites using ENCODE regulation data from the Broad Institute (ENCODE Project Consortium, 2012). The polymorphisms in these regions that were recognized by the SNP database (dbSNP2) were also analyzed with Regulomedb, a database that annotates SNPs with known or predicted interactions with regulatory elements in intergenic regions (Boyle et al., 2012). We determined that SNPs within selected regions at two genes, GPRIN3 and FER, were present at CTCF binding sites and were active in the brain (Table 2). The rank of an SNP represents the number of available datasets for that polymorphism, and the score is generated based on the integrated results from available datasets. In this analysis, the polymorphism listed at each gene was present in datasets from experimental transcription factor binding, matched transcription factor position-weight matrix (PWM), and DNase footprinting. We checked the HEK293 genome using the online database3 and did not find either variant in our cells (Lin et al., 2014). This information provides additional evidence that CTCF binding sites among common non-coding variants may be critical in disease pathogenesis.
Rotenone Modifies Epigenetic Patterns Across the Genome
We used the dot blot method to qualitatively assess changes in DNA methylation levels in cells exposed to rotenone. The total methylation level was detected with anti-5mC antibody and visualized using chemiluminescence. After 24 h, 5mC levels were strikingly reduced in genomic DNA (Figure 3A).
Figure 3. Global epigenetic patterns in rotenone treated HEK293 cells. (A) DNA methylation was visualized by dot blot method using the anti-5mC antibody. Three biological replicates are shown in this image. The standard ladder using 100% 5mC standard is shown in Supplementary Figure 1. (B) Global histone H3K27ac levels were measured from total extracted histones using Western blot. Total histone H3 was used as the loading control. Three biological replicates are shown in this image. This western blot was quantified using ImageJ software and is shown as the fold change in the amount of H3K27ac relative to the vehicle (DMSO) control. ∗p < 0.05 using a paired Student’s t-test with the ratio of H3K27ac/H3.
Next, we asked the extent to which rotenone exposure impacted histone acetylation. To investigate histone acetylation at active enhancers with PD-relevant SNPs, we selected H3K27ac mark (Wang et al., 2008) and examined H3K27ac levels from extracted histones using western blot. We measured a significant 1.3-fold increase in H3K27ac in rotenone-treated cells compared to the DMSO vehicle control (p < 0.05) (Figure 3B).
Rotenone Alters DNA Methylation Patterns at CTCF Binding Sites in Regulatory Regions of PD-Associated Genes
Because CTCF binding is methylation sensitive and changes in CG methylation correlate with disease risks (Wang et al., 2012), we next examined a total of 284 CG nucleotides from eight regions surrounding our selected genes. Our amplicon-sequencing results demonstrated that 233 of these nucleotides met our minimum requirement of 1000 × coverage (Supplementary Figure 2). We focused our analysis on 165 CG sites that met minimum coverage requirements and overlapped predicted CTCF binding motifs at seven of the selected genes. There were 45 differentially methylated CG sites, and 53% were hypermethylated (Table 3). Two of these CG sites, FER cg143 at chr5:102025097 and CHCHD2 cg217 at chr7:56174103, were significantly hypomethylated (FER cg143Δ = −4.4) and hypermethylated (CHCHD2 cg217Δ = 1.7) at the predicted CTCF binding sequence (Figures 4A,B). The DNA sequence of the CTCF motif was predicted with the CTCF prediction tool by Ziebarth et al. (2012). The different motif sequences, or PWMs, found within the enhancer of these two genes are listed in Figures 4A,B. The score of each PWM corresponds to the log-odds of the observed sequence being specific rather than randomly generated. Two genes, PARK2 and UBOX5, were significantly hypomethylated (PARK2Δ = −1.3) and hypermethylated (UBOX5Δ = 0.33) across the entire CTCF binding region with p < 0.05 but did not remain significant after multiple hypothesis testing (FDR > 0.05) (Figures 4C,D). Collectively, we conclude that methylation of CTCF binding motifs was vulnerable to rotenone exposure.
Figure 4. Differential methylation within CTCF motifs at Parkinson’s disease genes. (A,B) Two genes, FER and CHCHD2, had differential methylation at CG sites within their predicted CTCF binding motif. The amplified region at both genes covered the first exon and intron. The highlighted blue CG sites are all with significant differential methylation (>|1%|; q-value < 0.05). Red text emphasizes a common SNP. Orange text represents CTCF binding sites in the human genome identified by ENCODE. The bold orange text is the sequence motif predicted by the Cui Lab CTCF prediction tool (Ziebarth et al., 2012). The output of the CTCF binding prediction tool is listed in the table with the name of the position weight matrix motif, motif sequence, motif length, strand orientation, and the integrated output score. (C,D) Two genes, PARK2 and UBOX5, had a significantly different methylation across the region. Each dot represents a CG within the region. Delta indicates the change in the mean CpG methylation percentage and the associated p-value from Fisher’s exact test.
Figure 5. CTCF site H3K27 enhancer activation in response to rotenone. (A) The local abundance of H3K27ac within CTCF binding sites at Parkinson’s disease-associated genes was measured with ChIP-qPCR and expressed as the percent of total DNA input used for immunoprecipitation. (B) The negative control for ChIP analysis was Rabbit IgG. Significance was tested with paired Student’s t-test using the percent input of vehicle (DMSO) vs rotenone, and post hoc analysis for multiple hypotheses was done using the false discovery method. *FDR < 0.05.
Figure 6. CTCF binding at PD-associated genes in response to rotenone. (A) The local abundance of CTCF binding at Parkinson’s disease-associated genes was measured with ChIP-qPCR and expressed as the percent of total DNA input used for immunoprecipitation. (B) The negative control for ChIP analysis was Rabbit IgG. Significance was tested with paired student’s t-test using the percent input of vehicle (DMSO) vs rotenone and post hoc analysis for multiple hypotheses was done using the false discovery method. *FDR < 0.05.
Table 3. Differentially methylated CG sites at Parkinson’s disease-associated genes.
Rotenone Alters H3K27ac at CTCF Binding Sites in Regulatory Regions of PD-Associated Genes
We used ChIP-qPCR to test whether local H3K27ac enrichment overlapped CTCF binding sites in PD-associated genes. Our ChIP-qPCR results demonstrate that four genes (GPRIN3, UBOX5, FER, and BMP4) had significantly increased H3K27ac at CTCF binding motifs with FDR < 0.05. One gene, CNKSR3, had reduced H3K27ac at its CTCF binding motif but was not statistically significant (p = 0.07; FDR = 1) (Figure 5A). H3K27ac enhancer activity was correlated with gene upregulation in three genes (GPRIN3, UBOX5, BMP4). FER was downregulated despite increased H3K27ac in its first intron but may be related to changes in DNA methylation within that region, altered posttranscriptional regulation, or limitations of targeted ChIP analysis. Interestingly, the H3K27 region amplified in qPCR overlapped at least one differentially methylated cytosine for all four significantly enhanced genes (Table 4). Only one of the genes without H3K27ac enrichment, CHCHD2, also had differentially methylated cytosines within the amplified region. These genes had both increased and decreased changes in percent methylation.
Table 4. Differentially methylated CG sites within H3K27ac-enriched regions.
Rotenone Increases CTCF Binding at Three PD-Associated Genes
To determine whether altered DNA methylation and H3K27ac patterns would affect CTCF binding, we measured CTCF enrichment at its binding motif at seven PD-associated genetic loci. CTCF binding was increased at three genes (PARK2, GPRIN3, and BMP4) (Figure 6A). BMP4 had one hypomethylated CG and increased H3K27ac within our selected region. There was an increase in CTCF binding and mRNA expression. PARK2 had two hypomethylated CGs but no increase in H3K27ac in its CTCF binding domain. In this region, CTCF binding increased and mRNA expression decreased. GPRIN3, unlike the other two genes, had more hypermethylated CGs within its CTCF binding motif but the closest CG to its consensus sequence was also hypomethylated. There was increased H3K27ac enrichment at GPRIN3 and increased CTCF binding. GPRIN3 mRNA was significantly upregulated in response to rotenone. These data suggest that both DNA methylation and H3K27ac influence CTCF transcription factor binding and impact the expression of PD-associated genes (Table 5).
Table 5. Correlation of gene expression changes with CTCF binding at PD-associated genes.
Discussion
We selected rotenone based on its ability to model gene–environment interactions in rodents and non-mammalian models of PD (Cannon and Greenamyre, 2013; Johnson and Bobrovskaya, 2015). It is estimated that chronic exposure to concentrations of approximately 20–30 nM of rotenone is enough to cause degeneration of dopaminergic neurons in the midbrain (Greenamyre et al., 2003). While rotenone is a potent inhibitor of complex I in the electron transport chain of mitochondria, rotenone has been shown to cause neurodegeneration by mechanisms unrelated to its effect on complex I (Sherer et al., 2007; Choi et al., 2008). The transcriptome and its regulation have become a focus for understanding these mechanisms outside of the electron transport chain (Cabeza-Arvelaiz and Schiestl, 2012). We observed large-scale changes in gene expression profiles, and many of these genes were enriched in processes involved in gene regulation and chromatin organization (Figure 1). The pathway analysis of DEGs also revealed a large involvement in major intracellular transcription factor pathways. For instance, SMAD proteins are critical for transducing signals from the transforming growth factor (TGFβ) receptors at the plasma membrane which are essential for midbrain dopaminergic survival (Hegarty et al., 2014). Notch signaling is known to have an important role in regulating genes involved in nervous system development and synaptic plasticity (Ables et al., 2011). Lastly, the TP53 pathway is perhaps the most well-known of the toxicant-induced signaling mechanisms to control cell cycle progression and cellular survival. It is thus a critical regulator of programmed cell death in PD and rotenone-induced neurotoxicity (Venderova and Park, 2012).
To determine if the gene expression changes discussed above were due to changes in global levels of epigenetic patterns, we performed dot blots and Western blots to examine global levels of DNA methylation (5mC) and H3K27ac (Figures 3A,B). DNA methylation is the best-studied epigenetic modification, and small changes in methylation at regulatory regions of the genome can have substantial effects on genome integrity during aging (Horvath, 2013). As seen with other pesticide models, we observed a global decrease in total 5mC in response to rotenone. We also chose to look at global H3K27ac levels because it is tightly correlated with gene expression and vulnerable to environmentally driven enhancer activation (Wang et al., 2008). We observed a significant increase in H3K27ac across the genome. H3K27ac is not only an important mark to distinguish poised from active enhancers in bivalent chromatin but also a critical epigenetic modulator in post-mitotic neurons (Maze et al., 2015).
We searched the DEGs for non-coding risk variants associated with PD (fc = fold change). We discovered 13 genes with significant association to PD using GWAS meta-data. We focused on eight genes (ITGA8 0.6 fc, GPRIN3 1.5 fc, FER 0.6 fc, CNKSR3 1.6 fc, BMP4 1.5 fc, UBOX5 1.6 fc, PARK2 0.6 fc, and CHCHD2 1.9 fc) that remained significantly associated in at least two studies with their most significant variant lying in a non-coding region (PDGene; Lill et al., 2012). UBOX5 was the most significantly associated variant according to GWAS meta-data. Furthermore, UBOX5 was the only identified gene with its most significant non-coding variant having known interactions with regulatory elements such as CTCF (Regulomedb; Boyle et al., 2012).
UBOX5 is predicted to have a role in the ubiquitin proteasome system, a well-known PD pathway involved in protein quality control and cellular detoxification (McNaught and Jenner, 2001; McNaught et al., 2003). This pathway is involved in the function of multiple PD-associated genes most notably PARK2 which encodes a ubiquitin ligase. Mutations in PARK2 account for approximately 50% of familial early-onset PD, but the frequency of these mutations decreases with age (Bekris et al., 2010). These mutations generally occur at exon sequences, but other less penetrant but significantly associated polymorphisms with higher frequency in the population occur more often at intronic or regulatory sequences of the gene.
Each of the selected genes had a CTCF binding site determined by ENCODE and the CTCF prediction tool (ENCODE Project Consortium, 2012; Ziebarth et al., 2012). We used the online tool, Regulomedb, to investigate whether SNPs in these CTCF binding regions had evidence of an interaction with CTCF (Table 2). There was evidence of a CTCF interaction in the brain in two of the selected genes GPRIN3 and FER. This rank score generated by Regulomedb is indicative of the strength of the evidence for this interaction with one being the greatest.
CCCTC-binding factor binds regions with allele-specific methylation and preferentially binds the unmethylated allele (Wang et al., 2018). The methylation status of these allele-specific methylated regions is critical to functional CTCF binding and can explain as much as 41% of its variability (Wang et al., 2012). We have previously identified allele-specific methylated regions in the human genome and verified their sensitivity to rotenone exposure (Martos et al., 2017; Freeman et al., 2020). Therefore, we hypothesized that CTCF binding sites at PD-associated genes would also be vulnerable to rotenone. Out of 165 CG sites that met minimum coverage requirements and overlapped predicted CTCF binding motifs, we detected 45 differentially methylated cytosines (Table 3). In two of the genes, the cytosines within the CTCF consensus sequence were differentially methylated but not in any consistent direction (FER-hypomethylated; CHCHD2-hypermethylated) (Figures 4A,B). We saw a similar trend in PARK2 and UBOX5 which had differential methylation across the whole binding region but did not change in a consistent direction (PARK2-hypomethylated) and (UBOX5-hypermethylated) (Figures 4C,D). Overall, there was a slight increase in hypermethylated cytosines (53%) indicating a potential decrease in CTCF binding capacity.
One of the primary functions of CTCF is to act as an insulator by blocking enhancer–promoter interactions (Phillips and Corces, 2009). CTCF is thus tightly correlated with enhancer activity, and its interaction with active enhancers topologically is much greater than with silent regions of the genome (Ren et al., 2017). Histone acetylation patterns also determine chromatin structure, and the histone mark H3K27ac is correlated with active enhancer regions (Wang et al., 2008; Creyghton et al., 2010). p300 is one of the primary histone acetyltransferase enzymes and loads the acetyl group onto the lysine tail of histone 3 at active regions. CTCF binding sites are often located next to at least one p300 binding site and interacts with p300 at the chromatin with active acetylation (Ren et al., 2017).
Histone acetylation patterns are vulnerable to environmental factors and like DNA methylation are heritable (Chinnusamy and Zhu, 2009; Dai and Wang, 2014; Zhu et al., 2018). It is likely that histone acetylation patterns also contribute to the role of genetic variants in disease pathogenesis. We tested H3K27ac levels at CTCF binding sites to determine if acetylation patterns were also sensitive to rotenone at PD-associated genes. We saw an increase in H3K27ac at four of the eight identified genes suggesting strong chromatin interactions with CTCF (FDR < 0.05) (Figure 5A). Notably, all four genes with H3K27ac also overlapped a differentially methylated cytosine within the CTCF binding region (Table 4). This suggests a cross talk mechanism with DNA methylation patterns and H3K27ac enrichment at CTCF binding sites to control chromosomal organization and SNP impacted gene expression.
We observed increased CTCF binding at three differentially expressed PD-associated genes (Figure 6A). PARK2 is a well-known genetic factor in PD as described earlier. Increased CTCF binding at its upstream enhancer decreased its expression, thereby affecting its role in the ubiquitin proteasome system. BMP4 is a gene that encodes bone morphogenetic protein 4, and it regulates neurite outgrowth and axonal transport through the activation of the TGFβ/Smad pathway which is disrupted according to RNA sequencing reactome enrichment data (Supplementary Table 4). Increased CTCF binding was associated with increased BMP4 promoter expression which is essential for dopaminergic neuron differentiation and survival (Table 5; Hegarty et al., 2014). GPRIN3 encodes a protein involved in microtubule dynamics and neurite outgrowth which are both impaired in rotenone-induced neurotoxicity (Cabeza-Arvelaiz and Schiestl, 2012). Interestingly, increased CTCF binding occurred at GPRIN3 within an active transcription start site in the substantia nigra (Table 5; Boyle et al., 2012). It is also associated with a fully penetrant PD mutation causing a triplication of these loci and doubling of GPRIN3 mRNA transcripts (Devine et al., 2011). This observed increase in mRNA transcripts in a mature dopaminergic neuron differentiated from a patient was comparable to the observed 1.5-fold change increase in our rotenone-treated cells.
Data Availability Statement
RNA sequencing data has been submitted to the Gene Expression Omnibus (GEO) repository at NCBI Databank with record number GSE147617. Bisulfite-amplicon sequencing data is publicly available on the sequence read archive (SRA) at NCBI Databank with Bioproject ID PRJNA615220. All other data used to generate this report are available upon request from the corresponding author ZW (zwang47@jhu.edu).
Author Contributions
DF designed and performed experiments with supervision from ZW. Both authors analyzed and interpreted data and prepared the manuscript.
Funding
This project was made possible with support from the JHU Catalyst Award (ZW) and the student Frederik Bang Award (to DF). We acknowledge the National Institutes of Health T32 Pre-doctoral grant (5T32ES7141) for training support (to DF). Experimental disposals and salaries were supported by the U.S. National Institutes of Health (R01ES25761, U01ES026721 Opportunity Fund, and R21ES028351) to ZW.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We acknowledge JHU Center for Alternatives to Animal Testing for their expertise in cell culture techniques and treatments.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00627/full#supplementary-material
Footnotes
- ^ https://github.com/dpryan79/MethylDackel
- ^ http://www.ncbi.nlm.nih.gov/SNP
- ^ http://hek293genome.org/v2/
References
Ables, J. L., Breunig, J. J., Eisch, A. J., and Rakic, P. (2011). Not (ch) just development: notch signalling in the adult brain. Nat. Rev. Neurosci. 12:269. doi: 10.1038/nrn3024
Afgan, E., Baker, D., Van den Beek, M., Blankenberg, D., Bouvier, D., Èech, M., et al. (2016). The galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2016 update. Nucleic Acids Res. 44, W3–W10.
Akalin, A., Kormaksson, M., Li, S., Garrett-Bakelman, F. E., Figueroa, M. E., Melnick, A., et al. (2012). methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 13:R87.
Anders, S., Pyl, P. T., and Huber, W. (2015). HTSeq—a python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169. doi: 10.1093/bioinformatics/btu638
Bekris, L. M., Mata, I. F., and Zabetian, C. P. (2010). The genetics of parkinson disease. J. Geriatr. Psychiatry. Neurol. 23, 228–242.
Boyle, A. P., Hong, E. L., Hariharan, M., Cheng, Y., Schaub, M. A., Kasowski, M., et al. (2012). Annotation of functional variation in personal genomes using RegulomeDB. Genome Res. 22, 1790–1797. doi: 10.1101/gr.137323.112
Buniello, A., MacArthur, J. A. L., Cerezo, M., Harris, L. W., Hayhurst, J., Malangone, C., et al. (2018). The NHGRI-EBI GWAS catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res. 47, D1005–D1012.
Cabeza-Arvelaiz, Y., and Schiestl, R. H. (2012). Transcriptome analysis of a rotenone model of Parkinsonism reveals complex I-tied and-untied toxicity mechanisms common to neurodegenerative diseases. PLoS One 7:e44700. doi: 10.1371/journal.pone.0044700
Cannon, J. R., and Greenamyre, J. T. (2011). The role of environmental exposures in neurodegeneration and neurodegenerative diseases. Toxicol. Sci. 124, 225–250. doi: 10.1093/toxsci/kfr239
Cannon, J. R., and Greenamyre, J. T. (2013). Gene-environment interactions in Parkinson’s disease: Specific evidence in humans and mammalian models. Neurobiol. Dis. 57, 38–46. doi: 10.1016/j.nbd.2012.06.025
Chinnusamy, V., and Zhu, J. (2009). Epigenetic regulation of stress responses in plants. Curr. Opin. Plant Biol. 12, 133–139. doi: 10.1016/j.pbi.2008.12.006
Choi, W. S., Kruse, S. E., Palmiter, R. D., and Xia, Z. (2008). Mitochondrial complex I inhibition is not required for dopaminergic neuron death induced by rotenone, MPP+, or paraquat. Proc. Natl. Acad. Sci. U.S.A. 105, 15136–15141. doi: 10.1073/pnas.0807581105
Cicchetti, F., Drouin-Ouellet, J., and Gross, R. E. (2009). Environmental toxins and Parkinson’s disease: what have we learned from pesticide-induced animal models? Trends Pharmacol. Sci. 30, 475–483. doi: 10.1016/j.tips.2009.06.005
Coetzee, S. G., Pierce, S., Brundin, P., Brundin, L., Hazelett, D. J., and Coetzee, G. A. (2016). Enrichment of risk SNPs in regulatory regions implicate diverse tissues in Parkinson’s disease etiology. Sci. Rep. 6:30509.
Creyghton, M. P., Cheng, A. W., Welstead, G. G., Kooistra, T., Carey, B. W., Steine, E. J., et al. (2010). Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc. Natl. Acad. Sci. U.S.A. 107, 21931–21936. doi: 10.1073/pnas.1016071107
Dai, H., and Wang, Z. (2014). Histone modification patterns and their responses to environment. Curr. Environ. Health Rep. 1, 11–21. doi: 10.1007/s40572-013-0008-2
Dawson, T. M., Mandir, A. S., and Lee, M. K. (2002). Animal models of PD: pieces of the same puzzle? Neuron 35, 219–222.
de Lau, L. M., and Breteler, M. M. (2006). Epidemiology of Parkinson’s disease. Lancet Neurol. 5, 525–535.
Devine, M. J., Ryten, M., Vodicka, P., Thomson, A. J., Burdon, T., Houlden, H., et al. (2011). Parkinson’s disease induced pluripotent stem cells with triplication of the α-synuclein locus. Nat. Commun. 2, 1–10.
Do, C., Lang, C. F., Lin, J., Darbary, H., Krupska, I., Gaba, A., et al. (2016). Mechanisms and disease associations of haplotype-dependent allele-specific DNA methylation. Am. J. Hum. Genet. 98, 934–955. doi: 10.1016/j.ajhg.2016.03.027
Eden, E., Navon, R., Steinfeld, I., Lipson, D., and Yakhini, Z. (2009). GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. BMC Bioinform. 10:48. doi: 10.1186/1471-2105-10-48
ENCODE Project Consortium (2012). An integrated encyclopedia of DNA elements in the human genome. Nature 489:57. doi: 10.1038/nature11247
Fabregat, A., Korninger, F., Viteri, G., Sidiropoulos, K., Marin-Garcia, P., Ping, P., et al. (2018). Reactome graph database: Efficient access to complex pathway data. PLoS Comput. Biol. 14:e1005968. doi: 10.1371/journal.pcbi.1005968
Feng, Y., Liu, T., Dong, S.-Y., Guo, Y.-J., Jankovic, J., Xu, H., et al. (2015). Rotenone affects p53 transcriptional activity and apoptosis via targeting SIRT 1 and H3K9 acetylation in SH-SY 5Y cells. J. Neurochem. 134, 668–676. doi: 10.1111/jnc.13172
Franco, R., Li, S., Rodriguez-Rocha, H., Burns, M., and Panayiotidis, M. I. (2010). Molecular mechanisms of pesticide-induced neurotoxicity: relevance to Parkinson’s disease. Chem. Biol. Interact. 188, 289–300. doi: 10.1016/j.cbi.2010.06.003
Freeman, D. M., Lou, D., Li, Y., Martos, S. N., and Wang, Z. (2020). The conserved DNMT1-dependent methylation regions in human cells are vulnerable to neurotoxicant rotenone exposure. Epigenet. Chrom. 13:17. doi: 10.1186/s13072-020-00338-8
Greenamyre, J. T., Betarbet, R., and Sherer, T. B. (2003). The rotenone model of Parkinson’s disease: genes, environment and mitochondria. Parkinson. Relat. Disord. 9, 59–64. doi: 10.1016/s1353-8020(03)00023-3
Harrison, I. F., Smith, A. D., and Dexter, D. T. (2018). Pathological histone acetylation in Parkinson’s disease: Neuroprotection and inhibition of microglial activation through SIRT 2 inhibition. Neurosci. Lett. 666, 48–57. doi: 10.1016/j.neulet.2017.12.037
Hegarty, S. V., Collins, L. M., Gavin, A. M., Roche, S. L., Wyatt, S. L., Sullivan, A. M., et al. (2014). Canonical BMP–Smad signalling promotes neurite growth in rat midbrain dopaminergic neurons. Neuromol. Med. 16, 473–489. doi: 10.1007/s12017-014-8299-5
Horvath, S. (2013). DNA methylation age of human tissues and cell types. Genome. Biol. 14:R115. doi: 10.1186/gb-2013-14-10-r115
Hou, L., Zhang, X., Wang, D., and Baccarelli, A. (2012). Environmental chemical exposures and human epigenetics. Int. J. Epidemiol. 41 79–105. doi: 10.1093/ije/dyr154
Huang, M., Lou, D., Charli, A., Kong, D., Jin, H., Anantharam, V., et al. (2019). Mitochondrial dysfunction induces epigenetic dysregulation by h3k27 hyperacetylation to perturb active enhancers in Parkinson’s disease models. bioRxiv [Preprint]. doi: 10.1101/808246
Johnson, M. E., and Bobrovskaya, L. (2015). An update on the rotenone models of Parkinson’s disease: their ability to reproduce the features of clinical disease and model gene–environment interactions. Neurotoxicology 46, 101–116. doi: 10.1016/j.neuro.2014.12.002
Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14:R36.
Krueger, F. (2015). Trim Galore. A Wrapper Tool Around Cutadapt And FastQC to Consistently Apply Quality And Adapter Trimming to FastQ Files. Available online at: https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/ (accessed March 6, 2019).
Labbé, C., Lorenzo-Betancor, O., and Ross, O. A. (2016). Epigenetic regulation in Parkinson’s disease. Acta Neuropathol. 132, 515–530.
Lee, Y. C., Lai, C. Q., Ordovas, J. M., and Parnell, L. D. (2011). A database of gene-environment interactions pertaining to blood lipid traits, cardiovascular disease and type 2 diabetes. J. Data Min. Genom. Proteom. 2:2602. doi: 10.4172/2153-2602
Li, L., and Dahiya, R. (2002). MethPrimer: designing primers for methylation PCRs. Bioinformatics 18, 1427–1431. doi: 10.1093/bioinformatics/18.11.1427
Lill, C. M., Roehr, J. T., McQueen, M. B., Kavvoura, F. K., Bagade, S., Schjeide, B. M., et al. (2012). Comprehensive research synopsis and systematic meta-analyses in Parkinson’s disease genetics: the PDGene database. PLoS Genet. 8:e1002548. doi: 10.1371/journal.pone.1002548
Lin, Y., Boone, M., Meuris, L., Lemmens, I., Van Roy, N., Soete, A., et al. (2014). Genome dynamics of the human embryonic kidney 293 lineage in response to cell biology manipulations. Nat. Commun. 5, 1–12.
Martos, S. N., Li, T., Ramos, R. B., Lou, D., Dai, H., Xu, J., et al. (2017). Two approaches reveal a new paradigm of ‘switchable or genetics-influenced allele-specific DNA methylation’with potential in human disease. Cell Discov. 3:17038.
Maurano, M. T., Humbert, R., Rynes, E., Thurman, R. E., Haugen, E., Wang, H., et al. (2012). Systematic localization of common disease-associated variation in regulatory DNA. Science 337, 1190–1195. doi: 10.1126/science.1222794
Maze, I., Wenderski, W., Noh, K., Bagot, R. C., Tzavaras, N., Purushothaman, I., et al. (2015). Critical role of histone turnover in neuronal transcription and plasticity. Neuron 87, 77–94. doi: 10.1016/j.neuron.2015.06.014
McNaught, K. S. P., Belizaire, R., Isacson, O., Jenner, P., and Olanow, C. W. (2003). Altered proteasomal function in sporadic parkinson’s disease. Exp. Neurol. 179, 38–46. doi: 10.1006/exnr.2002.8050
McNaught, K. S. P., and Jenner, P. (2001). Proteasomal function is impaired in substantia nigra in parkinson’s disease. Neurosci. Lett. 297, 191–194. doi: 10.1016/s0304-3940(00)01701-8
Nalls, M. A., Pankratz, N., Lill, C. M., Do, C. B., Hernandez, D. G., Saad, M., et al. (2014). Large-scale meta-analysis of genome-wide association data identifies six new risk loci for parkinson’s disease. Nat. Genet. 46:989.
Nica, A. C., Montgomery, S. B., Dimas, A. S., Stranger, B. E., Beazley, C., Barroso, I., et al. (2010). Candidate causal regulatory effects by integration of expression QTLs with complex trait genetic associations. PLoS Genet. 6:e1000895. doi: 10.1371/journal.pone.1000895
Ning, Y., Yang, J., Ma, G., and Chen, P. (2011). Modelling rock blasting considering explosion gas penetration using discontinuous deformation analysis. Rock Mech. Rock Eng. 44, 483–490. doi: 10.1007/s00603-010-0132-3
Ong, C., and Corces, V. G. (2014). CTCF: an architectural protein bridging genome topology and function. Nat. Rev. Genet. 15:234. doi: 10.1038/nrg3663
Orth, M., Tabrizi, S., Schapira, A., and Cooper, J. (2003). α-Synuclein expression in HEK293 cells enhances the mitochondrial sensitivity to rotenone. Neurosci. Lett. 351, 29–32. doi: 10.1016/s0304-3940(03)00941-8
Otasek, D., Morris, J. H., Bouças, J., Pico, A. R., and Demchak, B. (2019). Cytoscape automation: empowering workflow-based network analysis. Genome Biol. 20, 1–15.
Park, G., Tan, J., Garcia, G., Kang, Y., Salvesen, G., and Zhang, Z. (2016). Regulation of histone acetylation by autophagy in parkinson disease. J. Biol. Chem. 291, 3531–3540. doi: 10.1074/jbc.M115.675488
Pedersen, B. S., Eyring, K., De, S., Yang, I. V., and Schwartz, D. A. (2014). Fast and Accurate Alignment of Long Bisulfite-Seq Reads. arXiv [preprint] arXiv:1401.1129.
Phillips, J. E., and Corces, V. G. (2009). CTCF: master weaver of the genome. Cell 137, 1194–1211. doi: 10.1016/j.cell.2009.06.001
Ren, G., Jin, W., Cui, K., Rodrigez, J., Hu, G., Zhang, Z., et al. (2017). CTCF-mediated enhancer-promoter interaction is a critical regulator of cell-to-cell variation of gene expression. Mol. Cell. 67, 1049–1058.
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Sadowski, M., Kraft, A., Szalaj, P., Wlasnowolski, M., Tang, Z., Ruan, Y., et al. (2019). Spatial chromatin architecture alteration by structural variations in human genomes at the population scale. Genome Biol. 20:148.
Sherer, T. B., Richardson, J. R., Testa, C. M., Seo, B. B., Panov, A. V., Yagi, T., et al. (2007). Mechanism of toxicity of pesticides acting at complex I: relevance to environmental etiologies of Parkinson’s disease. J. Neurochem. 100, 1469–1479.
Stepanenko, A., and Dmitrenko, V. (2015). HEK293 in cell biology and cancer research: Phenotype, karyotype, tumorigenicity, and stress-induced genome-phenotype evolution. Gene 569, 182–190. doi: 10.1016/j.gene.2015.05.065
Supek, F., Bošnjak, M., Škunca, N., and Šmuc, T. (2011). REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS One 6:e21800. doi: 10.1371/journal.pone.0021800
Tang, Z., Luo, O. J., Li, X., Zheng, M., Zhu, J. J., Szalaj, P., et al. (2015). CTCF-mediated human 3D genome architecture reveals chromatin topology for transcription. Cell 163, 1611–1627. doi: 10.1016/j.cell.2015.11.024
Teixeira, J., Basit, F., Swarts, H. G., Forkink, M., Oliveira, P. J., Willems, P. H., et al. (2018). Extracellular acidification induces ROS-and mPTP-mediated death in HEK293 cells. Redox Biol. 15, 394–404. doi: 10.1016/j.redox.2017.12.018
United States Environmental Protection Agency [USEPA] (2007). Reregistration Eligibility Decision for Rotenone (Case No. 0255). Available at: https://archive.epa.gov/ (accessed May 5, 2020).
Venderova, K., and Park, D. S. (2012). Programmed cell death in parkinson’s disease. Cold Spring Harb. Perspect. Med. 2:a009365. doi: 10.1101/cshperspect.a009365
Verstraeten, A., Theuns, J., and Van Broeckhoven, C. (2015). Progress in unraveling the genetic etiology of parkinson disease in a genomic era. Trends Genet. 31, 140–149. doi: 10.1016/j.tig.2015.01.004
Wang, H., Dong, X., Liu, Z., Zhu, S., Liu, H., Fan, W., et al. (2018). Resveratrol suppresses rotenone-induced neurotoxicity through activation of SIRT1/Akt1 signaling pathway. Anat. Rec. 301, 1115–1125. doi: 10.1002/ar.23781
Wang, H., Lou, D., and Wang, Z. (2019). Crosstalk of genetic variants, allele-specific DNA methylation, and environmental factors for complex disease risk. Front. Genet. 9:695. doi: 10.3389/fgene.2018.00695
Wang, H., Maurano, M. T., Qu, H., Varley, K. E., Gertz, J., Pauli, F., et al. (2012). Widespread plasticity in CTCF occupancy linked to DNA methylation. Genome Res. 22, 1680–1688. doi: 10.1101/gr.136101.111
Wang, Z., Zang, C., Rosenfeld, J. A., Schones, D. E., Barski, A., Cuddapah, S., et al. (2008). Combinatorial patterns of histone acetylations and methylations in the human genome. Nat. Genet. 40:897. doi: 10.1038/ng.154
Ye, J., Coulouris, G., Zaretskaya, I., Cutcutache, I., Rozen, S., and Madden, T. L. (2012). Primer-BLAST: a tool to design target-specific primers for polymerase chain reaction. BMC Bioinform. 13:134. doi: 10.1186/1471-2105-10-134
Zhu, Y., Li, Y., Lou, D., Gao, Y., Yu, J., Kong, D., et al. (2018). Sodium arsenite exposure inhibits histone acetyltransferase p300 for attenuating H3K27ac at enhancers in mouse embryonic fibroblast cells. Toxicol. Appl. Pharmacol. 357, 70–79. doi: 10.1016/j.taap.2018.08.011
Keywords: CTCF, Parkinson’s disease, single nucleotide polymorphisms, rotenone, DNA methylation, histone modifications
Citation: Freeman DM and Wang Z (2020) Epigenetic Vulnerability of Insulator CTCF Motifs at Parkinson’s Disease-Associated Genes in Response to Neurotoxicant Rotenone. Front. Genet. 11:627. doi: 10.3389/fgene.2020.00627
Received: 28 March 2020; Accepted: 26 May 2020;
Published: 07 July 2020.
Edited by:
Yanqiang Li, Houston Methodist Research Institute, United StatesCopyright © 2020 Freeman and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhibin Wang, zwang47@jhu.edu