The effect of genotype and in utero environment on interindividual variation in neonate DNA methylomes

1The relationship between interindividual variation in the epigenome—especially DNA methylation—and disease risk is an area of intense research interest. Although the effect of fixed genetic variation on DNA methylation is apparent in studies of allele-specific methylation and genomic imprinting, there is also emerging evidence for environmental influences as a source of epigenomic variation. Perinatal cohort studies offer a unique opportunity to explore the origins of variation across the epigenome, and in particular the extent to which fixed genetic variation can moderate the relationship between prenatal environmental factors and epigenetic status at birth.

2Epidemiological data link disease risk directly to the in utero environment or to birth outcomes as a surrogate for the in utero environment. This phenomenon is often called fetal programming and defines, in part, the developmental origins of health and disease.

3Stable alterations to the epigenome are considered to be a putative molecular mechanism for fetal programming. Thus the environmental epigenetic hypothesis suggests that the in utero environment affects the epigenome and that resulting epigenetic marks alter physiology to affect later disease risk. Evidence for this hypothesis derives from studies documenting the relationship between specific DNA methylation marks, in utero environments, and later phenotypes. For example, maternal carbohydrate intake in early pregnancy and offspring adiposity at 9 yr of age associate with the DNA methylation level of the RXRA promoter in the umbilical cord. There are examples of persistent, environmentally induced DNA methylation states. For instance, variations in maternal care alter the methylation of the glucocorticoid receptor gene promoter in rats, and differential methylation at this same region is associated with childhood trauma in humans. Although many aspects of the epigenome such as histone modification are likely to be involved, current evidence concentrates mostly on DNA methylation.

4Multiple epigenome-wide association studies (EWAS) within birth cohorts have been initiated to search for epigenomic signatures of early life environment that may influence later life phenotype. However, at least some DNA methylation marks are specified by sequence context in cis, precluding exclusive environmental regulation. In humans, interindividual variation in DNA methylation could be wholly or in part a consequence of nucleotide polymorphism, as local genetic variants influence the propensity for methylation of neighboring cytosines. These polymorphisms can be defined as methylation quantitative trait loci (methQTLs). MethQTLs have been postulated as a link between GWAS SNPs and phenotype.

5showed methQTLs in human tissue in studies of multiple human brain regions and found between 4% and 5% of the 27,000 CpG sites studied had methylation levels that were significantly dependent on genotype (as measured in 537,411 SNPs). Although both cis (defined as the CpG and SNP pair being within 1 Mb of the chromosomal region) and trans effects were detected in this study, cis methQTLs were in massive abundance. The average distance between a SNP and CpG pair that made up a significant methQTL was 81 kb, but the peak enrichment across the cis methQTLs was 45 bp. examined human cerebellum and found 9% of CpGs tested (748/8590) were within methQTLs. Another study detected 180 (1%) methQTLs in the 22,290 CpGs and 3 million SNPs investigated across 77 HapMap lymphoblastoid cell lines, five of which were also reported by. Most methQTLs were found in cis and over short distances of <5 kb. Moreover, genotype at one SNP associated with methylation at multiple neighboring CpGs, as might be expected given the positional correlation previously noted in methylation data. Similar results were obtained by, who also discovered methQTLs in adult adipose tissue acting in cis, which explains 19% of the observed variance in methylation levels.

6Ethnicity can be used as a proxy for genotype and has been shown to influence the DNA methylome. African and European individuals have population-specific patterns of DNA methylation at ∼30% of CpGs measured. Methylation levels at ∼50% of these population-specific CpGs are explained by divergence in allele frequencies at cis-acting SNPs between populations. Studies of methylation differences using Illumina InfiniumHumanMethylation450 BeadChip array data from 133 lymphoblastoid cell lines from European and African HapMap samples found that 13% of analyzed CpGs showed significant methylation differences between the populations, >50% of which were in methQTLs with local SNPs. CpGs showing differential methylation levels across ethnicities are more likely to be driven by genotype than other CpGs. Nevertheless, these studies also found that some methQTLs are specific to one population with no correlation between genotype and methylation in other populations, suggesting possible gene × environment interactions.

7Methylomes are more similar in related than unrelated individuals, and concordance tracks degree of relatedness; the methylomes of monozygotic twins are more closely related than those of dizygotic twins. Methylation profiles from three different tissues of twin neonates generated on the InfiniumHumanMethylation27K BeadChip only clustered into twin pairs between 29% and 71% during unsupervised analysis, suggesting nongenomic influences on the newborn methylome.

8In this article, we focus on the examination of the relative influences of genotypic, environmental, and gene × environment interactive effects on the neonatal methylome. Recent studies describe evidence for gene × environment interactions (G × E effects) on DNA methylation. found that the LEPR genotype interacted with maternal smoking to associate with methylation of LEPR. A SNP within the IL4R gene combined with methylation at a CpG site within the same gene predicts the risk of childhood asthma. Moreover, found that interaction of the FKBP5 genotype and early childhood trauma affects methylation of FKBP5 intron 7, FKBP5 expression, and subsequent deregulation of glucocorticoid receptor signaling. The proportions of interindividual variation in methylomes that are driven by genotype, environment, or an interaction of gene and environment (G × E) are currently unknown. To clarify the relative influence of gene and in utero environment on epigenetic status at birth, we studied the variation in genome-wide DNA methylation patterns in umbilical cord samples from 237 Asian neonates using the InfiniumHumanMethylation450 BeadChip together with genotyping and extensive measures of in utero environmental conditions. We report that genotype, and in particular G × E interactions, explain substantial proportions of interindividual variation in the methylome at birth.

9Umbilical cord tissue DNA from 237 individuals (131 Chinese, 72 Malay, 34 Indian) in the GUSTO birth cohort were interrogated on both Illumina OmniExpress + Exome genotyping arrays and InfiniumHumanMethylation450 BeadChip; 708,365 SNPs (from the 958,178 assayed) varied in genotype and 301,468 CpGs (from the 411,107 assayed) varied in methylation levels by >5% across the 237 individuals under study (for study subject characteristics, see Table 1). When the genotype data were subjected to principal component analysis, the samples were cleanly separated by ethnicity. The Indian subjects clustered away from Chinese and Malay subjects on principal component 1, and the Chinese and Malay subjects separated on principal component 2 (Fig. 1A). In contrast, when the DNA methylation data were subjected to principal component analysis, the samples did not separate well by ethnicity on components 1 or 2 (Fig. 1B). The methylation and genotyping arrays do measure different subsets of the genome; however, each have probes in every coding gene in the genome and so are broadly comparable. If methylation levels were specified wholly by genotype, we would expect ethnicity to drive the methylome in the same manner as observed for the genotype.

10Methylation levels at the majority of CpGs in the neonate methylomes were very similar across individuals, with median absolute deviation scores (MAD) <0.1. To identify regions with appreciable interindividual variance that could reflect individual in utero experience or genotype, we defined interindividual variably methylated regions (VMRs) across the 237 individuals using a previously published methodology. This analysis returned 1423 VMRs; these were found on every chromosome in a manner roughly dependent on probe spacing (Supplemental Table 1; Supplemental Fig. 1), although they tended toward north shore (P = 1.77 × 10−4) and open sea (P = 9.11 × 10−12) regions (Supplemental Table 2).

11The CpG with the highest MAD score within each VMR was used to represent each of the 1423 VMRs and denoted as VMR-CpGs. For each VMR-CpG, we compared methylation levels to SNP genotype at all 708,365 heterologous positions by linear regression, controlling for sex. Strong sex effects on the autosomes have previously been noted in InfiniumHumanMethylation450 BeadChip data and are at least partially artifactual, driven by cross-reaction of probes with the sex chromosomes.

12The best-match (lowest P-value) SNP for each VMR-CpG was retained. The final data set of 1423 CpG–SNP pairs showed a range of associative P-values skewed toward the low end, suggestive of methQTLs in the data set (Supplemental Fig. 2). P-values obtained from the regression analysis were binned into 1000 equally spaced bins. Distributions were defined as skewed if the first bin contained more than the 708 P-values expected for each bin if P-values were distributed evenly. Nine hundred sixty-six (68%) of the individual VMR-CpGs had a skewed P-value distribution. In addition, 1037 (73%) had a –log10 P-value above 6.5, which is approximately the background noise level seen across the genome (Fig. 2B,D,F).

13Twelve of the 1423 CpG–SNP pairs (<1%) included a SNP located within the CpG that either creates or eliminates a CpG site. We denoted these instances as “disrupting pairs.” Eight hundred twenty-eight of the 1423 pairs (58%) included a SNP and CpG from the same chromosome. We denoted these as “cis pairs.” The remaining 583 pairs (41%) include a CpG and SNP located on different chromosomes and were denoted as “trans pairs.” The proportion of cis pairs was much greater than expected by chance (chance would predict pairs to be equally distributed across chromosomes at ∼4.5%) and tended to be more likely to come from a skewed P-value distribution and to have higher R2 values than the trans pairs (Fig. 2A). The influence of a genotype operating independently of environmental context was undetectable for some VMR-CpGs (e.g., Fig. 2B,C) and very clear for others (e.g., Fig. 2F,G). Even in a mid-range where methylation at the VMR-CpG was not strongly related to genotype on a scatter plot (e.g., Fig. 2E), there was still a tendency for the best hit to be in cis (e.g., Fig. 2D).

14Disrupting pairs had significantly stronger associations between genotype and methylation level than cis pairs, which in turn had stronger associations than trans pairs (P = 7.3 × 10−117 by Kruskal-Wallis test) (Fig. 3). Within the cis pairs, there was a strong inverse relationship for the strength of association between methylation levels and genotype, and the distance between the CpG and SNP (P = 8.84 × 10−5). The mode was within 0–10 bp or 50–60 bp without the disrupting pairs (Fig. 4). There seemed to be a continuum of genotypic influence on methylation levels, with ∼68%–73% of the VMR-CpGs showing an appreciable association with genotype. Defining a genuine methQTL would necessitate choosing arbitrary cut-offs. However, the trans CpG–SNP pairs tended to be on the lowest end of the distribution of association statistics and hence were less likely to represent genuine methQTLs (Fig. 2A).

15The GUSTO study has ascertained multiple environmental measures that are known to influence or reflect the quality of the in utero experience. We chose 19 parameters as surrogate measures of the uterine environment (see Table 1) and examined whether methylation at each of our 1423 VMR-CpGs was better explained by genotype, environment, or G × E. We ran 39 regression models (see Methods), including the best SNP hit from the previous methQTL analysis, described above. For the G × E models, we tested all the SNPs located on the same chromosome as the CpG, as well as the best hit from the methQTL analysis. All models contained sex and were compared using Akiake information coefficients (AICs). Methylation levels at ∼25% of the 1423 VMR-CpGs were best explained by genotype alone, while the rest were best explained by G × E models (Fig. 5A). The models containing environment alone were never the best explanation of methylation at the VMR-CpGs. The information loss experienced between the two top models (Δ) ranged from 0.0002–49.0 (calculated by delta in AICs). The genotype-only model tended to be a “narrow winner”; i.e., the Δ was low in comparison to the Δ shown by the cases where G × E was the best explanation (Fig. 5B). The majority of best models showed an adjusted R2 value above 0.12. The distribution of adjusted R2 values for the cases where G × E was the best model was shifted slightly to the right compared with that for the cases of genotype being the best model (Fig. 5C). When the results were restricted to models without substantial support for the next best model, with Δ > 2 (Fig. 5D), and an adjusted R2 > 0.4, the proportion of VMR-CpGs for which genotype only was the best model was still ∼24%.

16The VMR-CpGs that were best explained by genotype alone tended to be in open seas (P = 1.38 × 10−5), while VMR-CpGs that were best explained by G × E tended to be both in open seas (P = 2.12 × 10−6) and north shores of CpG islands (P = 0.012) (Supplemental Table 2; Supplemental Fig. 1), consistent with previous observations. Once again the majority of the SNPs in the best models were in cis with the VMR-CpGs (85% for the winning G × E models and 78% for the G winning models). When only models with high levels of support (Δ > 2 and adjusted R2 > 0.4) are included, only cis SNPs remain (Supplemental Fig. 3).

17The G × E models allow for methylation associated with environment in each of the genotypic subgroups present but with different slopes. However, we often saw that environment and methylation were very closely associated in one genotype and less so in the other two genotypes. To identify these CpG-VMRs, we segregated the subjects by the genotype in the best G × E model and ran regressions of methylation and the phenotype for the best G × E model. Fifty VMR-CpGs were significant after Bonferroni correction. Examples are shown in Figure 6.

18The 491 VMR-CpGs with differential methylation for ethnicity have significantly stronger associations with SNPs compared with all other 930 CpGs (P = 2.17 × 10−19). We reasoned that the inclusion of multiple ethnicities would increase the genotypic influence on the methylome. We thus conducted an analysis including only the 131 Chinese subjects in our data set. In the Chinese-only subgroup, we still found many methQTLs, but there was a subtle (downward) shift to less significant associations for the majority of CpG–SNP pairs (Supplemental Fig. 4A). The number of VMR-CpGs for which the methylation was best explained by genotype alone decreased slightly to 21% (Supplemental Fig. 4B).

19We noted that the CpG MAD score across samples or the range of methylation values across samples was related to the strength of association between the CpG methylation values and the genotype of the best SNP (P = 3.54 × 10−39) (Fig. 7). Our 1423 VMR-CpGs were chosen to lie within VMRs. This approach improves the specificity of the analysis. However, we also picked a further 1500 CpGs that had the most extreme MAD scores in the data set and that were not included in a VMR. Indeed, within this set the genotype-driven model is a better explanation for a higher proportion of the CpGs (Supplemental Fig. 5).

20We used DNA obtained from the umbilical cord at birth to provide evidence for the influence of genotype on interindividual variability in DNA methylation. MethQTLs were readily apparent in the methylome of neonates born at term with birth weights that spanned the normal range from a community sample of pregnancies that were not selected or weighted for any particular outcome. Our findings are consistent with a recent report of 14,189–32,318 methQTLs (3%–7%) found within Infinium450K data by searching the SNPs (genotype 2.5 million) in 5-kb flanking regions on either side of each CpG in umbilical cord tissue and cord blood. However, did not describe any phenotypic or environmental data on their subjects or attempt to show any environmental influence on methylation levels. In these ways, our studies differ in aims and scope. We found that the effect of genotype on methylation was more pronounced in multiethnic populations but was also apparent within the sample of homogeneous ethnicity (Fig. 5; Supplemental Fig. 4B). Moreover, there was a strong relationship between the range of methylation values and influence of genotype (Fig. 7; Supplemental Fig. 5). There is recent evidence that epigenetic states may serve to directly mediate the relation between certain genetic polymorphisms and phenotype. The mechanism by which sequence polymorphism affects local CpG methylation is unknown (assuming the SNP is not within the CpG, a situation we have labeled disrupting pairs). However, the effect is likely mediated by the sequence-specific binding of proteins and alteration of chromatin factors that subsequently affects binding of de novo methylases or demethylating mechanisms. A similar mechanism has been evidenced for other epigenetic marks, recently shown to be affected by polymorphism in a similar way to methylation. A study of these other epigenomic marks in a similar context to that we describe for DNA methylation would be fascinating.

21Despite the evidence for genotypic influence independent of environmental influences, we found that the majority of the VMR-CpGs were best explained by the interaction of genotype and maternal/fetal environment (Fig. 5). The environmental measures selected for analysis included proxy factors such as birth weight and gestational age, maternal smoking and depression, all of which are known to influence a broad range of developmental outcomes. Despite the importance of G × E models in explaining the majority of VMR-CpGs, we failed to find VMR-CpGs that were best explained by environmental conditions independently of genotype. This finding emerged despite the considerable evidence for the effects of environmental factors in pregnancy on both epigenetic states and health outcomes. Our data suggest that genotype exerts a moderating influence on such environmental effects.

22The regression models explained a high proportion of the variance in methylation levels for a minority of the VMR-CpGs. For most VMR-CpGs, the models explained only 10%–20% of the variance (Fig. 5D). One reason why the models might fail to explain the majority of the variance for more VMR-CpGs is that the association between methylation and a specific environmental condition may only be apparent in one genotypic subgroup (Fig. 6). This observation is consistent with the emerging view that genotype can determine the degree of environmentally induced phenotypic plasticity (i.e., so-called “plasticity genes”) and that epigenetic mechanisms serve to maintain environmentally induced phenotypic variation. This is a particularly important finding for studies of either population differences or the influence of environmental factors on phenotypic outcomes. The finding of G × E interaction effects on the epigenome suggests that the necessary level of interrogation extends beyond simple EWAS analysis to include genotype. The failure to include assessment of genotypic moderation of environment–epigenome relations might result in an underestimation of the potential for environmental impact among subpopulations. show that candidate SNP and CpG loci with marginal associations in GWAS and EWAS, respectively, can show strong associations with disease (in this case type 2 diabetes) when combined.

23Another possible reason for the modest explanatory power of certain models is that they are limited to only 19 proxy factors related to the in utero environment and that only single individual genotypic and environmental influences were examined within the models. Models containing multiple SNPs (epistasis) or combining different facets of environments may explain more of the variance in DNA methylation. Moreover, the population under study did not show extreme values in the environmental measures (Table 1), and stronger effects may emerge with the study of high-risk populations.

24Cellular heterogeneity may also compromise our ability to account for variance across the methylome. The umbilical cords that we studied include a mixture of cell types. Recent studies in blood suggest that some of the interindividual variation in methylation is accounted for by differences in cellular content. We were unable to separate the cell types in the GUSTO umbilical cords as the tissues were frozen at collection. Our study involved only a subset of the SNPs in the genome, which included the 1 million SNPs on the OmniExpress + Exome arrays. Although missing SNPs could have accounted for additional variance in the VMRs, the 1 million SNPs were distributed relatively evenly across the genome and represented the majority of haploblocks. Finally, the CpG sites assessed by InfiniumHumanMethylation450 BeadChip are biased toward gene bodies and flanking regions, and therefore we did not cover intergenic regions thoroughly.

25Our analysis was intended to estimate the relative influences of genotypic, environmental, and G × E interactive effects on the neonatal methylome. The models should not be considered as fully determinative of specific outcomes. Indeed, many of the environmental factors considered here are interdependent. Parity and maternal age, for example, are obviously correlated. Moreover, potential sources of environmental influence, such as socioeconomic status, were not considered. Similarly, the associations shown in Figure 6 derive from models that did not control for potentially important covariates. In addition, due to the number of predictive factors tested, it is difficult to ascribe significance to any individual association; instead, we report the models that best explain the variance in methylation values. Nevertheless our data serve to underscore the importance of G × E interactions and suggest that models of epigenetic variation should consider such interactive influences.

26To our knowledge, our report is the first attempt to quantify the relative influence of genotype and environment, as well as their interaction on the human epigenome. This quantification is important as many reports compare DNA methylation to phenotype independently of genotype. Our results strongly suggest that genotype is an essential factor in these relationships. In particular, it is an important question to address in neonates because the influence of prenatal environment on future disease risk is intensely studied with respect to subsequent risk of illness. Our findings suggest that such studies should include an assessment of the degree to which environmental influences are moderated by genotype.

27Briefly, 244 umbilical cord samples from healthy babies who were part of the GUSTO birth cohort study were selected. Subject characteristics can be found in Table 1. Genotyping was performed on the Illumina OmniExpress + Exome array and processed in a standard fashion. DNA methylation profiling was performed on the InfiniumHumanMethylation450 BeadChip. Data were processed as described previously. Sex chromosome data were removed. A temporal batch effect was observed and removed using empirical Bayes methodology (COMBAT). All probes identified as cross-hybridizing in either or were removed from the data set; 301,468 probes remained. (For detailed sample and data acquisition, see Supplemental Material.)

28VMRs were detected as previously described. A candidate VMR was defined as at least two spatially contiguous probes within 1 kb of each other and with MAD values greater than the 95th percentile. We expanded candidate regions to contain more than two probes, as long as the distance between any two neighboring probes within the region was not larger than 1 kb. For each VMR, the CpG with the highest MAD score was utilized as the representative CpG. These 1423 CpGs were used as the data set for subsequent analyses. In addition, the 1500 CpGs with the highest MAD scores outside of VMRs were also captured for further analysis.

29Linear regression was performed for each of the 1423 VMR-CpGs, against all heterologous SNPs identified on the genome-wide arrays. For the purposes of these regressions, the genotype was coded as 1, 2, or 3 and treated as continuous whereby the heterozygote represents an intermediate state between the two homozygotes. In the model analysis below, genotype is categorical. The regressions were adjusted for sex. For each of the 1423 VMR-CpGs, we selected the CpG–SNP pair with the lowest P-value for further analysis. The linear regression analysis was performed in R.

30Methylation (Meth) at each CpG was subjected to the following models [(1), (2), (3)], and models were compared using Akiake information coefficients (AIC).where G1 is the best SNP from previous MethQTL analysis, treated as a categorical variable.where Envi is the phenotype (i = 1–19) (for its corresponding phenotype, see Table 1), which gave the lowest model AIC.where Envi is the phenotype (i = 1–19) (for its corresponding phenotype, see Table 1), which gave the lowest model AIC, and is the SNP on the same chromosome as Meth, which gave the lowest model AIC (also treated as a categorical variable).

31The model with the lowest AIC was declared the “winner,” i.e., the model that best explained the Meth compared with the alternative models.

32Akiake deltas were calculated as the difference between the AIC for the best model and the AIC for the next best model.

33Adjusted R2 was calculated by Equation 4:where is sample R2;, the number of predictors; and, the total sample size.

34For each of the CpGs for which a G × E model had the lowest AIC, the subjects were segregated by the genotype, and regression against phenotype was performed in each genotypic group. Bonferroni correction was performed for 2963 regressions that were run.

35CpG island shores were defined as up to 2-kb regions from the CpG island start or end as per convention. CpG island shelves were next defined as another 2 kb from the shore boundaries, as specified in the GenomeStudio Methylation Module v1.8 User Guide from Illumina (Table 10) (http://supportres.illumina.com/documents/myillumina/90666eaa-0c66-48b4-8199-3be99b2b3ef9/genomestudio_methylation_v1.8_user_guide_11319130_b.pdf). Open seas are regions that are not islands, shores, or shelves. The TSS SwitchGear track in the UCSC Genome Browser was used to delineate human promoters, and genomic coordinates of human enhancers were obtained from the VISTA enhancers track. For individual VMR lists, we determined the total number of probes belonging to each of the six genomic categories (CpG island, south shore, north shore, south shelf, north shelf, and open sea) and also for a background list of 55,003 regions that are possible on the InfiniumHumanMethylation450K BeadChip. The P-value for enrichment of the region lists with each genomic category is computed by a hypergeometric test (one-tailed).

36The data from this study have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE53816 and GSE54445.