Multigenerational epigenetic inheritance in humans: DNA methylation changes associated with maternal exposure to lead can be transmitted to the grandchildren

Contribution of covariates to detectable changes in DNA methylation in dried blood spot (DBS)

We performed two analyses with the DNA methylation data. First, we determined whether high BLLs in the MNBS can significantly affect the methylation profiles of the DNA isolated from the neonatal whole-blood of her children (“Effects of BLL in MNBS on 5 mC in CNBS”; Fig. 1B). Second, we determined whether similar DNA methylation changes occur in the current blood of her children. Using single-value decomposition (SVD) analysis, we explored the contribution of each covariate (maternal BLL, maternal smoking status, gender of offspring, and gestational age) on the DNA methylation profiles of our CNBS samples. SVD does not manipulate the data; rather it calculates the p-values based on the contribution of each component to the variability in the sample. The results indicate that majority of the variation in DNA methylation signal is determined by factors inherent to the HM450K array such as bisulfite conversion rate, hybridization efficiency, and target probe extension (p < 10⁻⁵) (Supplementary Figure 1A). Among biologically- relevant factors, smoking during pregnancy significantly contributed to the variation (p < 0.01) (Supplementary Figure 1A). This is consistent with previous studies that showed that prenatal tobacco smoke exposure affects global and gene-specific DNA methylation, but they used an earlier version of the HM450K assay that had approximately 27,000 probes¹³. However SVD was unable to detect any significant contribution of other factors such as gender, gestational age of the mothers, or age of the infant or sample BLL in our samples.

Studies have shown that gender can significantly impact the DNA methylation status of autosomal genes¹⁴. Therefore we looked at gender-specific differences at single-CpG sites in the DNA methylation profile, while controlling for covariates such as child’s neonatal and current BLL levels, mother’s neonatal BLL levels, gestational age, and age of mothers and children, using a mixed-effect model¹⁵. Similar to what was reported previously¹⁴, for CNBS, we found 179 autosomal CpG sites with gender-dependent DNA methylation effects at an FDR corrected P-value ≤ 0.05 (Supplementary Figure 1B). For CCBS we found almost the same number of CpG sites N = 144 with gender-dependent DNA methylation effects at mostly the same loci (Supplementary Figure 1C). The numbers of differentially methylated CpG sites associated with gender are modest in CNBS and CCBS suggesting that gender specific differences in DNA methylation are probably not established early in development. Studies have also reported that the DNA methylation levels in whole blood can be confounded by variability in the distribution of the blood-cell-type population¹⁶. Houseman et al., 2012 proposed a statistical model for estimation of blood-cell-type proportion using HM450K DNA methylation data¹⁶. We observed a significant variation in estimated granulocyte population across CNBS samples (Supplementary Figure 1D). However, generalized linear model (GLM) analysis showed no significant differences in the blood cell types when the samples were grouped by either the mother’s current smoking status or the mother’s neonatal BLL. This suggests the change in blood cell type population and gender has little contribution toward exposure associated variations in DNA methylation observed in neonates.

Lead (Pb) exposure dependent changes in DNA methylation profile modeled as CpG clusters

Sofer and colleagues suggested that adjacent CpG sites might show similar behavior in response to environmental exposure¹⁷. These sites can be assigned to clusters based on the methylation signature (β value) correlation across multiple samples and the pre-specified-base-pair windows. The effect of exposure on these CpG clusters can be tested using linear models such as the generalized-estimating equation (GEE)¹⁷. Using this powerful strategy for differential methylation analysis, we explored the methylation differences in DNA obtained from CNBS with grandmothers having shown either high BLL (≥5 μg/dl) or low BLL (≤5 μg/dl) during pregnancy (i.e., high or low BLL in the MNBS, assuming that the mother’s neonatal BLL corresponds to the grandmother’s BLL). This analysis was performed while controlling for covariates such as sample BLLs, age of the mother, gender of the child, smoking status and gestational age of the mother. At a significance level (FDR p-value < 0.05 and exposure effect cut-off = 0.02), we found 183 CpG clusters mapping to 564 CpG sites which were differentially methylated in CNBS with high BLLs in MNBS compared to low BLL in MNBS (Supplementary Table 2). We observed more hyper-methylated CpG clusters (n = 151) compared to hypo-methylated CpG clusters (n = 32) at an exposure effect cut-off of 0.02 (i.e., a 2% change in DNA methylation) (Supplementary Figure 2A). Further increasing the cut-off for differential methylation calls to 0.05 or 5% resulted in significant reduction in identified 5 mC clusters (n = 115). However, the trend in DNA methylation patterns was preserved, with a greater number of hyper-methylated regions (n = 98) compared to hypo-methylated regions (n = 17). A subset of these regions mapped to genes and enhancers which showed significant changes in the DNA methylation status from children (CNBS) who have grandmother’s with high BLL. (Table 1 and Supplementary Figure 3A–D). We note that 65% of our clusters had a change in DNA methylation greater than 5%, and over 10% had a change in DNA methylation over 10% (Supplementary Figure 4).

We were also interested in the Pb-exposure-dependent behavior of CpG sites in the children’s more recently obtained blood (CCBS). We hypothesized that the affected CpG sites in the child’s current blood will eventually re-establish the “normal” DNA methylation levels of the low BLL CCBS. In other words, the effects of acute high BLL on DNA methylation will “normalize” over time, as the blood cells undergo hundreds of mitotic divisions in the absence of Pb. To determine whether fetal germ-line exposure induced epigenetic changes “normalize” during the child’s development, we looked for CpG clusters in the CCBS which are associated with mother’s neonatal BLL, while controlling for covariates such as CCBS BLL, age and gender of the infants, smoking status and gestational age of the mothers. As expected, we found very little correlation between mother’s BLL and CCBS DNA methylation (14 CpG clusters mapping to 37 CpG sites at exposure effect size cut-off of 0.02 or 2%) (“Effects of BLL in MNBS on 5 mC in CCBS”; Fig. 1C). This suggests that, during postnatal development, a child’s changes in DNA methylation returns to the “normal” pattern (i.e., low BLL pattern) within the first 3–5 years of life. We also hypothesized that prenatal and post-natal Pb-exposure will have similar effects on DNA methylation at “Pb-sensitive” CpG sites. In a previous study, on the same cohort, we had identified 116 clusters mapping to 320 CpG sites which were sensitive to acute post-natal exposure to Pb in CCBS, independent of gender at a FDR corrected P-value cut-off of 0.05³. We checked for overlap of these sites with the 564 CpG sites which showed statistically significant changes (FDR cut-off 0.05, exposure effect size to in-utero exposure to Pb and found only 6 CpG sites in common (Fig. 1D).

Samples and sample classification

Methods were carried out in accordance with guidelines that were approved by the Wayne State University (WSU) Internal Review Board (IRB), the Michigan Department of Community Health (MDCH) IRB, and the Michigan Neonatal Biobank (MNB) IRB. Experimental protocols were approved by the WSU IRB, the MDCH IRB, and the MNB IRB. Informed consent was obtained from all subjects. For the study, we selected 35 dried blood spots (DBS) collected from mother-infant pairs from Health Fairs ran in three Detroit communities, Rosa Parks, Chene, and Kettering-Butzel, because they have a high prevalence (8–11%) of high BLL in children. The study only included mothers born after January 1, 1987 in Michigan with biological children ages, 3 month to 5 years also born in Michigan. The final sample consisted of 25 male children and 18 female children. We also collected the neonatal DBS and mother neonatal DBS for these mother-infant pairs from the Michigan Neonatal Biobank. Blood Pb measurement was done using Atomic absorption spectroscopy. The corresponding Pb measurement and covariate information can be found in Supplementary Table 1. Methods were carried out in accordance with guidelines that were approved by the Wayne State University Internal Review Board (IRB), the Michigan Department of Community Health (MDCH) IRB, and the Michigan Neonatal Biobank IRB.

Extraction, shearing and denaturation of DNA

DNA was isolated from dried blood spots with Qiagen EZ1 Advanced® using the DNA Investigator® reagents and protocol card. The “Stains on Fabric” preprocessing and Trace® (tip-dance) instrument protocol was used for isolation. The Quantifiler Human DNA Quantification Kit® (Applied Biosystems, Inc.) was used to determine the amount of amplifiable DNA. Approximately 3 μg genomic DNA was diluted in 130 μl of buffer TE (10 mM Tris (pH 8.0), and 1 mM EDTA (pH 8.0)) and sheared into ~200–600 bp fragment using microcavitation (Covaris, Inc, setting: Duty Cycle = 5%, Intensity = 3, Cycles/burst = 200, Time = 75 seconds run at 6–8 °C). 125 μl of the sheared DNA samples were mixed with 330 μl of buffer TE. The sheared DNA was denatured by boiling it in the Thermomixer at 95 °C and 700 rpm for ten minutes and left on ice for 10 min.

HM450K bead chip array

The Infinium HM450K array is a high-throughput technology which assesses the methylation status of approximately 450,000 sites on a sample of human genome^42,43. In the Infinium technique, genomic DNA is sheared, denatured, and bisulfite treated. Bisulfite treatment, if applied appropriately, converts cytosine to uracil, but leaves ^5mC essentially untouched⁴². The DNA is then denatured again and annealed to beads containing DNA oligonucleotides that are complementary to the bisulfite-treated DNA with a 3′ end adjacent to the cytosine whose methylation is being queried. Detection is facilitated by two different probe types (Type 1 and Type 2 probes). The Type 1 probes or the Infinium 1 (Inf1) probes consist of the methylated bead and an un-methylated bead⁴⁴. If the probe for methylated DNA matches the target site there is a single base extension which results in detection which signals into the red channel. Similarly if an un-methylated probe binds to the DNA it signals into the green channel. The type 2 probes or the Infinium 2(Inf2) queries both methylated and unmethylated DNA on a single bead, and the ratio of incorporation of two differently-colored fluorescent nucleotides (Signals A and B) determines the methylation signal. The results are represented in the form of β values, specifically, the average β value (AVG_Beta), representative of the average methylation level of the CpG dinucleotide, and a delta β value which signifies the difference in methylation levels between the control and the experimental group. The Beta or β for the ith interrogated CpG nucleotide is:

where y_i,methy and y_i,unmethy are the intensities measured by the i^th methylated and unmethylated probes, respectively. Illumina recommends adding a constant offset α (by default, α = 100) to the denominator to regularize β value when both methylated and unmethylated probe intensities are low. The Beta-value statistic results in a number between 0 and 1, or 0 and 100%⁴⁵. The raw data was retrieved from Genome Studio methylation module version 1.8™ in the form of 2 files; a sample methylation profile and control probe profile. Quality control, signal correction and normalization of the data was carried out using the HM450K BeadChip data processing pipeline proposed by Teschendorff et al., 2013 in R environment (R > 3.0.0)⁴⁶.Several studies have indicated that the Infinium 1 and 2 probes differed in chemistry, henceforth the HM450K are two separate experiment combined as one. The Infinium 1 probes were shown to have a more stable signal and extended dynamic range compared to the Infinium 2 probes⁴⁷. Therefore, a 3-state beta mixture model is utilized to assign methylation values to specific methylation states. Then the probability of assignment to particular state is divided in quantiles and finally a methylation dependent dilation transformation is performed to preserve sample monotonicity⁴⁶. Prior to analysis the beta values were corrected for Batch effect using Combat function in R and potential snp-containing probes were removed from the analysis (>2.15)⁴⁸.

Statistical analysis

For studying the gender specific effects we CpG association analysis¹⁵, which analyzes DNA methylation data using a mixed effect model at single CpG site level. For determining the Pb-dependent changes in DNA methylation we used adjacent site clustering algorithm, A-clustering to detect sets of correlated CpG sites and then tested the clusters for multivariate response to environmental exposure to Pb using the generalized estimation (GEE) equation approach. The aforementioned approach is efficiently implemented using the R-package Aclust¹⁷. For determining the differentially methylated clusters we used the recommended Aclust parameters; Spearman correlation, for calculating the distance between adjacent sites (dist_(i,j) = 1 corr_(i,j)), average clustering type, which require that mean distance between two sites be at least 0.25, 1000 bp distance restriction for merging of clusters, which ensures that clusters located far away from each other are not merged together based on correlation. The clustering approach is implemented with a 999 bp merge initiation step, which clusters all sites wedged between 2 high correlated sites within 999 bps of each other together, to reduce the complexity of data and the analysis time for the Aclust step. Finally the data was analyzed using a generalized estimation equation approach and filtered for significant DMCs using FDR corrected p-value cutoff = 0.05 and exposure effect size . To determine the genomic locations of the probes belonging to individual DMCs, they were annotated using the publicly available Illumina Human Methylation 450k annotation data in R (>2.15). The target genes mapping to DMCs were individually visualized using UCSC genomic browser. The Delta beta or the beta difference between the median of the beta values for each probe for low BLL samples and high BLL samples were mapped by the chromosomal location of the probes. If A-clustering is an effective technique for differentially methylated identification we hypothesized that the change in methylation status visualized using UCSC genome browser and Integrative genome viewer (IGV)will correspond to the exposure effect (i.e. increase or decrease in methylation) predicted by GEE in the respective regions and might serve as a useful tool for visualization. Single value decomposition (SVD) was implemented using the package ChAMP and estimation of cell distribution was implemented using the package minfi. All statistical analysis was implemented in R. The enhancer regions where downloaded as .bed files from Andersson et al., 2014⁴⁰, and the distance from the middle of the robust enhancer sites to the middle of the cluster was calculated using R package ChIPpeakAnno. Then the clusters where subsetted by Maximum distance <301 and exposure effect size (Supplemental Table 4).

This research was supported by R01 ES012933 and R21 ES021893 to D.M.R., the WSU-NIEHS Center (P30 ES020957), and a Michigan Bloodspot Environmental Epidemiology Project (BLEEP) pilot grant from the Michigan University Research Corridor (URC) to D.M.R.