Genome-wide analysis of DNA methylation identifies novel differentially methylated regions associated with lipid accumulation improved by ethanol extracts of Allium tubersosum and Capsella bursa-pastoris in a cell model

Hepatic steatosis is the most common chronic liver disease in Western countries. Both genetic and environmental factors are known as causes of the disease although their underlying mechanisms have not been fully understood. This study investigated the association of DNA methylation with oleic acid-induced hepatic steatosis. It also examined effects of food components on DNA methylation in hepatic steatosis. Genome-wide DNA methylation of oleic acid (OA)-induced lipid accumulation in vitro cell model was investigated using reduced representation bisulfite sequencing. Changes of DNA methylation were also analyzed after treatment with food components decreasing OA-induced lipid accumulation in the model. We identified total 81 regions that were hypermethylated by OA but hypomethylated by food components or vice versa. We determined the expression of seven genes proximally located at the selected differentially methylated regions. Expression levels of WDR27, GNAS, DOK7, MCF2L, PRKG1, and CMYA5 were significantly different between control vs OA and OA vs treatment with food components. We demonstrated that DNA methylation was associated with expression of genes in the model of hepatic steatosis. We also found that food components reversely changed DNA methylation induced by OA and alleviated lipid accumulation. These results suggest that DNA methylation is one of the mechanisms causing the hepatic steatosis and its regulation by food components provides insights that may prevent or alleviate lipid accumulation.

Introduction Nonalcoholic fatty liver disease (NAFLD) is a chronic liver disease caused by fat accumulation in the liver due to imbalance between triglyceride (TG) acquisition and removal without alcohol consumption [1]. Progress of NAFLD ranges from simple hepatic steatosis to non-alcoholic steatohepatitis (NASH), fibrosis, and even hepatic cancer. NAFLD is associated with obesity, dyslipidemia, and insulin resistance, which are also known as characteristics of metabolic syndromes [2]. Although the pathogenesis of NAFLD is not fully understood, it has been shown that hepatic de novo lipogenesis is increased by activation of lipogenic factors such as SREBP-1c, PPARγ, and fatty acid synthase (FASN) [2][3][4][5]. Subsequently, accumulation of free fatty acids (FFAs) in the liver causes lipotoxicity and oxidative stress, which lead to hepatocyte injury and progress to NASH and fibrosis [2][3][4]. It is of interest that dietary factors affect de novo hepatic lipogenesis via the crucial factors FASN and PPARγ, and can thereby mitigate NAFLD and obesity, based on a cell and an animal model [6,7]. However, underlying mechanisms of the regulation have not been clearly elucidated.
Substantial emerging evidence has demonstrated that the development and progression of NAFLD is regulated by epigenetic mechanisms including DNA methylation, histone modification and non-coding RNAs [8][9][10][11][12][13]. In addition, it was reported that both DNA methylation and histone modification are regulated by dietary factors in animal models [14,15]. However, the relevance between DNA methylation and histone modification has not been clearly elucidated.
Over the past three decades, it has been shown that various dietary factors, including methyl donors, protein, polyunsaturated fatty acid, sugar, and bioactive components, modulate epigenetic status and affect gene expression in various cell and animal models of human diseases including NAFLD [16,17].
In animal models of fatty liver, deficiency of methyl-donors such as betaine, choline, and folate affects one-carbon metabolism, and consequently progression to NASH [18,19]. In high-fat-sucrose diet-induced obesogenic mice, dietary methyl-donor supplements improved fatty liver by regulating DNA methylation of FASN and its expression [20]. Consistent with this, Chang et al. showed that berberine modulated DNA methylation of the promoter of microsomal triglyceride transfer protein, which is a key gene in lipid homeostasis [21]. Lingonberries prevent hepatic steatosis through regulation of DNA methylation of genes associated with inflammation and lipid synthesis in a high-fat diet-induced animal model [22]. These show that not only methyl-donors but also dietary components affect DNA methylation in an animal model of fatty liver.
Modification of histones by dietary components is also involved in prevention and/or attenuation of NAFLD. In previous study, it has been demonstrated that hepatic steatosis was improved through inhibition of histone acetylation by extract of Allium tuberosum (EAT) containing sulfur and phenolic compounds [23]. In addition, extract of Capsella bursa-pastoris (ECB) containing flavonoids decreased lipid accumulation through inhibition of histone acetyltransferase in an in vitro cell model [24].These results suggest that dietary components including EAT and ECB may be applicable for reducing lipid accumulation and improving hepatic steatosis. These studies also showed that 200-400 μg/mL EAT or ECB affect lipid accumulation and epigenetic status in HepG2 cells without toxic effects. However, little is known about global effects on DNA methylation by treatment with dietary components EAT or ECB in hepatic steatosis.
In this study, we performed reduced representation bisulfite sequencing (RRBS) to investigate changes in genome-wide DNA methylation by EAT and ECB in an OA-induced hepatic

Reduced representation bisulfite sequencing (RRBS) library preparation and sequencing
To construct RRBS libraries with MspI and ApeKI, 500 ng of input genomic DNA in 50 μl was digested with MspI (NEB, Ipswich, MA, USA) at 37˚C for 7 h. ApeKI (NEB) was then added and incubation was continued at 75˚C for 16-20 h. The digested products were purified with a MiniElute PCR Purification Kit (Qiagen, Venlo, Netherlands). After purification, dA was added to the digested products with blunt-ended ligation, followed by ligation of methylatedadapter. A slice of the 160-420 bp fraction was excised from 2% agarose gel. Bisulfite conversion was conducted using a ZYMO EZ DNA Methylation-Gold Kit (ZYMO Research, Irvine, CA, USA) following the manufacturer's instructions. The final libraries were generated by PCR amplification using PfuTurbo Cx Hotstart DNA polymerase (Agilent technologies, Santa Clara, CA, USA). RRBS libraries were analyzed by an Agilent 2100 Bioanalyzer (Agilent Technologies). Before sequencing the samples, the quantity of sequenceable library fragments was determined via qPCR. Samples were then diluted to 10 nM with elution buffer (QIAGEN

RRBS data analysis
We performed FastQC v0.11.2 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to control the quality of raw reads, and trimmed adaptor sequencing using trim galore v0.4.1 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Trimmed sequences were aligned to the human reference genome (hg19) using BS-seeker2 v2.0.10 (Guo et al., 2013) with Bowtie2. We built double enzyme MspI (CCGG) and ApeKI (GCWGC) fragments with length range 30-500 bp in silico to cover MspI and ApeKI fragments of RRBS libraries. We aligned the reads with Bowtie2 in local alignment mode allowing four mismatches per read. Unmapped reads were remapped in paired-end mode to improve mapping rates. Where two paired-end mates overlapped, we called methylation levels of each CpG site after removing one mate.
To avoid low mapping efficiency due to adapter contamination in the sequencing data, size selection (160-420 bp) was performed. It was found that mappability (>70%, S2 Table) and depth (>40 x, S2 Table) were better than those observed in previous studies, although these studies analyzed different cells and tissues [25,26]. This suggested that our sample preparation, generation of DNA methylomes, processes of sequencing, and mapping analysis had no critical problems. However, physical coverage could not be calculated in this analysis because C to T is the most common substitution (~65%) in all single nucleotide polymorphisms (SNPs) and could not be distinguished from C to T conversion by bisulfite treatment [25]. In general, less than four million CpGs out of 29 million in the genome were physically covered by our RRBS screening [26].

Differentially methylated region (DMR) analysis
We used a custom Perl script to identify DMRs (100 bp) between groups. Briefly, DNA methylation levels on the genome were profiled by sliding a fixed-size window (100 bp) in 50 bp increments through the reference genome (hg19). DNA methylation ratios (0 to 1) of all CpG sites in a given window were compared between two groups (control vs OA, OA+EAT vs OA and OA+ECB vs OA) using the Mann-Whitney U test (p < 0.01). To filter out unreliable DMR candidates, regions covered by less than 10 reads or showing mean difference of < 0.2 between groups were discarded. Identified DMRs were annotated using HOMER (v5.7) with the UCSC reference gene annotation (hg19).

Western blot analysis
HepG2 cells were harvested and homogenized in a cell lysis buffer (Cell Signaling Technology, Beverly, MA, USA) containing a Xpert phosphatase and protease inhibitor cocktail solution (GenDEPOT, Barker, TX, USA). Lysates were centrifuged at 10,000 g for 15 min at 4˚C. Total cellular proteins (20 μg) were loaded on SDS-PAGE and transferred onto nitrocellulose (NC) membranes (GE Healthcare Life Science, Pittsburgh, PA, USA). Blocking buffer contained 5% skim milk in TBS-T at room temperature. Blots were incubated with primary antibody against FASN (Cell Signaling Technology, Beverly, MA, USA) overnight at 4˚C. Secondary antibody conjugated with horseradish peroxidase was complexed with primary antibody and developed with an ECL detection kit (DoGEN, Seoul, Korea).

Quantification of gene expression using real-time PCR
RNA was extract from the treated cells with an RNeasy Mini kit (Qiagen), according to the manufacturer's instructions. A total of 500 ng RNA was reverse-transcribed with reverse transcriptase (TOYOBO, Osaka, Japan) at 30˚C for 10 min, 42˚C for 20 min, and 99˚C for 5 min. Relative quantification of gene expression was determined with the cDNA and primers listed in S1 Table. The reaction was carried out using SYBR green super mix (TOYOBO) and a thermal cycler (Bio-Rad, Hercules, CA, USA). Amplification conditions consisted of 40 cycles of 95˚C for 10 sec, 58˚C for 10 sec, 72˚C for 20 sec, and a final melting curve step.

Statistical analysis
All results were shown as the mean ± S.D. Statistical significances between groups were assessed using unpaired t-tests, using GraphPad Prism 5 Software (San Diego, CA, USA). Statistical significance was accepted at p < 0.05, p < 0.01 and p < 0.001.

Establishment of the cell model of hepatic steatosis and RRBS analysis
Consistent with previous reports showing that oleic acid stimulated lipid accumulation in HepG2 cells and increased expression of FASN [23,24], OA induced lipid accumulation up to 2 times and 200 μg/mL EAT (III) or ECB (IV) decreased the lipid accumulation ( Fig 1A). We assessed cytotoxicity of EAT or ECB in HepG2 using MTT assay. Treatment with EAT or ECB (0, 200, 400, 800 μg/mL) did not induced cytotoxicity in HepG2 cells (S1 Fig). Since it is known that FAS is a lipogenic enzyme which regulates fatty acid synthesis [27], we further examined the beneficial effects of EAT and ECB on protein expression change of FASN. Increased FASN expression by OA was significantly attenuated by treatment with EAT or ECB in the hepatic steatosis model (Fig 1B and S2 Fig). This showed that our hepatic cell model system was adequate for further investigation of the underlying mechanisms of hepatic steatosis.
As digestion of genomic DNA with double restriction enzymes has been found to increase CpG coverage [28], MspI and ApeKI, were used in this study to efficiently determine genomewide DNA methylation. In addition, to decrease the likelihood of false positives, less than ten sequenced reads were excluded from our analysis.
Although emerging data shows that single CpGs can be important in regulation of gene expression [29,30], roles of single CpGs in gene regulation are still debatable. To avoid selection of single differentially methylated CpGs, averages of all CpGs in 100 bp were calculated and statistically analyzed to identify DMRs. Non-CpGs, CHH and CHG, were excluded from the analysis since the mechanisms underlying whether non-CpG methylation plays a role in gene regulation have not been clearly elucidated [31,32].

Genome-wide methylation analysis in a cell model of hepatic steatosis
To investigate the underlying mechanisms of lipid accumulation during hepatic steatosis, global DNA methylation pattern, was analyzed using reduced representation bisulfite sequencing (RRBS). A total of twelve samples (n = 3 for each group) were prepared for RRBS as listed in S2 Table. In total, 517 million reads were sequenced and 378 million of these were mapped to the human reference genome. More than 70% of reads were successfully mapped. Sequencing depth ranged from 42 to 97 reads throughout the reference genome (S2 Table).
To identify significant DMRs, we selected regions where the difference in methylation level was more than 20% between groups (Mann-Whitney test, p < 0.01). In the OA group (II) compared to the control group (I), there was a total of 406 DMRs, including 215 hypermethylated and 191 hypomethylated ( Table 2). In the OA+EAT-treated group (III) compared to the OA group (II), 532 DMRs were identified, with 296 as hypermethylated and 236 hypomethylated. In the OA+ECB group (IV) compared to the OA group, there was a total of 265 DMRs of which 109 were hypermethylated and 156 were hypomethylated. It is of interest that about 60% of the identified significant DMRs were located in CpG islands (18-20%), exons (15-25%), TSSs (11-13%), and 5' UTRs (> 1%), and therefore more likely to be involved with gene expression (Table 3).
We further selected 22 DMRs that were hypermethylated in the OA group (II) compared to the control (I), but hypomethylated in the OA+EAT (III) compared to the OA (II), and 39 DMRs showing the converse methylation pattern in the same group comparison (Fig 2A). In addition, 11 DMRs were hypermethylated in the OA (II) but hypomethylated by the OA+ECB (IV), and nine DMRs conversely methylated between the same groups ( Fig 2B). A total of 81 regions that were hypermethylated in OA (II) but hypomethylated in OA+EAT (III) and OA +ECB (IV) or vice versa were identified to investigate the effects of EAT and ECB on DNA methylation during hepatic steatosis. As shown in Fig 3, it was evident that the selected DMRs between groups (IIvs I, III vs II, and IV vs II) were clearly clustered and methylation levels were significantly different.

Selection of putative genes proximally located at the selected DMRs
A total of 77 putative DMRs, excluding regions not assigned by the HOMER program, were identified as regions affected by EAT or ECB in the cell model of hepatic steatosis (Table 4). Among them, 31 DMRs were hypermethylated by OA while 46 DMRs were hypomethylated by OA. The level of DNA methylation in 72 DMRs was reversely changed by treatment with EAT or ECB while only five DMRs were affected by both EAT and ECB. Interestingly, 37 DMRs were located at functional genomic structures such as TSSs, exons, CpG islands, and introns. In view of the known link between hepatic steatosis and the metabolic syndrome [17], we summarized potentially relevant functions of the genes nearest to selected DMRs (Table 4). Among the annotated genes, 26 were found to be related to the metabolic syndrome, including obesity, diabetes, hypertension, cardiovascular diseases, inflammation, and stroke.

Gene expression affected by modulation of DNA methylation by EAT and ECB
To investigate the effects of DNA methylation on expression of genes proximal to DMRs, we examined the association of expression levels of genes with the selected DMRs in hepatic steatosis. As shown in Fig 4, seven genes (WDR27, GNAS, DOK7, EDN3, MCF2L, PRKG1, and CMYA5) were selected based on genomic location and relevance to metabolic syndrome. Expression of WDR27, GNAS, MCF2L, and PRKG1 was increased by OA but decreased by treatment with EAT or ECB (Fig 4A, 4B, 4E and 4F). Expression of DOK7 and CMYA5 was decreased by OA but increased by EAT (Fig 4C and 4G). Expression of EDN3 was decreased by OA but not changed by EAT or ECB (Fig 4D). DMR location of each gene was marked as thick red line in S3 Fig. Methylation

Selection of DMRs associated with hepatic steatosis
More than 60% of our selected DMRs were located at functional genomic regions such as TSSs, exons, 5' UTRs, CpG islands, and introns, while less than 40% of the non-selected DMRs were located in these regions (Table 3). This suggested that the selected DMRs may be more likely to be involved in regulation of genes expression.
The number of hypermethylated DMRs was not substantially different from that of hypomethylated DMRs in the analysis of OA (II) vs control (I), or in our other comparisons (III vs II and IV vs II) ( Table 4). This suggests that both up-and down-regulation of DNA methylation are involved in lipid accumulation and may stimulate or suppress gene expression, and is consistent with reports showing that some genes (FASN, PPARγ, and SREBP1) are increased but others (SIRT1, FOXO1, and ATGL) are decreased in cell models of hepatic steatosis [24,33].

Effects of dietary components on DNA methylation
Allium tuberosum (AT) and Capsella bursa-pastoris (CB) have been widely consumed as food ingredients in Korea. It has been known that AT exerts various health benefits in inflammation, diabetes, and cardiovascular diseases, as does CB in inflammation and cancer [34,35]. However, their underlying mechanisms are not fully understood. Recently, it was suggested that histone modifications by EAT and ECB may be involved in alleviating hepatic steatosis and provide a therapeutic target for its treatment or prevention [23,24]. This study showed for the first time that EAT and/or ECB reversed DNA methylation induced by OA in an in vitro cell model of hepatic steatosis (Table 4). Many studies demonstrated that Allium tubersosum (AT) and Capsella bursa-pastoris (CB) consist of sulphur-containing compounds, phenolic compounds, acylated flavonol glucosides, flavonoids, organic acid, and other many compounds [36][37][38][39]. Among these compounds, both AT and CB contain same flavonoid compounds such as kaempferol and quercetin. Kaempferol and quercetin are flavonoid compounds having an antioxidant activity. It has been known that they improved NAFLD by reducing hepatic lipid accumulation and oxidative stress [40][41][42]. It was also reported that kaempferol and quercetin induced epigenetic modifications through regulating histone deacetylases (HDACs) and/or DNMTs [43][44][45]. These suggested that AT and CB may improve NAFLD by regulating DNA methylation.
It is generally known that hypermethylated DNA suppresses gene expression while hypomethylation stimulates transcription. Consistent with it, we showed hypomethylation at an intron of MCF2L and an exon of PRKG1 by OA (Table 4), and the hypomethylatoin was associated with increased expression of the genes (Fig 4E and 4F). EAT and ECB induced hypermethylation of MCF2L and PRKG1 and decreased their expression. In addition, hypermethylation at transcription start site of CMYA5 by OA decreased its expression while hypomethylated by EAT increased its expression (Table 4 and Fig 4G). These suggest that level of methylation of MCF2L, PRKG1, and CMYA5 may regulate expression of the genes.
Although the physiological function of WDR27 has not been fully demonstrated, an SNP in intergenic region adjoining WDR27 (rs924043) was associated with type 1 diabetes, which suggests that its expression may be involved in metabolic syndrome [46]. In addition, duplication of WDR27 has been seen in an obese patient, which suggests that WDR27 may be overexpressed in obesity [47]. Consistent with this, WDR27 its expression was significantly increased by OA while decreased by EAT and ECB (Fig 4A). In this study, we showed that an intron of WDR27 was hypermethylated in the OA group but hypomethylated after treatment with EAT and ECB (Table 4). It is important to note that DNA hypermethylation can increase expression of genes although it is generally known that hypermethylation suppresses gene expression. Recently, this was supported by a systematic analysis of binding of 542 transcription factors (TFs) to methylated or unmethylated CpGs [48]. For activation of gene expression by the TFs, 34% and 23% of the TFs preferred hypermethylated and hypomethylated CpGs respectively, while 33% of the TFs did not prefer CpGs. Together, this suggested that DNA hypermethylation can also stimulate gene expression.
It is known that GNAS regulates homeostasis of glucose and energy metabolism [49]. Interestingly, the methylation level of CpG sites located at the upstream of the GNAS TSS was significantly decreased after dietary intervention [50]. Consistent with this, we showed that this TSS region was hypomethylated by ECB. Significantly decreased gene expression was also   observed with ECB treatment (Fig 4B). Together, this suggested that DMRs in the GNAS TSS region were affected by dietary factors and associated with its transcription. It was reported that DOK7 plays a crucial role in the progress of metabolic disease in an animal model through regulation of DNA methylation at its promoter, affecting its expression [51]. We found that exonic and intronic CpG islands in DOK7 were hypomethylated by OA and hypermethylated by EAT (Table 4), and that expression of the gene was decreased by OA and elevated by EAT in our cell model of hepatic steatosis (Fig 4C). Together, these findings also suggest that DNA methylation at CpG islands in DOK7 are regulated by dietary factors and associated with its expression.
Although it is known that genetic variants in a region between GNAS and EDN3 are associated with hypertension and cardiovascular disease [52,53], DNA methylation may not involve in expression of EDN3 in cell model of hepatic steatosis since hypermethylation and hypomethylation at the TSS of EDN3 by OA and ECB, respectively (Table 4), decreased expression of EDN3 (Fig 4D).
Although this study showed that hepatic steatosis in cell model was affected by DNA methylation regulating expression of each gene by dietary factors, it did not exclude the possibility that the selected genes may synergistically play roles in hepatic steatosis. As previously described, NAFLD is caused by multi-factors such as SREBP-1c, PPARγ, and FASN. It was also demonstrated that other factors and different regulatory mechanisms were involved in the progression of hepatic steatosis [54,55,17]. These studies suggest that several factors instead of a factor may synergistically and/or spatiotemporally play roles during hepatic steatosis. Further study will be required to uncover whether all the selected genes exert their functions in a combinational manner during hepatic steatosis.
Since cell line systems do not reflect exact whole organisms such as interactions with other cell types/tissues, metabolic status and effect of hormones etc., it has been still controversial whether the significance of cell line data can be reproduced in in vivo studies. However, cell line systems are very efficient to select or narrow down targets through screening of compounds and will provide information for further studies. This study described regulation of DMRs during steatosis in a cell model and will help further investigate the animal or clinical studies.
In conclusion, this study showed, for the first time, that modulation of DNA methylation is one of the mechanisms during hepatic steatosis in a cell model. This study also showed the regulation of expression of genes by DNA methylation in hepatic steatosis model alleviated by EAT and ECB. The data present here provide a potential lead into further studies investigating hepatic steatosis and may give an insight to development of prevention or treatment of hepatic steatosis.
Supporting information S1