Aflatoxin Exposure during Early Life Is Associated with Differential DNA Methylation in Two-Year-Old Gambian Children

Background: DNA methylation is an epigenetic control mechanism that may be altered by environmental exposures. We have previously reported that in utero exposure to the mycotoxin and liver carcinogen aflatoxin B1 from the maternal diet, as measured using biomarkers in the mothers’ blood, was associated with differential DNA methylation in white blood cells of 6-month-old infants from The Gambia. Methods: Here we examined aflatoxin B1-associated differential DNA methylation in white blood cells of 24-month-old children from the same population (n = 244), in relation to the child’s dietary exposure assessed using aflatoxin albumin biomarkers in blood samples collected at 6, 12 and 18 months of age. HM450 BeadChip arrays were used to assess DNA methylation, with data compared to aflatoxin albumin adduct levels using two approaches; a continuous model comparing aflatoxin adducts measured in samples collected at 18 months to DNA methylation at 24 months, and a categorical time-dose model that took into account aflatoxin adduct levels at 6, 12 and 18 months, for comparison to DNA methylation at 24 months. Results: Geometric mean (95% confidence intervals) for aflatoxin albumin levels were 3.78 (3.29, 4.34) at 6 months, 25.1 (21.67, 29.13) at 12 months and 49.48 (43.34, 56.49) at 18 months of age. A number of differentially methylated CpG positions and regions were associated with aflatoxin exposure, some of which affected gene expression. Pathway analysis highlighted effects on genes involved with with inflammatory, signalling and growth pathways. Conclusions: This study provides further evidence that exposure to aflatoxin in early childhood may impact on DNA methylation.


Introduction
Epigenetic alterations such as DNA methylation changes have been proposed as important mediators of long-term health effects of environmental exposures [1,2]. DNA 2 of 17 methylation patterns are set during early development, where they participate in the establishment and maintenance of cell lineages, but these patterns may be modified during life [3]. Consistent with the role of DNA methylation as the key mechanism for control of gene expression [4], changes in DNA methylation patterns have been observed with aberrant transcriptional changes in cancer and other diseases [5][6][7][8]. There is increasing evidence that environmental exposures, including cigarette smoking, heavy metals and persistent organic pollutants can alter DNA methylation in humans [1,2]. Changes in DNA methylation during foetal development have been associated with adverse birth outcomes, with evidence that these changes can also influence disease in later life [9,10], although the precise underlying mechanism remains poorly understood.
Aflatoxin B1 (AFB1) is a mycotoxin that is highly toxic and carcinogenic [11]. The Aspergillus species of fungi that produce AFB1 commonly grow on food crops, including staples such as maize and peanuts, in warm humid climates, leading to frequent contamination of these crops. At particular high risk of exposure to AFB1 are subsistence farmers in sub-Saharan Africa, where climate conditions, poor storage practices and low socioeconomic status combine to increase intake of contaminated foods. Numerous studies using biomarkers of exposure have shown chronic and occasionally very high levels of AFB1 exposure among rural populations in sub-Saharan Africa [12][13][14][15][16]. Importantly, many studies have linked this exposure to adverse health outcomes, including liver cancer, acute liver disease, and child growth impairment [17,18]. AFB1 is known to induce mutagenic DNA adducts that are likely to be involved in the development and progression of cancer [19]. The mechanism for other health effects have not been determined, but we have previously shown that in utero exposure to aflatoxin, as assessed by levels of aflatoxinalbumin adducts (AF-alb) in serum of the mothers, was associated with differential DNA methylation in a number of genes in the white blood cells (WBC) of the infants [20].
In this study, we examined the impact of early-life exposure to aflatoxin on the DNA methylome patterns, using a new-generation bead-chip arrays with genome coverage. We established the methylome profiles of young children's DNA from Gambian children (who had been naturally exposed to aflatoxin in their diet, as assessed by measuring AF-alb in serum) at two years of age. Further analyses revealed that exposure to aflatoxin during the first two years of life is associated with differential DNA methylation in a range of genes and pathways, some of which directly relate to known health effects of aflatoxin exposure.

Aflatoxin Exposure over Time in Early Life
Characteristics of the children and geometric mean values for AF-alb are given in Table 1. We observed an increase in median values of AF-alb with time (3.7 pg/mg albumin at 6 months, 22.7 pg/mg albumin at 12 months and 47.5 pg/mg albumin at 18 months; Figure 1B), although an increase between 6 and 12 months was not observed in all individuals (Supplementary file S1B).

DNA Methylation Markers of Aflatoxin Exposure and Associated Biological Pathways: Continuous Model
In the continuous model, the Q-Q plot ( Figure 1C) shows limited inflation (lambda close to 1.0), giving a low risk of false positives. Seven differentially methylated positions (DMPs) had a difference in DNA methylation associated with aflatoxin exposure, with an FDR-adjusted p-value < 0.5 ( Figure 1D; Supplementary file S2). The changes in DNA methylation according to AF-alb continuous values for these DMPs are shown in Figure 1E, which also shows that there was no difference according to sex of the children.
In addition to these DMPs, there were 509 DMRs showing differences in relation to continuous AF-alb levels at 18 months (Figure 2A; Supplementary file S3), with a large proportion of these encompassing at least three CpG sites with a consistent direction of effect. Ontology analysis of the results ( Figure 2B) showed that the DMRs are enriched in genes associated with acute hepatic failure, as well as others associated with developmental disorders including intrauterine growth retardation (Supplementary file S4; e.g., micropenis, short neck and intrauterine growth retardation). In terms of biological mechanisms, DNA methylation markers associated with aflatoxin exposure were enriched in inflammatory pathways (particularly the hypermethylated DMRs; e.g., IL4 signalling), nervous system signalling (particularly the hypomethylated DMRs; e.g., alpha adrenergic receptor signalling, EphrinB-EPHB pathway, adrenaline and noradrenaline biosynthesis, and neurotransmitter clearance) and signalling events affecting the growth hormone-insulin-obesity axis (enriched in both hyper-and hypo-methylated DMRs; e.g., leptin insulin overlap, IGF1R, insulin signalling, regulation of PGC-1a, growth hormone receptor signalling, and adipocytokine signaling; Figure 2B). In the continuous model, the Q-Q plot ( Figure 1C) shows limited inflation (lambda close to 1.0), giving a low risk of false positives. Seven differentially methylated positions (DMPs) had a difference in DNA methylation associated with aflatoxin exposure, with an inflammatory pathways (particularly the hypermethylated DMRs; e.g., IL4 signalling), nervous system signalling (particularly the hypomethylated DMRs; e.g., alpha adrenergic receptor signalling, EphrinB-EPHB pathway, adrenaline and noradrenaline biosynthesis, and neurotransmitter clearance) and signalling events affecting the growth hormone-insulin-obesity axis (enriched in both hyper-and hypo-methylated DMRs; e.g., leptin insulin overlap, IGF1R, insulin signalling, regulation of PGC-1a, growth hormone receptor signalling, and adipocytokine signaling; Figure 2B). Hyper-and hypo-methylated CpGs in relation to AF levels are shown in blue and orange circles, respectively, and the size of the circle is proportional to the effect size (being the regression coefficient). The middle zoomed cluster represents DMRs encompassing three or more CpGs showing a 100% consistent direction of effect. (B) Ontology (first row) and pathway (other rows) analysis using EnrichR was performed on all DMRs (first column) and then stratified by hyper-and hypo-methylated DMRs (second and third columns, respectively).
Hyper-and hypo-DMRs were defined as those in which at least 2/3 of the CpGs show a consistent direction of effect towards hyper-and hypo-methylation, respectively. Results are ranked by p-value (the most significant being on top), and only the two most significant pathways (p < 0.05) are shown for each database. ** FDR < 0.05; * FDR < 0.1. More results for Human Phenotype Ontology of All DMRs are shown in Table S3A as more than two significant (FDR < 0.1) results were obtained. In ontology analysis, only Human Phenotype Ontology is reported here as it yielded the most significant results; GO Biological Process did not yield significant ontologies (FDR < 0.1). A detailed report for the ontology and pathway analyses is available through the following links: All DMRs: https: //amp.pharm.mssm.edu/Enrichr/enrich?dataset=b4e3a404c15ef001a7c30c1a13e1d249#; Hypometh DMRs: https://amp.pharm.mssm.edu/Enrichr/enrich?dataset=9357fa44b29ebb8dc326803ebf26925e; and Hypermeth DMRs: https://amp.pharm.mssm.edu/Enrichr/enrich?dataset=fcbe92f8726688e4 eee8961de6e89d03 (Accessed on 23 March 2020).

DNA Methylation Markers of Aflatoxin Exposure and Associated Biological Pathways: Categorical Time-Dose Model
In the categorical analysis, the Q-Q plots also show low inflation, with lambda values being close to 1.0 ( Figure 3A,C,E). The largest number of hits was seen for the recurrent vs. low exposure model, followed by the recurrent vs. intermittent and intermittent vs. low models ( Figure 3B,D,F), suggesting a time-dose response. The DMPs (FDR < 0.05) of the recurrent vs. low ( Figure 3B) and recurrent vs. intermittent ( Figure 3D) exposure models are further detailed in Supplementary files 5 and 6, respectively. No DMPs (FDR < 0.05) were obtained in the intermittent vs. low exposure model ( Figure 3F). The methylation levels for each category (low, recurrent and high) are shown in Figure 3G, visualizing the effect of the time-dose response of exposure on methylation. Some sex-specific effects are seen for DMPs: in UBQLN2, WIPI2, PRIM2 and HRH3, the effect was only seen in girls, and in SULF2, boys showed an inverse (though minor) direction of effect relative to girls.  There were 714 DMRs obtained in the Recurrent vs. Low AF-alb model ( Figure 4; see also Supplementary file S7). 51 of these encompass at least one CpG with an effect size ≥ 5% (Supplementary file S8). The top four DMRs are shown here ( Figure 4A). The DMRs in the MIPEPP3 and HLA-DPB1 genes represent the two DMRs with the largest number of CpGs, each having an effect size ≥ 5%. Among the remaining DMRs that pass this criterion, those in the CCKBR and LOC494141 genes are shown because they encompass at least three consecutive CpGs with an effect size ≥ 5%. Except for HLA-DPB1, the other three DMRs show a time-dose response, with the effect sizes of intermittent exposure being in between those of the recurrent and low exposures. For HLA-DPB1, intermittent and recurrent exposures had almost identical methylation levels highlighting a potentially irreversible hypomethylation effect of AF on this gene, relative to low exposure. Other DMRs with less stringent criteria are shown in Supplementary file S9 (DMRs encompassing at least two consecutive CpG with effect size ≥ 5%, a max of one CpG being non-unimodally distributed, and 100% consistent direction of effect among its CpGs that have an effect size ≥ 3%). The 714 DMRs, particularly the hypermethylated regions, were enriched in genes associated with circulatory and heart development as well as cardiac-related pathways ( Figure 4B; see also Supplementary file S10; e.g., circulatory system development, heart development, hypertrophic cardiomyopathy, cardiac hypertrophic response, and cardiac progenitor differentiation). Similar to the pathways enriched in DMRs from the continuous AF-alb model, the DMRs for the recurrent vs. low model were also enriched in nervous system-(particularly the hypomethylated DMRs; e.g., axon guidance and GABA-R receptor II signalling), inflammatory-(particularly the hypomethylated DMRs; e.g., IL8-and CXCR1-mediated signalling) and adipose-(enriched in both hyper-and hypo-methylated DMRs; e.g., genes targeted by miRNAs in adipocytes, visceral fat deposits and metabolic syndrome, eicosanoid biosynthesis, and regulation of fat cell differentiation) related pathways. The most enriched ARCHS4 kinase pathway, LIMK1, associates with both nervous system (stimulate axon outgrowth) and cardiac-(predisposing to Supravalvular Aortic Stenosis) related pathways (GeneCards) (see also Supplementary file S11).     Figure 5B). DMRs from either model are also similarly enriched in promoter, 1-5 Kb to TSS, 5 UTR and first exon regions and depleted from exons, introns and intergenic regions relative to the 450K reference set. The DNA methylation at certain CpG loci affects RNA expression of nearby genes (FDR < 0.05), with the direction of effect being positive for some eQTMs and negative for others ( Figure 5C,D). ATP6V1G2 expression was associated with differential methylation of cg23965019 (HSPA1B gene) in both the continuous ( Figure 5C) and categorical time-dose ( Figure 5D) models. CpG methylation and RNA expression of nearby genes (within ± 250 Kb). The CpGs yielded from the DMR continuous and Recurrent vs. Low exposure models were used in (C,D), respectively. In both models, the eQTMs were adjusted for AF levels in order to eliminate any potential direct effect of AF on expression.

Study Population and Blood Sampling
The current work was part of a study on the impact of aflatoxin and child growth, which was a sub-study embedded within the Early Nutrition and Immune Development (ENID) trial (ISRCTN49285450). Pregnant mothers were recruited into the ENID trial from the rural, West Kiang region of The Gambia. Full details of the ENID trial protocol [21] and aflatoxin and growth sub-study [18] have been described elsewhere.
Samples used in the current DNA methylation study were from 243 children born into ENID between May 2011 and December 2012. Children were followed up over two years, with blood samples taken at 6, 12, 18 and 24 months of age. For the purpose of the current analysis we used serum extracted from samples at 6, 12 and 18 months and WBC extracted from samples taken at 24 months.

AF-alb Analysis
The analysis of AF-alb followed the method first described by Chapot and Wild [22]. Albumin was extracted from serum samples, hydrolysed using pronase and purified with Sep-pak C-18 cartridges (Waters, Milford, MA, USA). AF-alb was quantified using a competitive ELISA [22,23] using a rabbit anti aflatoxin-Cl2-BSA polyclonal antibody (a gift from Christopher Wild, IARC) as the primary antibody. The detection limit (LOD) was 3 CpG methylation and RNA expression of nearby genes (within ± 250 Kb). The CpGs yielded from the DMR continuous and Recurrent vs. Low exposure models were used in (C,D), respectively. In both models, the eQTMs were adjusted for AF levels in order to eliminate any potential direct effect of AF on expression.

Study Population and Blood Sampling
The current work was part of a study on the impact of aflatoxin and child growth, which was a sub-study embedded within the Early Nutrition and Immune Development (ENID) trial (ISRCTN49285450). Pregnant mothers were recruited into the ENID trial from the rural, West Kiang region of The Gambia. Full details of the ENID trial protocol [21] and aflatoxin and growth sub-study [18] have been described elsewhere.
Samples used in the current DNA methylation study were from 243 children born into ENID between May 2011 and December 2012. Children were followed up over two years, with blood samples taken at 6, 12, 18 and 24 months of age. For the purpose of the current analysis we used serum extracted from samples at 6, 12 and 18 months and WBC extracted from samples taken at 24 months.

AF-alb Analysis
The analysis of AF-alb followed the method first described by Chapot and Wild [22]. Albumin was extracted from serum samples, hydrolysed using pronase and purified with Sep-pak C-18 cartridges (Waters, Milford, MA, USA). AF-alb was quantified using a competitive ELISA [22,23] using a rabbit anti aflatoxin-Cl2-BSA polyclonal antibody (a gift from Christopher Wild, IARC) as the primary antibody. The detection limit (LOD) was 3 pg AF-lysine equivalents per mg albumin, and each ELISA included three positive and one negative control samples for quality control. All samples were measured in duplicate for each ELISA batch and repeated at least two times on separate days to confirm the results. Results were accepted when values within each ELISA had a %CV below 10% and samples tested on separate occasions had a %CV below 15%. If samples lay above the linear part of the standard curve, the ELISA was repeated at a suitable dilution. This ELISA has been validated for measuring aflatoxin albumin biomarker levels against dietary aflatoxin intake, and also against an LC-MS for measuring aflatoxin lysine. Whilst the ELISA consistently gives Af-alb values that are approximately three-fold higher than aflatoxin lysine values measured by LC-MS the results between methods correlate for the same samples [24]. The AF-alb values are, therefore, suitable for use in analyses to compare variations in aflatoxin exposure in relation to health outcomes or other variables such as differential DNA methylation.

Bisulfite Conversion
Bisulfite conversion was performed using the EZ-96 DNA Methylation kit (Zymo Research Corporation, Irvine, CA, USA), and its efficiency was tested using PCR primers that were specific either to bisulfite-converted or to non-modified DNA and spanning the NBR2 region.

Genome-Wide Methylation Array
Sample preparation and hybridization to the HM450 BeadChip arrays (Illumina, San Diego, CA, USA), which covers~450,000 CpG sites, were performed at the International Agency for Research on Cancer (IARC). From the 600 ng of bisulfite-converted DNA per sample, 200 ng was used for hybridization on the array. Each array consisted of 96 samples distributed equally among 8 chips in a completely randomized manner so that batch effects (e.g., sample position and intra and inter-variability in arrays and chips) do not completely confound with biological covariates of interest. The analysis on the bead array was conducted following the recommended protocols for amplification, labelling, hybridization and scanning [25].

DNA Methylome Data Pre-Processing
Data pre-processing and analysis were performed using R/Bioconductor packages (http://www.bioconductor.org/, accessed on 27 December 2019). Raw data of the samples constituted of 485,512 CpG probes from the HM450 array and were loaded from IDAT files and handled in the minfi package [26] to calculate the methylation level at each CpG as the beta-value (β = intensity of the methylated allele (M)/(intensity of the unmethylated allele (U) + intensity of the methylated allele (M) + 100)). Data quality was assessed using box plots for the distribution of methylated and unmethylated signals, unsupervised clustering showing close agreement between technical replicates, array-wide multi-dimensional scaling plots and unsupervised clustering to identify sample outliers, and multidimensional scaling plots clustering samples into two groups representing infant sex. Filtering of cross-reactive probes (n = 14,549) [27] and low quality probes (n = 1183) (more than 5% missing values per probe) resulted in a total of 469,853 CpGs to analyze. Beta values were normalized with Functional normalization (Funnorm, minfi R package) that was shown to perform equally good or outperform existing normalization methods [28]. It removes unwanted variation by regressing out variability explained by the control probes present on the array. Finally, we use surrogate variable analysis (SVA) [29] for batch effect correction and adjustment for latent variables, a choice validated by the findings of our benchmarking [30]. SVA also increases statistical power by removing (unwanted) variability through aggregating information at the data level and constraining the data's variability to the phenotype of interest [31]. Of 262 samples initially loaded onto the DNA methylation array, five data files were missing and one AF-alb result was missing. After visual inspection of density plots, one sample was removed for displaying an abnormal distribution. In addition, six samples were replicates and six samples were found to be mismatched for sex. This left 243 samples for further analysis.

Assessment of Cell Mixture Distribution
Quantile-normalized data were used to predict WBC subtypes based on a peripheral blood reference panel [32] using the minfi package [26]. Proportions of six cell types (B, CD4+ T, CD8+ T, granulocytes, monocytes and NK) were determined. SVA was also used as a reference-free method to adjust for differences in WBC composition [33].

Statistical Analysis, Dimension Reduction and Region-Level Methylation Analysis
Methylation beta-values (i.e., the ratio of the methylated probe intensity and the overall intensity), ranging from 0 to 1, were used for data plotting and (biological) interpretation. A zero-centered logarithmic transformation (M-value, calculated as the log2 ratio of the intensities of methylated probe vs. unmethylated probe) was used for most statistical analyses since they are more statistically valid, particularly in differential methylation analysis [34]. M-values were interrogated for association with Af-alb levels by modelling the study variables and potential confounders together with latent surrogate variables in a robust linear regression, using the R software. The distribution of Af-alb values at 6, 12 and 18 months was not Gaussian (p-value < 2.2 × 10 −16 for all 3 times points) and we performed a logarithmic transformation obtaining data conforming more closely to the normal distribution (p-value = 0.07 and 0.002 for time point 12 and 18 months respectively) ( Figure S1A). Statistical models for association of AF-alb with DNA methylation were adjusted (methylation = AF-alb levels + child sex + WBC types + Sample_Plate + Ethnicity + Season of sample collection), wherein Sample_Plate represents the HM450 batch of 96 samples each ( Figure S1C). A dimension reduction approach (DMRcate) was implemented to reduce the~450,000 individual sites into clusters of genetically proximal and correlated CpGs to enhance statistical power [31] and aid biological inferences (as single CpG sites often have more subtle functional relevance than CpG clusters), as per our previous work [35,36]. We used the recommended DMRcate proximity-based criteria and a minimum of two methylation sites per region with a maximum gap of 500 bp. To reduce the number of false positive results, we controlled for multiple testing using Benjamini Hochberg (FDR) procedure [37]. Identified differentially methylated regions (DMRs) were are reported if their FDR-adjusted p-values were < 0.05.
Two approaches to analysis of data were taken (Supplementary files S1C,D). Firstly, the continuous model assessed the association between DNA methylation and AF-alb levels from 18 months (at which the highest levels of exposure were observed and which is closest in time to the sampling of white blood cells for methylation analysis). Secondly, the levels of AF-alb were categorised depending on time and dose, specifically based on results for measurements at 6, 12 and 18 months, on whether levels were above or below the median, and whether this changed between time points. This categorical time-dose model consists of three categories: "low", where AF-alb was low at each timepoint; "recurrent", where AF-alb was high during the last two time points (12 and 18 months); and "intermittent", where AF-alb was low or high at different time points (Supplementary files S1D).

Discussion
In this manuscript, we report changes in DNA methylation associated with aflatoxin exposure in WBC of Gambian children who are naturally exposed to aflatoxin in their diet. AF-alb is an established biomarker of aflatoxin exposure that reflects dietary intake of aflatoxin [12,38]. Chronic aflatoxin exposure in this population is known to be high [39], largely due to the inclusion of home-farmed peanuts as a dietary staple, and the impact of the exposure on child growth and liver disease has been reported in this population and others [18,[40][41][42]. We have previously shown that exposure to aflatoxin in utero, as assessed by measurement of maternal AF-alb during pregnancy, is associated with differential DNA methylation in the WBC of Gambian infants at 6 months old [20]. Here, we have observed changes in DNA methylation in WBC of infants from the same population (but not the same children) at 2 years old, associated with levels of AF-alb during early life. The changes in DNA methylation were observed across a number of genes and pathways, with associations reported according to two methods of modelling the aflatoxin exposureusing a continuous model of AF-alb at 18 months (six months prior to blood sampling for epigenetic analysis) and a categorical model that reflects both time and dose response for aflatoxin exposure from six to eighteen months. In each case, analyses revealed DMPs and DMRs associated with AF-alb. Reassuringly, there is a notable amount of overlap between differentially methylated sites and regions across the analytical models.
Most of the associations observed were for small differences in DNA methylation, so care must be taken in attributing biological significance to the findings. Nevertheless, it is of interest to consider the genes and pathways in which differential DNA methylation was observed. In the model using AF-alb as a continuous variable, seven DMPs were observed. Four of these genes code for enzymes (SPR, an aldo-keto reductase; SESN1, a member of the sestrin family that play a role in the cellular response to oxidative stress; ME3, malic enzyme, which catalyzes the oxidative decarboxylation of malate to pyruvate using either NAD+ or NADP+ as a cofactor; CAPN5, a calcium-dependent cysteine protease involved in signal transduction in a variety of cellular processes), of which three are involved in response to oxidative stress. The others are BIRC8, an IAP-like protein 2 that is involved in pathways in apoptosis signalling and TNF signalling pathways; CBLN2, which may play a role in induction of synaptogenesis; and MAP7, a microtubule associated protein predominantly expressed in cells of epithelial origin and thought to be involved in microtubule dynamics.
A large number of DMRs were obtained (509), of which a large proportion encompassed at least three CpG sites, strengthening the reliability of these results. It is of interest that many DMRs were enriched in genes associated with hepatic failure, as the liver is the main target organ for the effects of aflatoxin, including acute poisoning and liver cancer [43]. This also highlights the ability of epigenetic markers to capture the biological mechanisms by which aflatoxin exposure leads to liver disease. Aflatoxin has also been associated with child growth impairment in several studies [18,38] and immune system effects [44]. In our data, aflatoxin-associated DMRs were enriched in inflammatory pathways and signalling events affecting the growth hormone-insulin-obesity axis, as well as Wnt signalling (involved in embryogenesis). Some of the DMRs relate to developmental disorders associated with reduced growth, such as brachydactyl syndrome, short neck and intrauterine growth retardation. Indeed, aflatoxin association in utero has been associated with impaired birthweight and growth in the infant [45]. Other DMRs observed from our data are in pathways involved in nervous system signalling events, blood clotting and morphogenesis.
The categorical time-dose model was used to explore whether changes in recurrence of exposure level over time affected the DNA methylation patterns. The children in Gambia are largely exclusively breast fed until five to six months of age [46]. During this period there may be some exposure to aflatoxin in the form of the metabolite aflatoxin M1 excreted in breast milk, which is an oxidised, and much less toxic metabolite of aflatoxin B1 [47]. Consequently, AF-alb levels tend to be very low in infants at 6 months, increasing with age as the children eat more family food [48]. This pattern is reflected in our data. However, as aflatoxin contamination of food is heterogeneous and will vary from family to family over time, and by season, levels of AF-alb in blood may be high in children at one sampling point and low at another, or vice versa. This complicates the interpretation of the role of the exposure in modifying DNA methylation, and it is not clear which time points may have a greater effect. The categorical time-dose model, based on the categorical split of AF-alb levels into high vs. low (based on the median value) attempts to take this variation into consideration. A small number of these DMPs showed sex related differences, but in general the results for boys and girls were similar in the continuous and categorical time-dose models.
The categorical time-dose model for recurrent vs. low AF-alb resulted in over 700 DMRs, of which 51 involved more than one CpG site with an effect size of more than 5%. Ontology analysis showed that DMRs were enriched in genes associated with cardiac and circulatory development. These genes were characterised by hypermethylation, whereas DMRs in genes from inflammatory and nervous system genes were mainly hypomethylated. When comparing the results between the models (Figure 5), the fact that the recurrent vs. low model (with the highest difference in AF-alb) yields a larger number of CpGs (FDR < 0.05) than the Recurrent vs. Intermittent or the Intermittent vs. Low models, is indicative of a time-dose response, which is further supported by a time-dose effect size observed for many of the yielded DMRs.
In our previous paper exploring epigenetic changes in relation to aflatoxin exposure in The Gambia, we compared DNA methylation in WBC of infants with AF-alb levels in serum of the mothers during pregnancy. We found differential DNA methylation in 71 DMPs associated with levels of maternal aflatoxin exposure. Of these, a small number were also found in this study. PRKAR1B, coding for a regulatory subunit of cAMP-dependent protein kinase A, was among genes in the DMRs for the Recurrent vs. Low AF-alb model and was also a DMP in this model. Mutations in the PRKAR1B gene cause neurodegenerative disorder [49]. Other genes that were DMP hits in the earlier paper that were found from one of the analytical models in the current study were TGFB1 (encodes a secreted protein involved in cell growth and proliferation), HAND2 (encodes a transcription factor essential for cardiac morphogenesis), ZBTB16 (encodes a zinc finger transcription factor involved in cell cycle progression), TNXB (encodes a extracellular matrix glycoprotein), AGAP3 (encodes an essential component of the NMDA receptor signalling complex) and DNAJB6 (encodes a molecular chaperone involved in protein folding and assembly). Although only 10% of the genes from the earlier paper were hits in this study, it is interesting that these genes have shown differential DNA methylation in both 6-month-old children associated with aflatoxin exposure in pregnant mothers and in 24-month-old children, associated with aflatoxin exposure during early life.
A number of previous studies have reported changes in DNA methylation in relation to AFB1 treatment of human cells in vitro [50][51][52][53]. When primary human hepatocytes were treated with AFB1, Rieswijk et al. found six genes that showed persistent hypomethylation associated with up-regulated gene expression. One of these genes, TXNRD1 (encodes an oxidoreductase) was a DMP in the recurrent vs. low AF-alb model. Wang et al. reported that increased methylation of the RUNX3 gene was integral to the transformation of AFB 1 treated LO2R cells. Methylation of a DMR including this gene was also seen to be associated with AF-alb in our Recurrent vs. Low AF-alb model. Modulation of the expression of the DNA methyl transferases, DNMT1 and DNTM3a, is a possible mechanism by which AFB1 may influence DNA methylation, leading to altered gene expression [54].
A limitation of the current study is that data on aflatoxin exposure of these infants in utero was not available. Our previous study on 6-month-old infants from the same Gambian population showed that maternal exposure to aflatoxin is associated with differential DNA methylation in the children, so it is possible that maternal exposure patterns may have contributed to DNA methylation in these children as well. Exposures at particular time points, such as in utero, may have a larger effect on DNA methylation changes induced by AFB1 than at other times. Nevertheless, it is interesting that 10% of genes with DMPs identified in the earlier study also showed up as affected in this dataset. We have taken into account several potential confounding factors in our analytical models. However, there will be other variables that could have impacted on DNA methylation of which we were not aware, or for which we do not have data. This is one of the reasons why we additional corrected for potential unknown effects using SVA (see Methods). Despite this, there is good evidence that aflatoxin B 1 exposure in early life does modulate DNA methylation, and with several genes showing up in previous studies in vitro and in vivo.
The high levels of aflatoxin exposure in children in Gambia reported here and in other populations in low-and middle-income countries, highlight the need for aflatoxin reduction strategies, which have been discussed elsewhere [55].
In summary, our results reveal that dietary exposure to aflatoxin B 1 during the first two years of life is associated with differential DNA methylation in a range of genes and pathways, some of which directly relate to known health effects of aflatoxin exposure. It is not clear whether the differences in DNA methylation associated with aflatoxin exposure are biologically significant. There was an apparent time-dose relationship in response and overlap of hits between analytical models. This provides further insight into potential mechanisms underlying the impact of exposure to aflatoxin on human health in early life periods. Further research into the role of DNA methylation on health effects of aflatoxin is warranted.