Arachidonic and oleic acid exert distinct effects on the DNA methylome

ABSTRACT

Abnormal fatty acid metabolism and availability are landmarks of metabolic diseases, which in turn are associated with aberrant DNA methylation profiles. To understand the role of fatty acids in disease epigenetics, we sought DNA methylation profiles specifically induced by arachidonic (AA) or oleic acid (OA) in cultured cells and compared those with published profiles of normal and diseased tissues. THP-1 monocytes were stimulated with AA or OA and analyzed using Infinium HumanMethylation450 BeadChip (Illumina) and Human Exon 1.0 ST array (Affymetrix). Data were corroborated in mouse embryonic fibroblasts. Comparisons with publicly available data were conducted by standard bioinformatics. AA and OA elicited a complex response marked by a general DNA hypermethylation and hypomethylation in the 1–200 μM range, respectively, with a maximal differential response at the 100 μM dose. The divergent response to AA and OA was prominent within the gene body of target genes, where it correlated positively with transcription. AA-induced DNA methylation profiles were similar to the corresponding profiles described for palmitic acid, atherosclerosis, diabetes, obesity, and autism, but relatively dissimilar from OA-induced profiles. Furthermore, human atherosclerosis grade-associated DNA methylation profiles were significantly enriched in AA-induced profiles. Biochemical evidence pointed to β-oxidation, PPAR-α, and sirtuin 1 as important mediators of AA-induced DNA methylation changes. In conclusion, AA and OA exert distinct effects on the DNA methylome. The observation that AA may contribute to shape the epigenome of important metabolic diseases, supports and expands current diet-based therapeutic and preventive efforts.

KEYWORDS:

• β-oxidation
• DNA methylation
• epigenomics
• fatty acid
• PPAR
• sirtuin

Abbreviations

450K array	=	Infinium HumanMethylation450 BeadChip
5mdC	=	5-methyl-2′-deoxycytidine
AA	=	arachidonic acid
AACpG	=	AA-induced differentially methylated CpG
CGI	=	CpG island
CpG	=	5′-CG-3′ dinucleotide
DB-CpG	=	CpG scored by Δβ
EPA	=	eicosapentaenoic acid
FACpG	=	differentially methylated CpG common to the AACpG, OACpG and PACpG sets
FDR	=	false discovery rate
MEF	=	mouse embryonic fibroblasts
N-CpG	=	number of CpGs normalized by dividing the number of differentially methylated CpGs by the total number of CpGs represented in the 450K array for a given gene
OA	=	oleic acid
ACpG	=	OA-induced differentially methylated CpG
PA	=	palmitic acid
PACpG	=	PA-induced differentially methylated CpG
PPAR	=	peroxisome proliferator-activated receptor
SIRT	=	sirtuin. Official acronyms were used for all genes

Introduction

An organism's fatty acid (FA) pool reflects both composition of the diet and endogenous synthesis. Furthermore, abnormal FA levels have long been recognized to participate in metabolic diseases such as diabetes.¹ Since metabolic diseases are frequently associated with aberrant DNA methylation profiles,² it is conceivable that FAs are among factors that mediate remodeling of the epigenome in response to dietary composition, the cell's metabolic status and pathogenic signals. In humans, only a handful of studies have investigated a link between DNA methylation and fat intake, focusing on specific genes associated with FA metabolism, inflammation and regulation of circadian rhythms.^3-8 Furthermore, a recent genome-wide study comparing peripheral whole blood methylation profiles to both quantity and quality of dietary fat intake in adolescents, found a larger number of methylation changes associated with the latter.⁹ In that study, DNA methylation profiles associated with polyunsaturated:saturated FA ratio were related to pathways regulated by the peroxisome proliferator-activated receptor α (PPAR-α) and adipogenesis. Conversely, no pathways were identified of DNA methylation profiles associated with the monounsaturated:saturated FA ratio. In addition, numerous mouse model studies have addressed short and long term epigenetic effects of diets containing variable combinations of specific FAs.^8,10

Although diet-based studies are important as they recognize the complexity of the organism's response to nutrients, mechanistic insights can be gained by a simplified complementary approach based on the manipulation of individual FA levels. Accordingly, the effects of short FAs, such as butyrate, on DNA methylation have been long recognized.¹¹ In addition, very low-density lipoprotein (VLDL) elicited a global DNA hypermethylation response that is markedly stronger than the one induced by low- or high-density lipoprotein in cultured human THP-1 macrophages.^12,13 The fact that VLDL is characteristically rich in triglycerides suggests that FAs might be mediators of the epigenetic responses to VLDL. Accordingly, arachidonic acid (AA) induces DNA hypomethylation in human umbilical vein endothelial cells at 3 μM concentration, globally and of angiogenesis regulator kinase insert domain receptor gene promoter.^14,15 Another study demonstrated that 100 μM eicosapentaenoic acid (EPA) directly modifies the methylation status of CCAAT / enhancer-binding protein delta gene.¹⁶ As for cellular disease models, palmitic acid (PA) was shown to induce global DNA hypermethylation in primary human myocytes and ex vivo human pancreatic islet cells at a 500 μM and 1 mM dose, respectively, affecting targets such as the PPAR-γ co-activator 1A gene.^17,18 Furthermore, a recent study in a cellular model of hepatic cancer has shown that a mixture of oleic acid (OA) and PA elicits hypermethylation of selected imprinted gene promoters.¹⁹ Possible mechanisms of epigenetic regulation by FAs include binding to PPARs, a family of transcription factors that regulate numerous metabolic processes via ligand-dependent transcriptional activation and repression.^20,21

Currently, it is unknown whether the above-described epigenetic effects are FA-specific, as is the contribution of FAs to disease-related methylation profiles. To understand those issues, we focused on the 2 long-chain unsaturated FAs, OA and AA, which are known to exert generally opposite cellular inflammatory responses.^22-24 We studied the effects of these particular FAs on the epigenome and transcriptome of THP-1 cells, a widely accepted human monocyte model,²⁵ and compared our results to available DNA methylation data of several human diseases and normal tissues. The implications of our results are discussed in the context of current knowledge of epigenetic regulation by lipid components and dietary factors, and their contribution to disease risk.

Results

Effects of the pure AA and OA on global DNA methylation in cultured cells

We first examined the effects of the pure FAs, AA, and OA, on global DNA methylation, i.e., total normalized 5mdC content in cultured human THP-1 monocytes. Stimulation experiments were carried out for 24 h using FAs in the 0–200 μM concentration range. These concentrations are below or within the reported circulating FA range—see, for example, Higashiyama et al.²⁶ The rationale for using a 24-hour stimulation is that epigenetic responses to lipoproteins were observed in THP-1 macrophages, a differentiated version of THP-1 monocytes.^12,13 In accordance with a similar study of AA-stimulated THP-1 cells, where cell proliferation was scored based on ³H-thymidine incorporation,²⁷ FAs did not affect cell proliferation as assessed by cell counting. Overall, AA and OA elicited distinct responses. AA induced a dose-dependent DNA hypermethylation peaking at the 100 μM dose and amounting to a ∼10.5% increase in 5mdC content at 100 μM relative to the 1 μM dose (Fig. 1A). In turn, OA induced a weaker response, with an overall DNA hypomethylation at 100 μM relative to the 1 μM dose. Noticeably, the effect of OA was significantly different from the one of the vehicle BSA only at doses >100 μM. Neither OA nor BSA elicited statistically significant responses relative to unstimulated cells or cells stimulated with the 1 μM of any FA, up to the 50 μM dose. To validate the divergent AA and OA dose responses, a 24-hour co-stimulation experiment, in which one FA was held constant at 100 μM concentration while the other varied between 1–100 μM, was performed. The results confirmed the distinct DNA methylation responses to AA and OA in THP-1 cells [Fig. 1B; note that the respective responses at 100 μM were not different (P = 0.08)]. Importantly, the observed AA- and OA-induced DNA methylation changes were not specific for THP-1 monocytes, as human embryonic kidney 293 cells also displayed distinct responses to these FAs following a 24-hour stimulation (Supplementary Fig. 1). These experiments were exhaustively repeated and their results were consistent across time (2004–2013), cell culture laboratories, and THP-1 cell stocks (Sweden, Mexico, and Spain for either), HPLC platforms (Mexico and Spain), and total 5mdC assays (HPLC-based or the ELISA-based MethylFlash system).

Figure 1. Effects of pure FAs on global DNA methylation in THP-1 monocytes. A, FA dose-response following a 24-hour stimulation. B, co-stimulation with AA and OA, in which each FA was held constant at the 100 μM dose (symbols in graph A) and the other varied between 1–100 μM (indicated as “variable FA” in the horizontal axis legend). Data in A were obtained by the HPLC-based 5mdC quantitation, data in B by the MethylFlash assay (see text for details). Data points represent averages and SD values. The horizontal dashed line indicates the value corresponding to unchallenged (control) cells. The latter value differs in panels A and B because of the different assays used. C, ALU methylation levels. Symbols are as in A. In all cases, asterisks indicate the significance of the difference in comparison with the respective 1 μM dose. Horizontal brackets in A indicate the AA or OA doses at which the response is significantly different from the response to the vehicle BSA. *, P<0 .05; **, P<0 0.01 (Scheffé post hoc test).

As an independent validation of the total 5mdC data, we performed an ALU-specific methylation assay. ALUs are an abundant transposon family and have often been used as a surrogate of total 5mdC determination.²⁸ In accordance with total 5mdC data, AA induced a significant dose-dependent increase in ALU methylation while OA showed a tendency to hypomethylation, although not statistically significant in the 1–100 μM range (Fig. 1C).

Effects of AA and OA on specific CpG site methylation

To detail genomic regions undergoing FA dose-dependent differential DNA methylation at the CpG level resolution, we analyzed THP-1 monocytes stimulated with 1, 10, or 100 μM AA or OA, using the Infinium HumanMethylation450 BeadChip arrays (450K array). The platform represents >485,000 intergenic and intragenic, non-repeated element CpGs.^29,30 The rationale for using 100 μM as maximal dose was that it yielded the most divergent effects between the 2 FAs in the total 5mdC assay.

Unsupervised clustering analysis of 450K array data showed that the FA-induced methylation profiles were more divergent at 100 μM FA dose compared to the relatively similar 1 and 10 μM-dose profiles (Fig. 2), mirroring global 5mdC data (compare with Fig. 1A). Next, we sought evidence for an AA dose-dependent DNA hypermethylation response and an opposite one for OA in 450K array data. We restricted our analysis to cells challenged with increasing—1, 10 and 100 μM—FA doses, as the unstimulated cell data reflect a non-physiological FA-free situation. A further reason for not including BSA-stimulated cells in the 450K array analysis, was that global 5mdC data indicated that the global trend of OA-induced DNA hypomethylation was statistically different from the corresponding BSA data set. Before conducting a detailed analysis of the 450K array data, we validated the methylation profiles of 11 CpGs corresponding to 6 genes by pyrosequencing. Genes were randomly selected among CpGs with absolute Δβ>0.1 between 1–100 μM AA. In all cases, pyrosequencing reproduced the DNA methylation profile of the 7 CpGs represented on the 450K array and additional flanking 4 CpGs (Supplementary Fig. 2). Subsequently, we sought autosomal CpGs that showed dose-dependent DNA methylation profiles, i.e., CpGs with statistically significant FA dose-dependent DNA methylation changes. We initially set a Δβ>0.2 between the extreme 1 and 100 μM FA doses as a stringent criterion to identify FA dose-dependent profiles. Such a filtering yielded 283 and 135 CpGs for AA and OA, respectively, a figure that is inconsistent with the global 5mdC data shown above. We therefore reasoned that FAs imposed small and widespread CpG methylation changes. Consequently, CpGs displaying an absolute Δβ>0.005 between consecutive FA doses (i.e., with absolute Δβ>0.01 between the extreme 1 and 100 μM FA doses) were considered for further analysis. We first separated the CpG sets for each FA into 2 subsets that displayed a positive or negative Δβ between extreme FA doses, and performed a Spearman's correlation test separately on the CpGs subsets. This allowed us to perform a one-tailed test using an absolute r>0.985 threshold value corresponding to P<0.05, for each subset. Those CpGs will be referred to as AACpGs or OACpGs. The analysis identified 57,187 AACpGs and 48,917 OACpGs, corresponding to 14,483 and 13,914 genes, respectively (a complete list of AACpGs and OACpGs can be provided upon request). The absolute Δβs between 100 and 1 μM FA doses ranged between 0.29–0.01 for the 2 sets compounded. None of the AACpGs or OACpGs reached the genome-wide significance threshold of P<10⁻⁷. A minority of AACpGs (6,432 or 11.2%) overlapped with OACpGs, corresponding to 3,722 genes.

Figure 2. Unsupervised clustering analysis of high-density DNA methylation array (450K array)-based analysis of FA dose-dependent CpG methylation.

Biological significance of AACpGs and OACpGs

Although the AA- and OA-induced differential methylation did not reach genome-wide significance, a substantial number of converging findings support a pathobiological significance for AACpGs and OACpGs.

First, in agreement with global DNA methylation data, a general tendency for AA-induced hypermethylation and OA-induced hypomethylation was observed, i.e., the majority of AACpGs (76.4%) and OACpGs (74.7%) showed dose-dependent hyper- and hypo-methylation, respectively (P = 6.7×10⁻³², Chi-square test; Fig. 3). Likewise, the compounded weight of hypermethylated AACpGs was higher than the hypomethylated counterparts: the sum of Δβs between 100 and 1 μM AA was 1,661.3 and −408.7 for the 2 sets, respectively. Conversely, the net effect of OACpGs was hypomethylation, since the sum of Δβs between 100 and 1 μM OA was 381.6 and −1,308.7 for the hyper- and hypo-methylated OACpG set, respectively. This overall divergent DNA methylation profile included the 6,432 common CpGs between the AACpG and OACpG sets, as the M values between 100 and 1 μM FA doses showed a negative correlation between the 2 sets (r = −0.14, P<10⁻⁶).

Figure 3. - Distribution of hypermethylated and hypomethylated AACpGs and OACpGs. The AACpG and OACpG sets are ordered for increased Δβ between the 2 extreme FA doses. The horizontal dashed lines indicate the overall average Δβ for each CpG set.

Second, at the gene level, most AACpGs and OACpGs co-localized to gene bodies and promoters and, at the same time, were remarkably divergent for the direction of DNA methylation change. For all overlapping genes, AACpG and OACpGs showed hyper- and hypo-methylation, respectively, between extreme FA doses (P<10⁻⁵ in all cases). This tendency was particularly evident of genes that harbored >50 AACpG or AOCpG/gene, which included protein tyrosine phosphatase, receptor type, N polypeptide 2 (PTPRN2), mitotic arrest deficient-like 1 (MAD1L1), PR domain containing 16 (PRDM16), tenascin-XB (TNXB), regulatory-associated protein of mTOR (RPTOR), inositol-1,4,5-trisphosphate 5-phosphatase (INPP5A), adenosine deaminase, RNA-specific, B2 (ADARB2), sidekick cell adhesion molecule 1 (SDK1), ATPase, class VI, type 11A (ATP11A), tubulin folding cofactor D (TBCD), and disco-interacting protein homolog 2 (DIP2C). Furthermore, such divergence in FA-specific methylation profiles was also observed of at least 2 large intergenic regions spanning 605.7 kb and 758.0 kb in size on chromosome 4 and 8, respectively. Fig. 4 shows illustrative examples of reciprocal AA- or OA-induced methylation profiles of three genes that showed more >50 AACpG or AOCpG/gene and of the intergenic region on chromosome 8. All, or only overlapping AACpGs and OACpGs, are shown (left and right panels, respectively).

Figure 4. Examples of locus-specific methylation profiles of AACpGs and OACpGs. Graphs on the left show the changes in DNA methylation of AACpGs (solid circles) and OACpGs (open squares). P indicates the position and number of CpGs mapping to promoters. Graphs on the right show the Δβ values of common CpGs to the AACpG and OACpG sets in each gene or intergenic region. P-values refer to unpaired (left) and paired (right) non-parametric tests (Mann-Whitney U test and Wilcoxon test, respectively).

Third, functional annotation analysis using the top 150 genes ranked by number of differentially methylated AACpG and OACpGs, normalized by the number of probes present in the 450K array for a given gene (N-CpGs), revealed significant enrichments for the G protein-coupled receptor (GPCR) and olfactory signaling pathways (REACT_14797 and hsa04740, respectively; Supplementary Table 1). Conversely, no signaling pathways were identified for the top 150 genes ranked by Δβ (DB-CpGs) (Supplementary Table 1). Remarkably, functional gene analysis of differentially expressed genes between 100 and 1 μM AA and OA doses revealed a significant enrichment of the same GPCR signaling pathway (FDR < 10⁻⁴ in both cases; see paragraph below “FA-induced CpG methylation and gene expression”). Notably, the average absolute Δβ value of the top 150 N-AACpG and N-OACpGs were ∼4-fold smaller relative to DB-AACpG and DB-OACpGs (0.042 and 0.044 vs. 0.174 and 0.161 on average, respectively; Supplementary Table 2), again pointing to small but biologically relevant effects.

Fourth, we compared AACpG or OACpG to a study addressing the effects of PA on the DNA methylome of primary human pancreatic islets, a model of pro-inflammatory free FA-induced diabetes.¹⁸ Notably, reminiscent of the effects of AA in THP-1 monocytes, PA induces a net DNA hypermethylation relative to pancreatic islets stimulated with FA-free BSA alone. In a first instance, we asked whether any differentially methylated CpGs were shared by the 2 studies, i.e., the dose-dependent AACpGs and OACpGs, and PA-induced differentially methylated CpGs (PACpGs). Remarkably, despite significant experimental design differences between the PA-stimulated pancreatic islet model and AA- or OA-stimulated THP-1 cells, including FA concentration (1 mM and 1–100 μM range), incubation time (48 and 24 hours) and cell type (pancreatic islets and THP-1 monocytes), a 10–12% overlap between CpGs was observed. A total of 445 CpGs corresponding to 430 genes, were common to the AACpG, OACpG and PACpG datasets (referred to as FACpGs; Supplementary Table 3). A functional annotation analysis revealed a significant enrichment for regulation by a number of transcription factors, including the biologically plausible CCAAT/enhancer binding protein (C/EBP), α (CEBPA), a critical activator of adipocyte-specific genes (FDR = 8.2×10⁻⁷) and PPAR-α (FDR = 0.03), but not the related PPAR-γ (FDR = 0.33).^31,32 Furthermore, we found that the M-values between AACpGs and PACpGs were positively correlated (r = 0.21, P = 1.20×10⁻⁵) and the majority (72.8%) were hypermethylated in both sets. Conversely, no correlation between the OACpGs and PACpGs was observed (r = 0.08, P = 0.10) and the majority (70.3%) of FACpGs were hypo- and hyper-methylated in the 2 sets, respectively. Furthermore, akin to AA and OACpGs, functional analysis of the top 150 normalized (number of differentially methylated CpG divided by the number of probes in the 450K array) PACpGs revealed enrichment in G-protein and olfactory signal transduction pathways (Supplementary Table 1). These data suggested a functional divergence between AA and PA on the one hand, and OA on the other hand with respect to DNA methylation.

The analogies between FA-stimulated THP-1 cells and a model of pancreatic inflammation raised the question whether similar coincidences exist within a range of human diseases. To address that issue, we compared AACpGs and OACpGs to the 450K array-based human DNA methylation profiling data sets that were available to date. These datasets covered a wide range of disease typology and conditions, i.e., metabolic diseases (atherosclerosis, obesity, type 2 diabetes),^33-35 psychiatric disorders (autism—Brodmann areas 10 and 24 profiling—and schizophrenia),^36,37 cancer (leukemia, colorectal carcinoma, hepatocellular carcinoma, blood from breast cancer patients),^38-41 and aging⁴² (Supplementary Table 4). The inclusion of the mentioned studies is justified by the general relevance of FAs as suppliers of energy and biologically active metabolites, and their recognized involvement in psychiatric disorders (reviewed in⁴³). In addition, the comparative analysis included 450K array-based DNA methylation profiles of 8 non-diseased tissues.⁴⁴ The portion of differentially methylated CpGs identified in each study, that overlapped AACpG or OACpGs ranged between 12–17%. Furthermore, despite differences in Δβ cut-off between arrays (Supplementary Table 4) we observed no significant correlation between the extent of overlap and Δβ cut-off values. Excluding the breast cancer blood data set, which did not represent the disease target tissue, we observed a notably lower overlap in cancer compared to metabolic, psychiatric and normal tissue sets (P<0 0.01 in all comparisons, Scheffé post-hoc test; aging was not included in the analysis given that it represented a single study). Since we found no correlation between extent of overlap and number of differentially methylated CpGs yielded by each study (r = −0.12, n = 20, P = 0.62), our data suggested a non-random overlap despite the relatively high number of AACpGs and OACpGs. Gene function analysis for the top 150 N-CpG-harboring genes revealed a common significant enrichment in the GPCR signaling pathway among the cancer and autism (BA24) profiles (FDR<10⁻² in all cases), but not in normal tissue samples, suggesting an underlying biological specificity (in italics in Supplementary Table 1).

Next, we compared the methylation state of AACpGs and OACpGs to that of differentially methylated CpGs identified in diseased or normal tissue. Clustering analysis of CpG methylation levels (Δβ) averaged by gene, clearly showed that the AA- and PA-induced profiles were more similar to each other than to OA-induced profiles. Furthermore, the former clustered with metabolic disease, aging and autism profiles, whereas cancer, schizophrenia and normal tissue-specific profiles grouped away from FAs (Fig. 5A). These patterns were readily observed in the clustering analysis of the top 50 genes, ranked by the number of overlapping CpGs (Fig. 5B). For example, the majority of genes were consistently hypermethylated following stimulation with AA or PA but hypomethylated following stimulation with OA.

Figure 5. Clustering analysis of CpGs that are common to AACpG or OACpG or both, and to the indicated disease and normal tissue DNA methylation profiles. A, Overall data clustering. B, Top 50 genes ranked by the number of harbored CpGs included in the clustering analysis. Asterisks indicate the top 11 genes ranked by the number of AACpG/gene. Red and blue, hypermethylation and hypomethylation, respectively. The tissue origin or relation is indicated for non-normal tissue samples as follows: Ad, adipose tissue; Ao, aorta; Bl, blood; Br, brain; In, intestine; Li, liver; Pa, pancreas.

Fifth, we asked whether AACpGs significantly overlap with a set of 450K array methylation profiles obtained in human aortic atherosclerotic lesions of varying histological severity.⁴⁵ That study identified 1,985 autosomal CpGs, which significantly change their methylation status with lesion progression (grade-CpGs), mostly by drifting toward hypermethylation as lesion grade increases. Grade-CpGs were significantly enriched in AACpGs, as 510 CpGs overlapped between the 2 sets and 324 showed the same methylation trend, mostly (311 out of 324) toward hypermethylation both in response to increasing AA doses and with lesion progression (P = 0.011 for enrichment, hypergeometric test). A list of the 324 overlapping and concordant CpGs between the 2 sets is shown in Supplementary Table 5.

Sixth, AA- and OA-induced DNA methylation profiles showed a degree of specificity for sequences with specific epigenetic and functional signatures. As for preference for hypo- or hyper-methylated sequences, we analyzed the distribution of the AACpG and OACpG across β-values. We observed a clear enrichment of the initially hypermethylated group, i.e., at 1 μM FA dose for both AACpGs and OACpGs, that was independent of the direction of methylation change (Supplementary Figure 3). As for gene compartments and GC content, hypermethylated AACpGs preferentially mapped to sites within open sea (i.e., >4 kb from the nearest CpG island), whereas hypomethylated AACpGs mapped preferentially to CpG islands (P = 0.031 and P = 0.009, respectively, compared to the 450K array probe distribution; Supplementary Figure 4). No other CpG set showed any significant mapping preference.

Overall, these findings support a broad biological significance and instructive mechanism(s), rather than random seeding, for the overall comparatively small AA- and OA-induced DNA methylation profiles.

FA-induced CpG methylation and gene expression

To assess the impact of dose-dependent changes in DNA methylation on gene expression, RNA was extracted from the same cells used to produce 450K array data and analyzed using the Affymetrix Human Exon 1.0 ST platform. To probe for an association between expression and DNA methylation, we focused on AACpGs and OACPGs that colocalized to either exonic or promoter regions. For promoter and exonic CpGs, we correlated dose-dependent changes in DNA methylation with average expression values of all exons or the corresponding exon probe, respectively. While DNA methylation changes in any genomic context showed both positive and negative associations with transcription, context-dependent correlations were uncovered; in particular, a negative correlation in the promoter, 5′UTR and first exon and a positive correlation within gene body and 3′UTR regions. Importantly, average exonic Δβ values between the 100 and 1 μM FA dose were higher of AACpGs compared to OACpGs (0.023 and −0.021, respectively). Since this mirrored the overall tendency of all AACpG and OACpGs, we concluded that gene body methylation was more frequently associated with an increase, rather than a decrease, in transcription. RT-PCR analysis of genes that showed a high number of opposite dose-dependent changes in DNA methylation following stimulation with AA and OA, also supported a tendency for a positive association between gene body methylation and expression for some of the genes analyzed (Fig. 6B,C). Furthermore, functional gene analysis revealed a significant enrichment in GPCR signaling (FDR<10⁻⁴ in both cases), mirroring the functional enrichment observed in FA-induced methylation profiles (see above).

Figure 6. Effects of dose-dependent changes in DNA methylation on gene expression. A, Correlation between 100–1 μM AA or OA Δβ and expression array-based corresponding expression levels. B, RT-PCR analysis of selected genes following AA and AO stimulation of THP-1 cells for 24 h at 1 and 100 μM concentration. C, GAPDH RNA-normalized expression. The bars for each gene indicate, left to right: 1 μM AA, 100 μM AA, 1 μM OA, 100 μM OA, respectively.

Factors mediating FA-induced changes in DNA methylation

To explore potential mechanisms underlying AA- and OA-induced changes in DNA methylation, we focused on PPAR-α and PPAR-γ, given that both are expressed in THP-1 cells and that FAs, or derivatives thereof, are known ligands of these nuclear receptors.^46,31 First, because the above presented functional annotation data suggested a potential role for PPAR-α, but not PPAR-γ in regulating FACpG methylation and, second, due to the link between G-protein cannabinoid receptor signaling pathway and PPAR-α transcription pathway,^47,48 we examined the effects of the PPAR-α and PPAR-γ antagonists GW6471 and GW9662, respectively, on the FA-induced methylation responses. GW6471, significantly and substantially, inhibited global DNA hypermethylation induced by 100 μM  ΑΑ ιν a doσϵ-dependent fashion, while no effect was observed in cells stimulated with OA (Fig. 7A, B). As for GW9662, neither AA- nor OA-induced methylation responses were significantly affected. Carnitine palmitoyltransferase-1 (CPT1), an enzyme that regulates the import of long-chain FA transport into the mitochondria, is a known target of PPAR-α. Akin to the PPAR-α antagonist, etomoxir, an inhibitor of CPT1, significantly inhibited global DNA hypermethylation induced by 100 μM  ΑΑ ιν a doσϵ-dependent fashion (Fig. 7A), while no significant effect of cells stimulated with OA was observed (Fig. 7B).

Figure 7. Participation of PPARs and β-oxidation in AA- or OA-induced global DNA methylation changes. THP-1 monocytes were stimulated with 100 μM FA for 24 h following a 1 h pre-stimulation with inhibitors of β-oxidation (etomoxir), PPAR-α (GW6471) or PPAR-γ (GW9662). Data are represented as difference in DNA methylation between AA and AO-stimulated and unstimulated cells. The corresponding inhibitor concentration is indicated (μM). Statistical significance is in comparisons with cells grown in the absence of any inhibitor (C). The horizontal dashed line indicates the baseline (unstimulated cells) level. For statistics symbols and test, see legend of Fig. 1.

Numerous data lend support to complex interactions between PPARs and sirtuins, the latter of which also regulate transcriptional networks of critical metabolic processes via their NAD+ dependent histone deacetylase activity.⁴⁹ To understand whether sirtuins played a role in the FA methylation responses, mouse embryonic fibroblasts (MEFs) obtained from SIRT1-, SIRT2-, or SIRT6-null mice and their respective matched WT controls, were stimulated with 100 μM AA or OA for 24 hours. These particular sirtuins were chosen due to their facultative or preferential nuclear location and prior evidence for their regulatory roles in FA oxidation (reviewed in⁴⁹). Importantly, both AA and OA elicited methylation responses in WT MEFs, which were comparable to those observed in human THP-1 cells. However, among the analyzed mutant MEFs, AA-, and OA-specific DNA methylation responses were only significantly inhibited in SIRT1-null cells (Fig. 8A). As further validation of SIRT1 participation in AA-induced hypermethylation, we examined the effects of SIRT1 inhibitors, sirtinol and splitomycin, in THP-1 monocytes. Both inhibitors led to reduced levels of DNA methylation in cells stimulated with 100 μM AA, while TSA, a general histone deacetylase inhibitor with a preference for non-NAD+ dependent histone deacetylases, failed to elicit any effects (Fig. 8B).

Figure 8. Participation of histone deacetylases in AA- or OA-induced global DNA methylation changes. A, DNA methylation response in MEFs genetically null for SIRT1, −2, or −6 (indicated as “null”) and matched WT controls, stimulated with 100 μM AA or OA for 24 h, or unstimulated (U). Statistical significance is in comparison with corresponding WT cells. B, Effects of inhibitors of SIRT1 (sirtinol, splitomycin) and broad range histone deacetylases (TSA) on the 100 μM AA-induced response. For statistics symbols and test, see legend of Fig. 1.

Discussion

We have uncovered distinct DNA methylation profiles induced by AA and OA, 2 unsaturated FAs generally known for their pro-inflammatory and anti-inflammatory properties, respectively. Our global, ALU-specific and individual CpG methylation data obtained in various laboratories across a long period of time (see Results), consistently show divergent effects of the 2 FAs in the analyzed dose range, i.e., AA and OA induce DNA hypermethylation and hypomethylation, respectively, both in human and mouse cultured cells. A comparative analysis of publicly available human DNA methylome data, further confirm the biological relevance of these findings, despite the different criteria used to identify relevant CpG loci. The pro-inflammatory FAs AA and PA induce markedly similar DNA methylation profiles, compared to the corresponding OA-induced profiles. Yet, FA-induced methylation and expression profiles are both enriched in the GPCR signaling function, suggesting antagonistic effects on common cellular pathways. Those similarities are further supported when comparing FA-induced profiles with disease- and normal tissue-specific corresponding data. AA and PA cluster together with metabolic diseases (atherosclerosis, obesity, type 2 diabetes), aging, and autism. OA, on the other hand, clusters more weakly with the latter conditions. This clustering therefore closely reflects the known distinct pathobiological properties of AA and PA on the one side, and OA in the other. Furthermore, several of the identified target genes are relevant from a nutritional and evolutionary perspective. For example, a 12-week supplementation with AA and DHA in baboons is associated with downregulation of PTPRN2 in the cerebral cortex.⁵⁰ Likewise, an association between DNA methylation, a sequence variant of PTPRN2 and islet insulin secretion has also recently been identified.⁵¹ In turn, TNXB is one of the two genes showing population-specific polymorphisms associated with dietary adaption to fat, milk, and meat.⁵² These associations, although independent and converging, are not definite proof of causality. Also, the effects of FAs presented here are quantitatively modest, reflecting the 10% or less change in DNA methylation profiles that was detected in response to FA, high-fat diet and exercise.^18,53,54 Despite these limitations, the data provide testable hypotheses to better understand the regulation of the DNA methylome by FAs and possibly to improve diet-based strategies to improve human health.

These comparisons offer a number of insights into the pathophysiological relevance of the epigenetic responses to FAs. In the case of human atherosclerosis, we show concordance with 2 previously published studies. First, the histological grade-associated CpGs⁴⁵ are significantly enriched in AACpGs. Second, 16 of the most represented genes across the analyzed DNA methylation profiling studies, are densely hypermethylated in the early stages of the aortic atherosclerosis and display a direction of DNA methylation change that is concordant with the one of AA and PA.³³ In both atherosclerosis and AA-stimulated cells, hypermethylation is prominent within the body of these genes. Taken together, these comparative data suggest that a portion of the DNA hypermethylation that characterizes atherosclerosis in its initial phases and during the progression of the stable vascular lesion,^33,45 may be seeded by circulating pro-inflammatory FAs or their metabolites. Coupled with evidence that hypomethylating agents can slow the progression of the vascular lesion,⁵⁵ our observations may lead to testable nutritional strategies to reprogram or safeguard the vascular epigenome. As for autism, the involvement of AA in that condition has been established by early prospective studies documenting the effects of AA-depleted infant formulas on autism risk and by a recent comparative study of red blood cell AA content between autistic patients and controls.^56,57 Another finding of our study is the divergence of type 2 diabetes DNA methylation profiles from the ones of AA and, notably, PA. The hypermethylated loci identified in both AA- and PA-stimulated cells show a distant, if any, resemblance to the corresponding profiles in diabetes, thus revealing a set of CpGs that do not model the involvement of free pro-inflammatory FAs in diabetes. Additionally, we show that FA-induced profiles are relatively distant from normal tissue-specific ones. This result is expected, as due to the ubiquitous access of circulating FAs to organs and tissues, it is likely that FAs seed DNA methylation profiles that represent common features of the majority of the organism´s epigenomes. Incidentally, our findings confirm and expand the previous observation that normal blood displays a strikingly unique DNA methylation profile, compared to the other tissues and diseases tested.⁴⁴ In addition to all these considerations, it is noteworthy that the specificity of the overlap with a number of published studies constitutes a powerful independent validation of our 450K array-based results.

Pathways linked to β-oxidation play a pivotal role in determining FA-related methylation profiles. We show that AA-induced DNA hypermethylation is sensitive to inhibition of FA import to mitochondria and is downstream to PPAR-α and SIRT1, 2 central regulators of FA metabolism and β-oxidation.⁴⁸ THP-1 monocytes, in contrast to cancer cells in general, rely more on β-oxidation than glycolysis for energy production and AA has been shown to be a substrate for β-oxidation in THP-1 cells.^58,59 Both these observations support the effects of β-oxidation inhibition in our model system. Indeed, the comparatively small overlap between cancer and FA-induced DNA methylation profiles may be explained by the Warburg effect which states that cancer cells mainly produce ATP from glycolysis even in conditions of hypoxia.⁶⁰ Following this logic, the DNA hypomethylation observed in diabetes is a predictable consequence of the defective mitochondrial activity associated with that disease.⁶¹ In general, blockage of the mitochondrial machinery rapidly shifts ATP production to glycolysis to assure cell survival. We assume that such an adaptation takes place in THP-1 cells following exposure to β-oxidation inhibitor etomoxir, suggesting that inhibition of AA-induced DNA hypermethylation by etomoxir does not result from reduced ATP supply to chromatin modifiers. On the other hand, the OA-induced hypomethylation response is not affected by these inhibitors. This suggests that the distinct AA and OA-induced methylation profiles might relate to FA-specific differences in β-oxidation. Indeed, both in vitro and a recent in-vivo dietary study in humans show that OA relative to other FAs, such as PA and EPA, is a poor substrate of β-oxidation,^62-64 although the opposite has also been reported.⁶⁵ Alternatively, the increased levels of specific intermediate β-oxidation products observed in PA- relative to OA-stimulated human skeletal muscle cells could also explain part, or all of, the PA-induced methylation.⁶² Although our data contrast with the predominant view that SIRT1 is activated solely in condition of starvation, NAD+ independent SIRT1-mediated upregulation of oxidative metabolism has been reported in several other studies.^66-68 Finally, the observation that SIRT1 participates in the AA-induced DNA hypermethylation response in MEFs, indicates that the corresponding response observed in THP-1 cells is not a cancerous cell-restricted, p53- and MIR34a-dependent event.⁶⁹

The PPAR-α specific regulation of AA-induced DNA methylation, presumably reflects the observation that polyunsaturated FAs are more potent activators of PPAR-α than PPAR-γ.⁷⁰ The involvement of PPAR-α, but not PPAR-γ, is further supported by the results of the functional annotation analysis of FACpG-harboring genes. SIRT6 can also be activated by long-chain FAs, in particular OA, and contributes to regulation of β-oxidation.^71,72 However, knockdown of SIRT6 in MEFs had no effect on FA-induced changes in DNA methylation.

Our observation that AA and OA participate in shaping metabolic disease-specific DNA methylomes through β-oxidation, PPAR-α, and sirtuin 1 signaling, has potential implications for diet-oriented therapy and prevention. Indeed, it has been argued that the high-fat, low-carbohydrate ketogenic diet, which is successful in clinical management of patients with inborn errors of metabolism and pharmacoresistant epilepsy, may be beneficial in a much broader range of diseases than previously recognized.^73,74 Conversely, a low-fat, low-carbohydrate diet can lead to normalization of β cell function and insulin resistance in individuals with type 2 diabetes.⁷⁵ However, whether or not a shift in substrate preference toward fat oxidation lowers disease risk is still heavily debated. For example, prolonged dietary administration of the PPAR-α agonist WY-14643, is associated with a gradual decreased in global and LINE DNA methylation and hepatocarcinogenesis in mice.^76,77 Furthermore, the differential tumorigenic responses of xenotransplanted leukemia cells THP-1 or NB4 exposed to either a high-fat or high-carbohydrate diet,⁷⁸ underlines the importance of understanding energy metabolism of disease subtypes prior to dietary interventions.

Furthermore, we uncover a positive correlation between gene body methylation and expression in THP-1 cells. That correlation represents the behavior of a majority, but not all genes, underlining a complex association between transcription and gene body methylation. Similar conclusions have been drawn in several genome-wide studies of DNA methylation with recent reports suggesting that gene body methylation increases gene expression levels.^79,80

Methods

Cell culture

THP-1 monocytes were cultured in RPMI-1640 medium (Gibco) supplemented with 2mM L-Glutamine (Sigma), 10% fetal calf serum (Gibco) and 1% Penicillin/Streptomycin (Gibco). Cells were never allowed to grow above a 1–1.5×10⁶ cells/ml concentration. MEFs obtained from SIRT1- or SIRT6-null mice were a kind gift from Dr. Raúl Mostoslavsky (Harvard University, Boston, USA) and Dr. Eva Bober (MPI, Bad Nauheim, Germany), respectively. SIRT2-null MEFs were generated by standard procedures from SIRT2-null mice. WT MEFs were derived from littermates of each mutant mice strain. MEFs and HEK293T cells were grown in the same conditions as THP-1 monocytes but with DMEM as medium. Pure FAs (Sigma) were conjugated with cell culture-grade FA-free BSA (fraction V, FA-free, Sigma no. 820022) to achieve a FA:BSA 6:1 ratio, essentially as described.⁸¹ Typically, 5–6×10⁶ cells in 10 ml medium were stimulated with 100x BSA-FA mix in 2% FCS. Exclusion of trypan blue was used as a criterion for viability. The inhibitors etomoxir, GW9662, GW6471, sirtinol, splitomycin or trichostatin A (TSA) were used at the concentrations and experimental conditions reported in previous studies.^82,83 At least 3 technical and biological replicates were performed for each experiment, excluding microarray experiments.

Global DNA methylation

DNA was extracted from cultured cells by standard methods (DNeasy system, Qiagen). For the measurement of global DNA methylation levels, 5-methyl-2′-deoxycytidine (5mdC) and 2′-deoxyguanosine (dG) were determined by an HPLC-based method.⁸⁴ Normalized 5mdC levels were calculated as percentage 5mdC/dG. When indicated, global DNA methylation was measured by the MethylFlash system (Epigentek). The assay yielded 5mdC values as percentage of PicoGreen-quantified input DNA (wt/wt). Each value was multiplied by 4 to obtain a final value comparable to that obtained by the HPLC-based method. The MethylFlash assay yielded responses to AA and OA in THP-1 monocytes, that were consistent with the ones obtained by the HPLC-based 5mdC quantitation, except for baseline (unstimulated cells) values that were 5% lower (P>0 0.05, data corresponding to a total of 15 and 21 experiments, respectively).

ALU repeat methylation

ALU element methylation was determined according to Sirivanichsuntorn and collaborators.²⁸ Briefly, bisulfite-treated DNA was amplified by PCR with the ALU-specific primers AluF: 5′-GGY GYG GTG GTT TAY GTT TGT AA-3′ and AluR: 5′-TTA ATA AAA ACR AAA TTT CAC CAT ATT AAC CAA AC-3′ as follows: 95°C for 30s, 53°C for 20s (40 cycles), and 72°C for 20s. Subsequently, the 117 bp ALU amplicon was digested overnight at 65°C with TaqI (Invitrogen). The resulting products, ranging between 42 and 117 bp depending on the methylation states of ALU elements, were separated by gel electrophoresis, stained with GelRed and quantified using ImageLab (BioRad).

DNA methylation arrays

To identify targets of FA-induced DNA methylation we used the Infinium HumanMethylation450 BeadChip platform (Illumina). Three biological replicates were pooled in equal proportion for each FA concentration. Pooling effectively averages the data and reduces noise.⁸⁵ DNA quality checks, bisulfite modification, hybridization, data normalization with the GenomeStudio software (Illumina) and Beta value calculation were carried out as described elsewhere.^29,86 SNP-proximal probes and probes with a detection P-value >0.05, i.e., not statistically different from the background, were discarded as reported.⁴⁵ The methylation level for each cytosine, expressed as a β-value, was calculated as the fluorescence intensity ratio of methylated to unmethylated versions of each probe. Beta values ranged between 0 (unmethylated) and 1 (methylated). Logit-transformed β-values (M-values) were used in statistical tests unless indicated.⁸⁷ β or Δβ values were used for data description, due to their more intuitive nature. The annotation relating to CpG islands (CGIs) uses the following nomenclature: “shore,” each of the 2 kb-sequences flanking a CGI; “shelf,” each of the 2 kb-sequences next to a shore; “open sea,” DNA not included in shores, shelves or CGIs.²⁹ TSS200 or TSS1500 indicate the region between position −200 bp or −1,500 bp and the Transcription Start Site (TSS), respectively. When a cytosine mapped to different genic elements due to the presence of multiple alternative transcripts, all mapping variants were counted. Functional gene annotation was analyzed with the DAVID tool (http://david.abcc.ncifcrf.gov/). 450K array data are accessible in NCBI's Gene Expression Omnibus⁸⁸ with the accession number GSE67331. DNA methylation array data of selected CpGs were validated by pyrosequencing (PyroMark Q96 ID, Qiagen). Sequences of pyrosequencing primers are shown in Supplementary Table 6.

Expression arrays

Genome expression was analyzed by Human Exon 1.0 ST arrays (Affymetrix) and we used the Affymetrix Expression Console™ software (version 1.3) and Robust Multichip Average algorithm (RMA) to normalize and analyze array data. The library used was the HuEx-1_0-st-v2.r2.pgf with the probe set HuEx-1_0-st-v2.r2.dt1.hg18.extended.mps (gene level) and HuEx-1_0-st-v2.r2.dt1.hg18.full.mps (exon level). Expression array data are accessible in NCBI's Gene Expression Omnibus⁸⁸ with the accession number GSE57076. Validation of selected targets was performed by semi-quantitative PCR normalized to GAPDH expression levels.

Statistics

The Kruskal-Wallis test was performed to test for differences in DNA methylation levels between treatments or doses for continuous variables. If the overall comparison was significant, samples were re-analyzed by ANOVA followed by Scheffé's post-hoc test to identify specific different groups. When comparing 2 groups, the Mann-Whitney U test (unpaired samples) or the Wilcoxon test (paired samples). The Chi-square test was applied to compare percentages. Correlations were tested by calculating the Pearson's correlation coefficient and the associated p value. Tests were performed with the Statistica software (StatSoft) or Excel statistics tools (Microsoft Office 2011 for Macintosh). The significance of CpG enrichment was estimated with the hypergeometric test using an online tool (www.geneprof.org/GeneProf/tools/hypergeometric.jsp).

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgments

We thank Drs. Holger Heyn and Anna Marazuela-Duque for advice and assistance; Dr. Luis Delaye and Alejandro Gonzalez for granting us access to their computer cluster facility (funded by the State of Guanajuato's Council for Science and Technology (CONCyTEG), project “Laboratorio de aprendizaje e investigación en cómputo biológico LAICBIO” no. GTO-2012-C02-187442). This work was generously funded by: CONCyTEG grant no. 08-03-K662-020-A01, Mexican Council for Science and Technology (CONACyT) Basic Science (Ciencia Básica) grant no. 83401 to G.L; CONACyT Sabbatical Fellowships no. 166209 to G.L. and no. 166058 to S.Z; CONACyT Ph.D. Fellowship no. 334370 to G.A.S.-M; the Spanish Ministry of Economy and Competitiveness (MINECO) grant SAF2011-25619 to A.V.