Nuclear and Mitochondrial DNA Methylation Patterns Induced by Valproic Acid in Human Hepatocytes.

Valproic acid (VPA) is one of the most widely prescribed antiepileptic drugs in the world. Despite its pharmacological importance, it may cause liver toxicity and steatosis through mitochondrial dysfunction. The aim of this study is to further investigate VPA-induced mechanisms of steatosis by analyzing changes in patterns of methylation in nuclear DNA (nDNA) and mitochondrial DNA (mtDNA). Therefore, primary human hepatocytes (PHHs) were exposed to an incubation concentration of VPA that was shown to cause steatosis without inducing overt cytotoxicity. VPA was administered daily for 5 days, and this was followed by a 3 day washout (WO). Methylated DNA regions (DMRs) were identified by using the methylated DNA immunoprecipitation-sequencing (MeDIP-seq) method. The nDNA DMRs after VPA treatment could indeed be classified into oxidative stress- and steatosis-related pathways. In particular, networks of the steatosis-related gene EP300 provided novel insight into the mechanisms of toxicity induced by VPA treatment. Furthermore, we suggest that VPA induces a crosstalk between nDNA hypermethylation and mtDNA hypomethylation that plays a role in oxidative stress and steatosis development. Although most VPA-induced methylation patterns appeared reversible upon terminating VPA treatment, 31 nDNA DMRs (including 5 zinc finger protein genes) remained persistent after the WO period. Overall, we have shown that MeDIP-seq analysis is highly informative in disclosing novel mechanisms of VPA-induced toxicity in PHHs. Our results thus provide a prototype for the novel generation of interesting methylation biomarkers for repeated dose liver toxicity in vitro.

Introduction

In recent years, regulatory toxicology has adopted the concept of adverse outcome pathways, which are defined as causal linkages between a series of key events and the induction of a certain end point of toxicity.[1] This has further increased the demand for better understanding underlying mechanisms of action, and this is facilitated by the application of the so-called omics technologies, which have enabled global investigations of genomic entities in response to toxic challenges as well as in vitro models for predicting human toxicity.[2,3]

Within this context, the epigenome has recently gained considerable interest, in particular because epigenomic changes induced by toxicants may be persistent and even inheritable.[4] Thus, it may be hypothesized that epigenomic modifications present relevant biomarkers for repeated dose toxicity, in particular in in vitro models.

Among multiple epigenomic processes, DNA methylation refers to the conversion of a cytosine residue into 5-methylcytosine (5mC) by the enzyme DNA methyltransferase (DNMT), which uses S-adenosyl methionine (SAM) as the methyl donor.[5−8] DNA hypermethylation is expected to result in gene silencing.[9,10] Interestingly, oxidative stress, which is considered a causal factor in a multitude of diseases and drug-induced organ injuries, has been shown to play a role in epigenetic processes such as DNA methylation.[11]

Recently, it has been debated whether not only nuclear DNA (nDNA) but also mitochondrial DNA (mtDNA) may be subjected to epigenetic modifications (e.g., methylation).[12−14] Such mitochondrial epigenetic modifications in mtDNA may modulate nDNA, whereas nuclear epigenetic modifications may affect mtDNA as well.[13,14] Overall, whereas drug-induced oxidative stress is known to play a role in liver toxicity,[15] a possible epigenetic mechanism linking oxidative stress and changes in patterns of 5mC in nDNA and mtDNA to drug-induced liver injury has not been investigated yet. Therefore, we decided to investigate changes in DNA methylation patterns induced by the prototypical liver toxicant valproic acid (VPA) in primary human hepatocytes (PHHs).

VPA is one of the most widely prescribed antiepileptic drugs in the world. Despite its pharmacological importance, it may cause liver toxicity, particularly steatosis via an abnormal retention of lipids within the hepatocyte.[16,17] VPA is a simple fatty acid that is almost entirely metabolized by the liver.[17] Treatment with VPA has been shown to lead to a significant increase in the production of reactive oxygen species (ROS) in isolated rat liver mitochondria,[16] rat hepatocytes,[18] and mouse embryos.[19] It is thus hypothesized that VPA exposure induces hepatotoxicity via oxidative stress, possibly mediated by the overproduction of ROS and consequent mitochondrial dysfunction.[16,18,20] Interestingly, VPA may also affect the epigenome: prolonged VPA treatment (3 days) induces a decrease in 5-hydroxymethylcytosine (5hmC) in isolated mitochondria of mouse fibroblast 3T3-L1 cells.[21] In addition, Dong et al. showed that VPA induces DNA hypomethylation in the nuclear extracts of adult mouse brains.[22] However, the effect of VPA on the epigenome in human liver cells is not well understood. Therefore, this study focuses on VPA-induced changes in methylation patterns in nDNA and mtDNA. We furthermore investigated whether, in particular, mtDNA epigenetic modifications are persistent upon terminating VPA treatment.

For this, PHHs were exposed daily for 5 days to a subcytotoxic dose of VPA, after which DNA was isolated. In order to evaluate the potential persistence of induced methylome changes after terminating the VPA treatment, some of the cells were kept in medium for 3 days. Identification of methylated regions in the DNA was carried out by means of the methylated DNA immunoprecipitation–sequencing (MeDIP-seq) method. Data analysis comprised (1) quality control, (2) mapping the reads against the UCSC human genome 19 (hg19), and (3) selecting significant differentially methylated regions (DMRs) according to the protocol described in the MEDIPS Tutorial.[23−26] DMR-derived genes were fed into ConsensusPathdB for biological interpretation and into Cytoscape for network visualization.

Experimental Procedures

Cell Culture and Treatment

Cryopreserved PHHs (Invitrogen) from three individuals (Hu8084, Hu4197, and Hu4227) were thawed for 1 min at 37 °C in a water bath. Next, these PHHs were pooled in order to bypass interindividual variability in susceptibility to toxicants and cultured in 6-well plates in a collagen sandwich,[27] according to the supplier’s protocol (Invitrogen). In brief, three vials (one per individual) were thawed in 50 mL of CHRM medium (CHRM CM7000, Invitrogen) and centrifuged for 10 min at 100g at 4 °C. After removing the supernatant, the cell pellet was suspended in CHPM medium (CHPM CM9000, Invitrogen) at a concentration of 1.0 × 106 cells/mL. Pooled PHHs were cultured in collagen-precoated 6-well plates (Gibco) (density: 2.0 × 106 cells/well), and cells were attached to the wells in the incubator for 4 h at 37 °C. After the incubation, the debris was removed by shaking and washing the cells twice with William’s medium E (WME CM6000, Invitrogen). Subsequently, cells were covered with 150 μL/well of a 1.0 mg/mL collagen mixture (containing 10× DMEM, 1.0 mg/mL Collagen type 1 (BD), 0.2 M NaOH, and MQ) and incubated at 37 °C for approximately 30 min until the collagen was fixed. After that, culture medium (WME + 20 mL of Hepatocyte Supplement Pack (CM4000, Invitrogen) substituted with 1% penicillin/streptomycin (Gibco)) was added.

After 3 days, the PHHs were exposed daily to VPA or 1% EtOH (vehicle control) in culture medium for 5 days. A screening experiment examining an incubation concentration range of 0–30 mM VPA revealed that 15 mM VPA created a diffusion gradient across the collagen layer sufficiently high to induce observable steatosis in PHHs, as confirmed by clear fat accumulation visualized by means of standardized BODIPY 493/503 staining after 24 h of VPA treatment (Figure S1). Simultaneously, this incubation concentration of 15 mM VPA did not yield any notable signs of PHH cytotoxicity, as evaluated by MTT assay, after 24 or 48 h. Culture medium was changed daily, thereby administering a new incubation concentration of VPA or EtOH to the cells. After the exposure period of 5 days, PHHs were lysed for DNA isolation. Another well of PHHs was maintained in culture for 3 consecutive days without VPA treatment (from now on called the washout (WO) period); the culture medium was changed every day. All experiments were performed in biological triplicates.

DNA Isolation

PHHs were collected in 500 μL of digestion buffer (1 mM EDTA, 50 mM Tris-HCl, pH 8.0, 5% SDS) and proteinase K (1 mg/mL) (Ambion). After incubation for 1 h at 55 °C, the proteinase K was inactivated at 80 °C. Treatment with RNase A (400 μg/mL) (Qiagen) and 1% collagenase (Sigma) was performed for 1 h at 37 °C. An equal amount of phenol–chloroform–isoamylalcohol (PCI; 25:24:1) (Sigma) was added, and the samples were shaken manually for 5 min. After centrifugation, the upper phase was again treated with PCI until protein was no longer visible at the interphase. The upper phase was precipitated using 50 μL of 3 M NaOAc pH 5.6 and 1250 μL of cold 100% ethanol. The DNA pellet was washed using cold 70% EtOH, dissolved in 50 μL of nuclease free water, and quantified spectrophotometrically using a NanoDrop 1000 (Thermo Scientific, Waltham, MA). The total amount of DNA obtained was at least 10 μg, the 260/280 ratio was between 1.7 and 1.9, and the 260/230 ratio was higher than 1.6.

MeDIP-Seq Protocol

MeDIP-seq was performed, with all the biological triplicates after DNA isolation, according to the protocol of Taiwo et al.,[28] with minor adjustments.

DNA Fragmentation to a Size of ∼200 bp

For DNA fragmentation, 300 ng of isolated DNA was sonicated with a bioruptor (Diagenode) using instrument settings of 15 cycles, each consisting of 30 s on/off periods. After fragmentation, the genomic DNA size range was assessed using an Agilent 2100 Bioanalyzer and high-sensitivity DNA chips (Agilent Technologies), according to the manufacturer’s instructions.

Library Preparation and Size Selection

Libraries were prepared using 300 ng of fragmented DNA (∼200 bp) and the NEBNext ultra DNA library prep kit for Illumina (NEB), according to the manufacturer’s protocol. In detail, 300 ng of fragmented DNA was mixed with 3.0 μL of EndPrep enzyme mix and 6.5 μL of 10× end repair reaction buffer. The mixture was incubated in a thermocycler for 30 min at 20 °C and subsequently for 30 min at 65 °C.

Adaptors were ligated by adding 15 μL of blunt/TA ligase master mix (NEB), 2.5 μL of NEBNext adaptor for Illumina (NEB), and 1 μL of ligation enhancer (NEB) to the end prep reaction mixture. These samples were handled according to the manufacturer’s instructions for adaptor ligation in the NEBNext ultra DNA library prep kit for Illumina (NEB). Thereafter, adaptor-ligated DNA was size-selected using size selection of adaptor-ligated DNA in the NEBNext ultra DNA library prep kit for Illumina (NEB). In short, first, AMPure XP beads were used to remove the unwanted large DNA fragments. Second, AMPure XP beads were used to bind the DNA fragments of interest. The DNA fragments from the beads were diluted into 28 μL of 10 mM Tris-HCl.

MeDIP Analysis

The purified adaptor-ligated DNAs were used for methylated DNA immunoprecipitation (MeDIP), according to the manufacturer’s instructions from the MagMeDIP (Diagenode) and IPure (Diagenode) kits. In brief, IP incubation mix (containing magbuffer A, magbuffer B, 1.5 μL of methylated DNA positive control, and 1.5 μL of nonmethylated DNA negative control) was added to the 300 ng of purified adaptor-ligated DNA samples, denaturated at 95 °C for 3 min, and quickly chilled on ice; 7.5 μL (10%) was kept aside as an InPut sample (IP sample) and stored at 4 °C. Beads were washed three times using ice-cold diluted magbuffer A (11 μL of magbeads was needed per IP sample) and collected in diluted magbuffer A (20 μL per IP sample). The remaining sample, 90% of the 300 ng denatured and purified adaptor-ligated DNA sample, was immunoprecipitated (MeDIP-sample) overnight on a rotating wheel at 4 °C using 5 μL of the monoclonal antibody against 5-methylcytosine (MagMeDIP, Diagenode) and 20 μL of magbeads (MagMeDIP, Diagenode).

Next, the MeDIP samples were washed twice with wash buffer (WB) 1 and once with WB2 (MagMeDIP, Diagenode) and kept on ice. DNA was eluted two times in 50 μL of elution buffer (EB) (IPure, Diagenode) for MeDIP samples and once in 92.5 μL of EB for IP samples (every sample has at this moment a volume of 100 μL). DNA from all samples was bound to the magnetic beads (IPure, Diagenode) by incubation for 1 h on a rotating wheel at room temperature. Thereafter, bound DNA was washed once using WB1 and once using WB2 (IPure, Diagenode), eluted in buffer C (IPure, Diagenode), and precipitated with ethanol (IPure, Diagenode). The DNA pellets were washed using 70% EtOH, air-dried, and resuspended in 30 μL of DNase-free water.

Quality Control

Quantitative PCR (qPCR) was used for controlling DNA methylation enrichment. qPCR was performed by measuring the cycles to threshold (Ct) values of 1 μL of purified DNA sample and 24 μL of qPCR mixture (1 μL of the provided primer pair (reverse and forward), 12.5 μL of SYBR Green PCR master mix, and 10.5 μL of water) using the following temperature profile: 95 °C for 7 min, 40 cycles at 95 °C for 15 s and 60 °C for 1 min, followed by 1 min at 95 °C. The efficiency of MeDIP was calculated by performing qPCR and using the following formula: %(meDNA-IP/total input) = 2[( × 100%. The efficiency for methylated DNA fragments was good (>50%) for all samples. More interestingly, the efficiency for nonmethylated DNA fragments was overall lower than 1.0%. This implies that only methylated DNA fragments have been selected in our MeDIP assay and that nonmethylated DNA fragments have been discarded in the process.

PCR Amplification and Size Selection

PCR was used to amplify the MeDIP adaptor-ligated DNA fragments. In brief, 25 μL of NEBNext high fidelity 2× PCR master mix (NEB), 1 μL of index primer (NEB) that was used as a barcode for each sample (Table S1), and 1 μL of universal PCR primer (NEB) were added to 23 μL of the MeDIP adaptor-ligated DNA fragments. PCR was performed by using the following temperature profile: 98 °C for 30 s, 15 cycles at 98 °C for 10 s, 65 °C for 30 s, and 72 °C for 30 s, followed by 5 min at 72 °C and hold at 4 °C as described before.[28]

Thereafter, PCR-amplified DNAs (libraries) were cleaned using cleanup of PCR amplification in the NEBNext ultra DNA library prep kit for Illumina (NEB). Fragmented DNA size and quality were checked using the Agilent 2100 Bioanalyzer and high-sensitivity DNA chips (Agilent Technologies). In addition, generated libraries were size-selected on a 2% TAE low melting point (LMP) agarose gel; fragments of 250–350 bp were excised and the MinElute gel DNA extraction kit (Qiagen) was used to extract and purify the DNA libraries. Libraries were quantified on a Qubit fluorimeter (Invitrogen) by using the Qubit dsDNA HS assay kit (Invitrogen). All kits and chips were used according to the manufacturer’s protocol.

Sequencing

The 12 amplified libraries, each sample having its own index primer, were pooled at an equimolar concentration of 2 nM, based on Qubit measurements. From the 2 nM stock solution, 10, 15, and 20 pM were then loaded onto three separated channels of a 1.4 mm flow cell (Illumina) and cluster amplification was performed on a cBot (Illumina). Clusters were generated on cBot (Illumina) using the TruSeq PE cluster kit V3, according to the manufacturer’s instructions (Illumina), and the paired-end libraries were sequenced using 2 × 100 cycles TruSeq SBS kit v3 paired-end by sequencing by synthesis (SBS) on the Illumina HiSeq 2000. Base calling was performed by using Casava 1.8.2 (Illumina), and demultiplexing was performed using bcl2fastq 1.8.4 (Illumina). Sequence reads were aligned against the human reference genome called UCSC hg19. This alignment results in FASTQ files for each barcoded library. MeDIP-seq raw data are available on arrayExpress (accession number: E-MTAB-4437).

Data Analysis

MeDIP-Seq Analysis

FastQC was applied to check the quality of the 100 bp reads pairs of the 12 sequenced samples. Paired-end sequencing reads were aligned against hg19 using Bowtie2 software. The MEDIPS package (version 1.16.0, Bioconductor) was used for the analysis of the MeDIP-seq data.[26,29,30] The default parameters described in the MEDIPS guideline (version 1.16.0)[23] were applied to all data from individual chromosomes, including mitochondrial DNA (chrM). The data set was divided into four different groups of triplicates: (1) Co_MeDIP includes the sequencing data of PHHs exposed daily for 5 days to the vehicle control; (2) VPA_MeDIP includes the sequencing data of PHHs exposed daily for 5 days to VPA; (3) CoWO_MeDIP contains the sample exposed daily for 5 days to the vehicle control followed by a WO period of 3 days; and (4) VPAWO_MeDIP includes the sequencing data of PHHs treated daily with VPA for 5 days followed by a WO period of 3 days. Statistical analysis was performed by applying the default parameters of MEDIPS, using the edgeR module (empirical analysis of digital gene expression data in R) that uses Bayes estimation and exact tests based on the negative binomial distribution.[31] Notably, raw count data was normalized using the weighted trimmed mean of M-values (TMM normalization). Regions were considered significantly methylated if the edgeR p-value was below 0.01 and if the number of reads, of a specific region, in one of the samples was higher than the mean of reads of all regions (the whole genome), which is the background correction. This p-value was derived from other studies where MeDIP-seq analysis was performed.[32−34]

In addition, copy number variations were the same between control samples and treated samples, since we used the same pool of donors and PHHs, which are nongrowing cells. Annotation of DMRs into different genomic locations was achieved by using HOMER software, mapping each DMR to a 5′-untranslated region (5′UTR), 3′-untranslated region (3′UTR), exon, intergenic region, intron, noncoding region, promoter transcription start site (promoter-TSS), or transcription termination site (TTS). Furthermore, each DMR was annotated into a gene name also by using HOMER. Regions were merged if (1) the start of a region was consecutive to the end of the previous region and (2) if the HOMER annotations of these consecutive DMRs were the same.

Pathway Analysis

Significant selected DMR lists and unique gene lists were uploaded onto VENNY.[35] ConsensusPathDB[36] was used to identify and visualize the involvement of the unique genes in biological processes that may be affected on the level of pathways by selecting significant pathways with a p-value < 0.01 from a gene enrichment analysis. The liver steatosis adverse outcome pathway of Landesmann et al.,[37] describing the key events in liver steatosis development like the accumulation of triglycerides in the liver,[37,38] was used to select steatosis-related pathways from the ConsensusPathDB results.

Network Visualization

Methylated genes were then uploaded to Cytoscape. A circular layout was selected, and the network was analyzed as undirected. Fold change (FC) data were added, and nodes were colored (green = hypermethylation (positive FC) and red = hypomethylation (negative FC)).

The first neighbors of methylated hub genes were selected by using the tool first neighbors of selected nodes in Cytoscape. Then, a submolecular induced epigenome network with its first neighbors was prepared in Cytoscape.

Results

Methylation Changes in nDNA

Daily VPA treatment of PHHs for 5 days resulted in 3006 hypermethylated nDNA DMRs and 3077 hypomethylated nDNA DMRs (Figure A). This comprises 0.02% of the ∼31 million possible methylated 100 bp nDNA regions. After VPA treatment, DMRs annotated in nDNA exons were found to be significantly more hypermethylated than hypomethylated (278 DMRs vs 48 DMRs) (Figure B). Furthermore, in introns, noncoding regions, and promoter-TSS, hypomethylation of nDNA occurred more frequently than hypermethylation (Figure B). We observed no difference between the number of hyper- and hypomethylated nDNA DMRs at other genomic locations: 3′UTR, 5′UTR, intergenic, and TTS (Figure B). The relative number of DMRs (the density of DMRs per chromosome divided by the number of bases per chromosome) was determined to obtain a global overview of the methylation status of nDNA per chromosome. After VPA treatment, methylation changes in nDNA were present in all chromosomes but occurred relatively more frequently in chr07, chr16, chr17, chr19, chr20, and chr22 (Figure ).

Figure 1

(A) Number of significantly hyper- and hypomethylated nDNA DMRs after 5 days of daily administration of VPA to PHHs. Differentially hyper- and hypomethylated regions were obtained using the edgeR approach and mapped to regions using human genome 19. (B) Histogram of genomic locations after 5 days of daily treatment of PHHs with VPA. Significant hyper- and hypomethylated nDNA DMRs (p ≤ 0.01) are shown in different genomic regions: 3′UTR (3′), 5′UTR (5′), exon, intergenic, intron, noncoding, promoter-TSS, and TTS, as defined and mapped by HOMER using human genome 19. PHH, primary human hepatocytes; DMRs, differentially methylated regions; VPA, valproic acid; 5mC, 5-methylcytosine; TSS, transcription start site; TTS, transcription termination site; nDNA, nuclear DNA.

Figure 2

Histogram of the distribution per chromosome of the “relative number” of hypermethylated and hypomethylated nDNA DMRs after 5 days of daily exposure of PHHs to VPA. The relative number of DMRs was measured by dividing the observed number of significantly hypo- or hypermethylated DMRs per chromosome by the number of bases per chromosome. Chromosomes that show more 5mC modifications than general are indicated with an asterisk (*). PHHs, primary human hepatocytes; DMRs, differentially methylated regions; VPA, valproic acid; 5mC, 5-methylcytosine; nDNA, nuclear DNA.

The promoter-TSS nDNA DMRs (38 hyper- and 68 hypomethylated) were classified into one significantly enriched pathway, namely, semaphorin interactions. Furthermore, the exon nDNA DMRs (278 hyper- and 48 hypomethylated) were classified into two significantly enriched pathways: (1) FTO obesity variant mechanism and (2) cell junction organization. HOMER-annotated genes of the 3006 hypermethylated and 3077 hypomethylated nDNA DMRs were classified into 119 significantly enriched pathways, in which the steatosis-related pathway “regulation of lipid metabolism by peroxisome proliferator-activated receptor alpha (PPARalpha/PPARA)” was present in the top 10 of the list (Table S2). The overlap of differentially methylated genes that were present in all 119 significantly enriched pathways after VPA treatment formed a complex network module (Figure S2A,B). The hypomethylated gene 1994 (ELAVL1) demonstrated the largest number of neighbors (i.e., number of connections within the whole network module), namely, 65. ELAVL1 (FC = −3.76 and p-value = 0.0065) is involved in the control of generic gene expression. Second in line is 2033 (EP300 (FC = −5.2 and p-value = 0.00099)) with 33 neighbors (Figure S2B,C and Table S3). EP300 is of major interest since it regulates histone acetyltransferase (HAT) and its neighbors could be classified into significantly enriched pathways that comprise steatosis-related pathways (Table S3). Here, we have identified strongly hypomethylated genes that play a role in lipid regulation (such as CREB binding protein (CREBBP), retinoid X receptor alpha (RXRA), and nuclear receptor corepressor 1 (NCOR1)).

Additionally, we found 11 hypomethylated nDNA DMRs and 20 hypermethylated nDNA DMRs that persisted after the WO period (Table S4).

Methylation Changes in mtDNA

After 5 days of daily VPA administration, we observed significant hypomethylation of 5 DMRs in mtDNA, whereas hypermethylation of DMRs in mtDNA could not be demonstrated. After HOMER annotation of these DMRs, we identified 7 hypomethylated genes that are involved in protein synthesis during translation (such as mitochondrially encoded tRNA alanine (MT-TA), mitochondrially encoded tRNA asparagine (MT-TN), mitochondrially encoded tRNA cysteine (MT-TC), mitochondrially encoded tRNA serine 1 (MT-TS1), and mitochondrially encoded tRNA aspartic acid (MT-TD)) or that affect the third enzyme of the electron transport chain, cytochrome c oxidase (such as mitochondrially encoded cytochrome c oxidase I (MT-CO1) and mitochondrially encoded cytochrome c oxidase II (MT-CO2)). The top 2 mtDNA genes, with the lowest p-value, were MT-CO1 and MT-TS1.

After the WO period, none of these mtDNA DMRs remained significantly hypomethylated.

Discussion

Insights vary with regard to whether not only nDNA but also mtDNA may be subjected to epigenetic modifications such as DNA methylation.[12−14] This is of interest as recent work demonstrated that drug-induced DNA methylation changes can contribute to the development of hepatotoxicity.[39] Therefore, in this study, we investigated hepatotoxicant-induced changes in patterns of methylation in nDNA and mtDNA, as well as the involvement of these changes in mechanisms of steatosis, using VPA as a prototypical steatotic compound; we also evaluated whether VPA-induced methylome changes were persistent upon terminating VPA treatment.

For analyzing VPA-induced epigenomic alterations, we applied MeDIP-seq. Several methods for genome-wide methylation profiling have been extensively compared in the past, namely, MeDIP and methylCAP, as enrichment-based methods with two bisulfite-based methods (RRBS and Infinium),[27,40] as well as MeDIP and MBD assays with the bisulfite-based methods (MethylC and RRBS).[41] All methods appear to provide accurate DNA methylation detection, with a very high correlation between MeDIP- and bisulfite-based methods, whereas the MeDIP enrichment method has the advantage of being the more cost-effective approach.[27,40,41]

Methylation Changes in nDNA and mtDNA of PHHs Exposed to VPA

Methylation Changes in nDNA

It has been shown that hypermethylation of the promoter-TSS and the first exon lead to the strongest gene silencing.[42] In VPA-treated PHHs, 38 hypermethylated and 68 hypomethylated promoter-TSS DMRs were found. In addition, 278 significantly hypermethylated exon nDNA DMRs, which contain 31 first exon regions, and 48 hypomethylated exon nDNA DMRs, which contain 2 first exon regions, were found. Pathway analysis of HOMER-derived genes only from promoter-TSS DMRs or only from exon DMRs did not identify any pathway that could be linked to steatosis. The semaphorin interactions pathway was enriched by the promoter-TSS DMRs, whereas the exon DMRs were classified into the “FTO obesity variant mechanism” and “cell junction organization” pathways. Interestingly, upon evaluating all DMRs, we found an important steatosis-related pathway to be modulated, namely, “regulation of lipid metabolism by peroxisome proliferator-activated receptor alpha (PPARalpha/PPARA)” (Table S2). Therefore, we suggest that methylation changes in genomic regions other than promoter-TSS and exon regions are also important for VPA-induced toxicity. The exact function and role of methylation changes in the body of a gene (introns, intergenic, noncoding, TTS, 3′UTR, and 5′UTR regions) are not clear yet, but these have been suggested to play a role in splicing.[43]

In the network (Figure S2A), which was built by using HOMER-derived genes from DMRs of all regions, we found EP300 (strongly hypomethylated, chromosome 22:41535301–41535400, intron) interacting with the most neighbors linked to steatosis-related pathways, analyzed with ConsensusPathDB (Table S3 and Figure S2C). This is an interesting finding because the encoded protein, histone acetyltransferase (HAT) p300, is known to be involved in the induction of steatosis: the amount of p300-C/EBPα/β complexes were increased in the livers of mice that developed steatosis by activating five promotor genes that drive triglyceride synthesis.[44] VPA itself is known as a histone deacetylase (HDAC) inhibitor.[45] Reduced histone acetylation is associated with decreased DNA methylation.[46] The VPA-induced hypomethylation of EP300 may result in the induced expression of EP300 and the increased protein expression of p300, which may thus contribute to the VPA-induced development of steatosis in PHHs. However, Lamparter et al.[47] have shown a decrease in protein expression of p300 after 24 h of VPA treatment in mouse P19 embryonic cells.[47] Furthermore, we found CREBBP significantly hypomethylated, and it is known that the protein encoded by this gene, besides regulating lipid metabolism by peroxisome proliferator-activated receptor alpha, has intrinsic HAT activity.[48] This suggests that both the VPA-induced hypomethylation of CREBBP and EP300 increases HAT activity.

Overall, these epigenetic changes in the nDNA of PHHs after VPA treatment seem to be largely reversible. We identified 31 nDNA DMRs (11 hypomethylated nDNA DMRs and 20 hypermethylated nDNA DMRs) that persisted after the WO period. Among these, 5 zinc finger protein genes (3 hypermethylated and 2 hypomethylated) were methylated and point toward a persistent impact on a wide range of generic cell functions (Table S4).

Methylation Changes in mtDNA

Iacobazzi et al. and Ghosh et al. suggested that not only nDNA but also mtDNA may be subjected to epigenetic modifications related to disease development, environmental exposure, drug treatment, and aging.[13,49] We were able to confirm that after 5 days of VPA treatment a general hypomethylation of mtDNA in PHHs was observed. Specifically, we found 7 hypomethylated mtDNA genes. In particular, genes encoding for cytochrome c oxidase (such as MT-CO1 and MT-CO2) were hypomethylated. Cytochrome c oxidase is the final enzyme of the electron transport chain (ETC) and the proton pump of complex III in the ETC.[50] It is known that complexes I and III of the ETC are the major sites of ROS formation.[51−54] Bayir et al. showed that cytochrome c oxidase dysregulation induces the release and production of superoxide radicals (increased ROS production).[54] It has been shown that increased ROS production, in our case as a consequence of VPA metabolism, induces oxidative stress and has been hypothesized to play a role in the formation of liver toxicity (such as NAFLD) and therefore also in steatosis.[18,20,55] This is in line with observations in other studies.[56−58]

After the WO period, none of these mtDNA DMRs remained significantly hypomethylated.

The mitochondrial D-loop has been considered the third strand of mtDNA (next to the heavy (H) and light (L) strands); mtDNA synthesis starts from here.[59,60] Moreover, the D-loop plays a central role in mtDNA transcription and replication.[61] Therefore, it is likely that methylation in this D-loop region influences mtDNA gene expression.[61] However, we did not find any significant DMRs in this D-loop region. Apart from D-loop methylation, methylation of other subregions of the mtDNA, namely, the gene bodies, may have a profound effect on gene expression.[61] It has been shown that gene body methylation in mtDNA correlates with gene expression and the progression of nonalcoholic steatohepatitis.[61,62] More interestingly, van der Wijst et al.[60] have shown that induction of methylation by M.CviPI (mitochondria-targeted DNA methyltransferases, MLS-M.CviPI: Nys-1 chorella virus GpCme methyltransferase) decreased the expression of certain mitochondrial genes (such as MT-ND1, MT-CO1, MT-CYB, and MT-RNR2), which were measured by RT-PCR.[60] Furthermore, methylation in the MT-CO2 gene was also reported.[60] In line with van der Wijst et al.,[60] we have found methylation changes in the gene bodies of the genes MT-CO1 and MT-CO2 after 5 days of treatment.

Hypothetical nDNA and mtDNA Crosstalk of PHHs Exposed to VPA

In nDNA of PHHs exposed daily for 5 days to VPA, we found hypermethylated DMRs that were annotated by HOMER to MAT and DNMT. This hypermethylation of MAT may result in a repressed expression of SAM, which is used as the methyl donor for DNA methylation. This may decrease the 5mC content (hypomethylation).[5−8] Furthermore, hypermethylation of DNMT may lead to a decreased 5mC content.[5−8,63] We therefore suggest that the hypermethylation of MAT and DNMT in nDNA, which results in decreased 5mC content (hypomethylation), underlies the observed hypomethylation of mtDNA (genes such as MT-TA, MT-TN, MT-TC, MT-TS1, MT-TD, MT-CO1, and MT-CO2). After the WO period, this crosstalk between nDNA and mtDNA does not seem to be persistent.

Conclusions

The present study demonstrates that relatively short treatments with xenobiotic compounds can actually impact the human epigenome in vitro, changing methylation patterns in PHHs of all genomic regions in both nDNA and mtDNA and thereby contributing to the induction of steatosis. In particular, EP300 seems to be a promising candidate biomarker for the development of steatosis in PHHs in vitro. This hypothesis needs to be validated in the future, in which multiple hepatotoxic compounds should be considered. Interestingly, most nDNA DMRs, all mtDNA DMRs, and the suggested VPA-induced crosstalk between mtDNA and nDNA appear to be reversible upon terminating VPA treatment. Overall, our study has demonstrated that MeDIP-seq analysis of VPA-treated PHHs is highly informative in disclosing novel mechanisms of VPA-induced toxicity. Future research should focus on seeking confirmation of the relevance of the observed epigenomic events on the global transcriptome and proteome levels and on validating the hypothetical crosstalk between the methylation changes in nDNA and mtDNA.