Epigenome-Wide Association Study of Mediterranean Diet on Blood DNA Methylation Patterns during Pregnancy: A Pilot Study

1. Introduction

The Mediterranean Diet (MD), characterized by a high intake of fruits, vegetables, whole grains, legumes, nuts, olive oil, and moderate consumption of fish and dairy products, has received attention for its health-promoting effects. Epidemiological studies consistently demonstrate that adherence to the MD is associated with a significantly reduced incidence of chronic non-communicable diseases (NCDs) such as cardiovascular disease, type 2 diabetes, and various cancers, and a lower all-cause mortality [1,2,3,4,5]. These protective effects are largely attributed to the MD diet’s anti-inflammatory properties, the ability to promote favourable gut microbiota, and the potential to modulate metabolic pathways [6]. Promoting adherence to a MD during pregnancy may be a promising prenatal strategy to support healthy pregnancies and for NCD risk reduction in both mother and the offspring. Inflammation in pregnancy can predispose to an elevated risk for pregnancy complications and NCDs in later life, and therefore reducing inflammation through dietary strategies could have beneficial effects on both mother and offspring. While specific mechanisms are still unclear, epigenetic changes may mediate MD diet effects on the host immune system.

Previous studies have shown that dietary metabolites typically associated with a MD can induce epigenetic changes in immune cells in proportion to the intake of various bioactive nutrients [7]. Flavonols, a class of phytochemicals, are found in high amounts in onions, broccoli and tomatoes and other fruits and vegetables commonly consumed in the MD and have been shown to influence miRNA expression patterns related to inflammation [8,9]. Sulforaphane (organosulfur compound) found in cruciferous vegetables, decreases DNA methylation levels in nuclear factor erythroid 2-related factor 2 (NRF2) promoter and has been shown to increase the antioxidant and anti-inflammatory capacity in neurodegenerative diseases [10,11]. Lignans, found in high quantities in extra-virgin olive oil, have been shown to influence epigenetic regulation, including DNA demethylation, and may indirectly affect histone modifications such as H3K9me3, which are associated with heterochromatin formation and chemo-preventive effects [12]. Adherence to MD in adults participating in the PREDIMED-Navarra randomized controlled trial was associated with differential methylation of CpG base modifications in COL18A1 gene encoding for the alpha chain of type XVIII collagen, and PPARGC1B transcription factor stimulant implicated in type 2 diabetes [13]. These diet-induced epigenetic modifications can impact the function and responsiveness of immune cells at the nuclear level, thereby influencing systemic inflammation.

BioMood is an observational study exploring how consumption of a MD during pregnancy influences the gut microbiome, metabolome and inflammation. Previous work by our group on this cohort has demonstrated that adherence to a MD during pregnancy is characterized by shifts in amino acids, increased concentration of gut-microbial metabolites and is inversely associated with biomarkers of systemic inflammation, GlycA and GlycB [14]. This study also demonstrated distinct differences in the inflammatory marker Supramolecular phosphocholine composite (SPC) , derived from phospholipid contained within the cardioprotective HDL-4 subfraction, which were modified in the HMDA group compared to the LMDA group [14].

The present study builds upon this work to explore whether genome-wide blood DNA methylation changes were detectable in women with high- versus low- adherence to MD throughout pregnancy. We conducted both hypothesis-free and hypothesis-driven analyses, examining whether serum metabolite changes observed in this cohort were associated with epigenetic changes in circulating immune cells. Our findings provide preliminary evidence for blood epigenetic changes that warrant further studies.

2. Results

2.1. Description of the Cohort

The demographic characteristics of those included in the study are depicted in Table 1. The ethnic profile of the cohort was mixed consisting of mostly Caucasian and South Asian ancestries. All women were taking common pregnancy supplements for iron, folate and vitamin D. Inflammatory morbidities (allergies, asthma, inflammatory bowel diseases) were generally low (less than 10% of the population) and equally distributed between groups.

2.1. Sequencing Quality Control and Exploratory Analysis

Whole blood DNA methylation sequencing depth averaged +196 million reads per sample. An average of 89.49% base calls per sample were of a quality score greater than or equal to Q30 (99.9% base call accuracy). During quality control analysis, four samples within the same target-capture pool were identified as poor quality based on visual inspection of Beta density plots and multidimensional scaling (MDS) analysis (Figure 1). These were removed from the data set and we retained 48 samples for analysis, with 7.7 million CpG measures. We derived principal components to understand major sources of variation in this data set and correlated the top five components with clinical measures of subject age, BMI and MDA group. Participant age was significantly correlated with component 2 (r=0.34, P = 0.02) and MDA group had weak evidence of anti-correlation with component 5 (r= -0.25, p= 0.08). Regression models were thus adjusted for the first four components as covariates.

2.2. Epigenome-Wide Association Analysis of MD Diet Group

We applied a genome-wide logistic regression analysis comparing MDA groups to initially search for associations across all chromosomes. No individual CpGs reached the threshold for genome-wide significance (p-value < 5e-8), however, 54 CpG sites met the threshold for a suggestive association (p-value < 1e-5), albeit the effect sizes were very small (Figure 2,Table 2).

2.3. Gene-Centric Analysis

To restrict the search space we also conducted a gene-centric analysis focusing on CpGs within gene promoter regions, as these are critical for gene regulation. Using the same regression model, we detected 16 CpGs located within promoter regions with suggestive associations with MDA (p value < 1e-5). A possible differentially methylated region (DMR), comprising 7 CpG units within a 44-bp window within the IFI30 gene encoded on Chromosome 19 exhibited a consistent gain of methylation in the HMDA group (Figure 3).

2.4. Candidate Gene Analysis of Previously Identified MD Genes

We also conducted a candidate gene analysis focusing on CpGs mapping to genes previously identified by Arpón et al. [13] as associated with MD adherance. We found moderate evidence (p value <= 1e-5 ) for differential methylation at two genes, COL18A1 and PPARGC1B however effect sizes were small (Figure 4, Table 3).

2.5. Association Between Inflammatory Markers and Methylation Levels

2.5.1. Glycoprotein Analysis

To supplement the binary logistic regression analysis of diet groups, we also conducted a more hypothesis-driven approach, employing linear regression to examine whether the circulating inflammatory markers - GlycA, GlycB, SPC₁, SPC₂ and SPC₃ - that were previously shown to vary between the LMDA and HMDA groups [14], were associated with any changes in site-specific DNA methylation. To address potential co-linearity, we performed a correlation analysis of GlycA and GlycB which highly positively correlated (r = 0.9549, p-value < 0.01) and therefore we used GlycA as the independent variable. Genome-wide linear regression models revealed 57 CpG sites with suggestive associations (p-value <1e-5), 13 CpGs of which were located in promoter regions, but none surviving stringent multiple testing adjustment. A potential DMR with suggestive associations was identified within an intergenic region of chromosome 13 comprising four CpGs located within 37 bases of each other (Figure 5, Table 4).

2.5.2. SPC Lipoprotein Analysis

A Pearson correlation analysis between SPC identified that SPC₂ and SPC₃ levels are highly positively correlated (r = 0.8459, p-value < 0.01) and subsequent analysis focused on SPC₃ only. Genome-wide linear regression detected one significantly different methylated CpG site within gene TSC2 with strong evidence for association with SPC₃ circulating levels (p-value < 5e-8, Figure 6). A further 87 CpGs sites exhibited suggested associations with changes in SPC₃ levels (Supplementary Table 2). A possible DMRs within ZFPM1 and an intergenic region of chromosome 6 were also identified (Table 5). We also repeated this analysis on SPC₁ levels, which were independent of SPC₂ and 3. No genome-wide significant associations were identified and 61 CpG modifications exhibited suggestive associations (Table 6).

3. Discussion

Summary of key findings.

This pilot study aimed to explore the relationships between Mediterranean diet adherence in pregnancy whole blood DNA methylation using a multi-faceted approach. The study was designed to generate hypotheses and assess feasibility as a prelude to larger-scale dietary methylome-wide association testing. Despite the limitations of a pilot study, robust signals were identified in specific areas with a priori biological plausibility, while broader genome-wide suggestive associations provide valuable hypotheses for larger-scale studies.

Genome-scale Analysis Findings.

Using stringent thresholds for genome-wide significant we did not identify any site-specific or regional associations with Mediterranean diet group. Despite this, 54 CpG sites met the criteria for suggestive associations, exhibiting small effect sizes. Collectively this indicates a lack of statistical power to detect small effects, as commonly observed in complex disease studies, using mixed tissues such as blood, and suggests true biological effects of diet on methylation are subtle. Future studies will require hundreds of samples to be sufficiently powered for robust associations at genome-scale. Narrowing the breadth of our search to include only gene-promoter regions identified suggestive regional associations on 7 consecutive CpGs, a possible differentially methylated region. These potential signals warrant further investigation, particularly since the IFI30 gene plays a vital role in immune function and is involved in immune-mediated diseases. Differentially methylated regions (DMR) are often more biologically relevant since the represent coordinated methylation changes across a genomic locus, potentially impacting gene regulation. Our results position this DMR as a strong candidate for subsequent studies, which should be conducted in purified relevant immune cell types to substantially improve statistical power.

Candidate gene analysis.

Similar to previously published findings [13], we observed moderate evidence for differential methylation within the inflammation-related genes COL18A1 and PPARGC1B. Again, effects sizes were small and not entirely consistent in direction (loss/gain) methylation, likely due to the mixed cell composition of the specimen. These genes play a role in adipogenesis, metabolism, and diabetes and perform functions associated with inflammation, which can be related to the described favourable effects of the Mediterranean Diet on the inflammatory response [15]. Collectively, these data generally align with pre-existing studies and provide stronger evidence for a true biological link.

We also investigated associations with circulating metabolite levels that were previously found to be modified by HMDA in this cohort, and we detected small effects on site-specific methylation [14]. Association testing with Glycoprotein levels revealed a potential DMR with chromosome 13, a non-coding region of the genome. Non-coding methylation changes can influence chromatin structure and potentiate gene expression through chromatin interactions. SPC₃ levels were associated with significant genome-wide methylation changes within the gene TSC2, a tumour suppressor gene. Although the biological consequences of these epigenetic associations are unclear, they are strong candidates for further investigation and could aid in understanding the link between diet, metabolism and immune health.

Strengths and Limitations.

A novel aspect of this study was the methylation sequencing platform, which compared to microarrays that have typically be used in this field of study [13], offers inherently more quantitative count-level data, is less sparse and covers nearly five times more of the mammalian genome than current generation microarrays. A key strength of the study is the robustly selected cohort of women studied with prior evidence of metabolic changes associated with Mediterranean diet adherence scores. Despite this, the primary limitation of the study is the small sample size and mixed cell tissue type, which increases the change of false negative associations and small effect sizes. While our extremes of exposure design enriches for biological signals, the generalizability of the findings may be limited to the specific characteristics of the pilot cohort and replication in larger more diverse populations is needed. We were unable to disentangle effects of highly correlated phosphopipids or glycoproteins so further in-depth analysis is required. The cross-sectional design precludes causal inference, and we acknowledge the potential influence of unmeasured confounding factors.

Conclusion

This study provides valuable hypotheses for future exploration of epigenetic mechanisms associated with diet in pregnancy and suggests potential links between circulating inflammatory markers and blood methylation. Longitudinal studies with larger sample sizes, using purified cell populations would be needed to replicate associations and establish biological significance. Compelling preliminary findings around SPC₃ associations and candidate genes COL18A1, IFI30 and PPARGC1B provide a foundation to guide future research.

4. Materials and Methods

Study Cohort

Participants included in the BioMood study were nested within the ORIGINS parent cohort, a longitudinal study of family health outcomes, commencing in pregnancy [16,17]. The present study utilized available whole blood samples collected at the 36-week of pregnancy time-point from eligible mothers. The BioMood study entry criteria have been previously published [14]. The pregnant women who were recruited to the BioMood study comprised women who had completed a modified 13-item Mediterranean Diet Questionnaire (MDQ) twice during pregnancy with scores that varied by no more than 2 points. A total of 52 women met eligibility criteria, which included 25 control participants within the Low Mediterranean Diet Adherence (LMDA) group (an MDQ score equal to or less than 4 at 2 timepoints) and 27 case participants within the High Mediterranean Diet Adherence (HMDA) group (an MDQ score equal to or greater than 8 at 2 timepoints). Participants were excluded if they did not have available biological samples for analysis, were exposed to antibiotics at any time during pregnancy, reported use of anti-depressant or anti-anxiety medication or had a pre-pregnancy BMI over 40kg/m². Study groups were balanced for potential confounding variables which include age, parity, education, work status, pre-pregnancy weight and pre-pregnancy BMI.

Serum Inflammatory Marker Analysis

Metabolic profiles were generated using 1H Nuclear Magnetic Resonance (NMR) spectroscopy. Further details on the acquisition, processing and analysis of metabolic profiles are described in the BIOMOOD study [14]. The NMR signal intensities of interest were glycoprotein signals GlycA and GlycB, and supramolecular phospholipid composite (SPC) peaks, which correspond to small HDL phospholipids (SPC₁ peak), larger HDL phospholipids (SPC₂ peak) and LDL phospholipid particles (SPC₃ peak).

DNA Extraction and Library Preparation

Genomic DNA was extracted from whole blood using either chemagic™ DNA Blood 400 Kit H96 (Cat. No. CMG-1091) on the chemagic™ 360 instrument (n = 49), or the DNeasy Blood & Tissue Kit (QIAGEN, Cat. No. 69504; n = 2), following manufacturers’ protocols. DNA concentration and purity were measured using Qubit ™ dsDNA High Sensitivity (HS) Assay Kit (Cat. No. Q33231, ThermoFisher Scientific, USA) on a Qubit fluorometer (ThermoFisher Scientific, USA).200 ng of genomic DNA per sample were fragmented to 275- 325bp using the Covaris E220 sonicator (Covaris, Woburn, MA, USA). Control DNA (0.5% methylated pUC19 and unmethylated lambda DNA) was included to assess enzymatic conversion efficiency.

Indexed pre-capture libraries were prepared using the NEBNext® Enzymatic Methyl-seq Library Preparation Protocol (Twist Bioscience and NEB) following manufacturers protocol. Eight libraries were pooled equimolarly for capture, balanced by diet group and BMI. Target enrichment was performed using the Twist Human Methylome Panel (Code: TWB105520), capturing 3.98 million CpG sites. Hybridization was carried out for 16 hours at 60°C, followed by washing and 6-cycle PCR amplification according to the Twist Targeted Methylation Sequencing Protocol Library size and concentration were assessedusing the Agilent 4200 TapeStation system with the D5000 ScreenTape Assay (Agilent Technologies, Santa Clara, CA).

Sequencing

Final libraries were sequenced by Genomics WA accredited facility on an Illumina NovaSeq 6000 (Illumina, CA, USA) using the XP workflow on an S4 flowcell (2x 150 (bp) reads. with a loading concentration of 253pM and a 10%PhiX spike-in.

Bioinformatic analysis: Pre-processing

Raw methylation FASTQ files were processed using the nf-core/methylseq v2.6.0 pipeline. Methylseq pipeline was executed under Nextflow v23.10.0 [18] and the Human Genome Reference Consortium Human Build 37 (GRCh37) using the bismark workflow. Binary Alignment/Map (BAM) Index files were constructed using samtools version 1.13 [19]. The nf-core methylseq pipeline was run under the profile singularity and a custom configuration file adopted from the pre-configured Pawsey Setonix HPC settings, specifying the compute resources required for the job scheduler (slurm) [20].

Statistical analysis

All statistical analyses were conducted in R language v4.4.1 and RStudio 2023.03.0+386 using base R functions (unless otherwise specified). Processing of Bismark methylation calls was performed using the bsseq package version 1.40.0 [21]. Bismark files were processed to filter any CpGs that were not represented in all samples, and with a coverage of less than 5, or greater than 500 reads. Non-standard and off-target chromosome reads were removed from the set. Sample quality was assessed by examining genome-wide methylation patterns. Raw methylation data were extracted using the getMeth function from the minfi package (v1.50.0) and visualized with densityPlot [22]. Sex was predicted based on median methylation values of loci on chromosomes X and Y, calculated using chrSelectBSseq, getMeth, and colMedians, with females identified as those with less than 0.25 difference in median methylation between chromosomes.

Methylation data, stored as a bsseq object, were annotated to the Twist Human Methylome Panel Target BED File (Twist Bioscience). Promoter regions, CpG islands, and enhancer regions were identified using build_annotations (annotatr package v1.30.0) and subsetByOverlaps [23]. For each region, methylation ratios were extracted using getMeth (bsseq package v1.40.0), and principal component analysis (PCA) was performed using prcomp. Samples identified as outliers in the PCA were excluded from further analysis (4 samples removed from the same capture pool). Methylation beta values were extracted using getMeth. M-values were calculated by log2- transforming the ratio of methylated to unmethylated read counts obtained via getCoverage. Batch effects were assessed by PCA on beta values. Associations between the top five principal components and clinical traits were evaluated using the WGCNA package (v1.72-5) [24]. Spearman correlations and corresponding p-values were calculated and visualized using labeledHeatmap. Pearson correlations between GlycA and GlycB, and among SPC1, SPC2, and SPC3, were assessed and visualized. Normality of metabolite distributions was also assessed using ggplots. SPC metabolites were log-transformed prior to correlation analysis.

Hypothesis testing

Genome-wide differential methylation analysis was performed using the limma package (v3.60.4) [25]. To identify CpG sites differentially methylated between low- and high-MDA groups, logistic regression models were fit M-values and adjusting for the first four principal components. Empirical Bayes statistics were calculated, and top-ranked differentially methylated sites were identified. Effect sizes are reported on the Beta value scale. Differentially methylated CpG sites were deemed statistically significant at unadjusted p values < 5x10-8 and suggestive at unadjusted p values < 1x10-5.

A ‘gene-centric’ analysis was performed by repeating the genome-wide analysis on the methylation dataset after pre-filtering to include only gene promoter regions using the annotatr annotation package(v1.30.0)[23].

An additional candidate gene analysis of CpG sites annotated to genes associated with MD published by Apron et al. [13] was conducted using logistic regression. Genes of interest included EEF2, RUNX3, IL4I1, MAPKAP2, COL18A1, LEPR, PLAGL1, IFRD1, PPARGC1B. Specific CpG sites annotated to these genomic coding regions based on genomic location (hg19 build) were extracted for analysis. The specific regions of interest and corresponding CpGs are available in the published codebase. Statistical significance was declared at a Bonferroni-corrected p value threshold of less than 0.05.

Linear regression models were used to test associations between genome-side CpG and metabolite levels (GlycA, GlycB, SPC1, SPC2, and SPC3), adjusting for MDA status and the first four principal components. As related metabolites were exhibited co-linearity, models were further adjusted by incorporating correlated metabolites as adjustment variables. Differentially methylated sites were visualized using manhattan_plot from the ggmanh package (v1.8.0) [26].