Quality control and statistical modeling for environmental epigenetics: A study on in utero lead exposure and DNA methylation at birth

Abstract

DNA methylation data assayed using pyrosequencing techniques are increasingly being used in human cohort studies to investigate associations between epigenetic modifications at candidate genes and exposures to environmental toxicants and to examine environmentally-induced epigenetic alterations as a mechanism underlying observed toxicant-health outcome associations. For instance, in utero lead (Pb) exposure is a neurodevelopmental toxicant of global concern that has also been linked to altered growth in human epidemiological cohorts; a potential mechanism of this association is through alteration of DNA methylation (e.g., at growth-related genes). However, because the associations between toxicants and DNA methylation might be weak, using appropriate quality control and statistical methods is important to increase reliability and power of such studies. Using a simulation study, we compared potential approaches to estimate toxicant-DNA methylation associations that varied by how methylation data were analyzed (repeated measures vs. averaging all CpG sites) and by method to adjust for batch effects (batch controls vs. random effects). We demonstrate that correcting for batch effects using plate controls yields unbiased associations, and that explicitly modeling the CpG site-specific variances and correlations among CpG sites increases statistical power. Using the recommended approaches, we examined the association between DNA methylation (in LINE-1 and growth related genes IGF2, H19 and HSD11B2) and 3 biomarkers of Pb exposure (Pb concentrations in umbilical cord blood, maternal tibia, and maternal patella), among mother-infant pairs of the Early Life Exposures in Mexico to Environmental Toxicants (ELEMENT) cohort (n = 247). Those with 10 μg/g higher patella Pb had, on average, 0.61% higher IGF2 methylation (P = 0.05). Sex-specific trends between Pb and DNA methylation (P < 0.1) were observed among girls including a 0.23% increase in HSD11B2 methylation with 10 μg/g higher patella Pb.

Keywords:

• DNA methylation
• environmental exposure
• lead
• pyrosequencing
• quality control
• statistical methods

Abbreviations

ANOVA	=	analysis of variance
DMR	=	differentially methylated region
ELEMENT	=	early life exposures in Mexico to environmental toxicants
GEE	=	generalized estimating equation
GLM	=	general linear model
H19	=	H19, imprinted maternally expressed transcript (non-protein coding)
HSD11B2	=	hydroxysteroid (11-β) dehydrogenase 2
IGF2	=	insulin-like growth factor 2
K-XRF	=	K X-ray fluorescence
LINE-1	=	long interspersed element-1
OLS	=	ordinary linear regression
Pb	=	lead
PCR	=	polymerase chain reaction

Introduction

Environmentally-induced epigenetic alterations at susceptible life-stages, such as in utero development, may persist and impact disease risk later in life. Epigenetic changes have been associated with early-life exposures to chemicals including bisphenol A, genistein, persistent organic pollutants, arsenic, maternal smoking, and cadmium in mice^1-4 and in human birth cohorts.^5-10 Lead (Pb) exposure is of global public health concern, and recent studies suggest that Pb changes DNA methylation patterns among rodents following perinatal exposure^11,12 and is associated with differences in DNA methylation in humans.^13-17 However, the impact of in utero Pb exposure on the epigenome of neonates is understudied, and its impact on specific genes related to health outcomes associated with Pb exposure remains unknown.

DNA methylation is the most widely studied epigenetic modification in human cohorts due in part to its implications for gene regulation and development^18,19 and the multitude of technologies designed to examine DNA methylation. Epigenome-wide technologies interrogate methylation at thousands to millions of CpG sites,²⁰ and quantitative and semi-quantitative technologies such as pyrosequencing,²¹ Sequenom EpiTYPER,²² and methylation-specific PCR ²³ interrogate methylation at specific loci. While epigenome-wide platforms have many advantages, candidate loci approaches are economical and particularly beneficial for studies with hypothesized genes of interest or for validation purposes. As such, candidate loci approaches are well suited for human cohort studies that investigate specific health outcomes, and are preferable for large numbers of samples due to lower cost.

Technologies for epigenome-wide and candidate approaches differ not only in terms of study questions (discovery vs. testing of specific pathways) and cost, but also in data quality management. While statistical and bioinformatics pipelines for managing data from epigenome-wide technologies are constantly being evaluated, quality control practices for candidate loci approaches have been much less systematic. For example, data analysis tools and methods have been developed to address experimental batch effect, biases in array probe design, and dye bias in epigenome-wide data.^24-27 Normalization to internal controls and correction for batch effect are also standard to the data analysis pipeline for molecular biology methods including Real-Time quantitative PCR.^28,29 However, for candidate loci methylation technologies, such data quality control methods are rarely discussed. While in some instances, batch effects are considered and use of rigorous quality control is reported,³⁰ often the data are used as is and quality control consists of 2 known methylated controls (ambiguously called ‘low’ and ‘high’) on each experimental plate.

In environmental epigenetics studies involving humans, the biological significance of results is complicated by limitations including source tissue available for DNA, timing of sample collection for DNA extraction and exposure assessment, and small effect sizes observed between exposures and epigenetic modifications. While the field works to address these issues, there is also a need to improve and standardize quality control and statistical procedures for handling candidate gene methylation data. The reliability of the data and the ability to observe associations with small effect sizes can be improved by implementing appropriate quality control methods and by selecting statistical methods that will maximize statistical power.

Using data from the Early Life Exposures in Mexico to Environmental Toxicants (ELEMENT) birth cohort, this article pursues 2 goals. First, we compared approaches for quality control and data analysis relevant to environmental epigenetic studies utilizing quantitative candidate loci methylation technologies such as pyrosequencing. Models were compared using simulated data that has similar characteristics to the data observed in the ELEMENT cohort. Since data are simulated according to a known toxicant-DNA methylation association, analysis approaches can be objectively evaluated. In particular, we discuss the need for batch effect adjustments and modeling CpG site variance and dependence when studying the relationship between environmental factors and DNA methylation. Second, we evaluate the influence of in utero Pb exposure (measured via 3 biomarkers) on neonate DNA methylation (via bisulfite sequencing candidate regions in umbilical cord blood leukocyte DNA with the pyrosequencing platform). In this cohort, we found modest evidence for associations, often sex-and CpG site-specific, between in utero leukocyte DNA methylation at long interspersed elements 1 (LINE-1), as well as growth-related genes IGF2 and HSD11B2.

Results

We first give descriptive analyses of the ELEMENT study data that motivated our primary objective, followed by the results of the simulation studies that demonstrate best practices for analyzing this type of data, and finally show results of the Pb-methylation associations in the ELEMENT study.

Descriptive statistics of ELEMENT study data

Characteristics of the study population (n = 247 infant-mother pairs) can be found in Table 1. The majority of infants were full term (96.3%, gestational age ≥37 weeks) and male (58.1%). Few mothers smoked as assessed immediately following pregnancy (3.6% at one month post-partum). All subjects have data from at least one Pb biomarker and methylation from one genic region. 243 participants have methylation data for all genic regions, and 202 have Pb measures in all 3 exposure biomarkers.

Table 1. Characteristics and Pb biomarker levels of study population.

		n (%)	Mean ± SD
Infant	Gestational age (weeks)	245	39.26 ± 1.33
	Male	143 (58.1%)
	Female	103 (41.9%)
	Birth weight (kg)	246	3.173 ± 0.39
Maternal	Age at delivery (years)	247	24.49 ± 4.62
	BMI (kg/m^2, pre-pregnancy)	238	25.05 ± 3.66
	Smoking at 1 month postpartum	9 (3.6%)
	Folate (μg intake/day)	247	366 ± 144
Pb Biomarkers	Umbilical cord blood Pb (μg/dL)	227	6.49 ± 3.52*
	Maternal tibia Pb (μg/g)	243	9.82 ± 9.45
	Maternal patella Pb (μg/g)	226	14.4 ± 14.8

* Median umbilical cord blood Pb = 5.80 μg/dL

Lead exposure biomarkers reveal a wide range in cord blood Pb (1.2-29.9 μg/dL), maternal mid-tibia shaft Pb (−20.8-51.0 μg/g of bone), and maternal patella Pb (−29.9-77.7 μg/g). Negative point estimates for bone Pb can result when the true value is close to zero due to normalization to bone density. Negative values passing other quality control checks are retained here as this has been shown to be a better use of this type of data in epidemiological studies.^31,32 Tibia and patella Pb levels quantified at one month postpartum (Pearson r = 0.26, P < 0.001) and patella and log-transformed cord blood Pb (r = 0.19, P = 0.006) are significantly correlated while cord blood Pb and tibia Pb do not correlate strongly (r = 0.11, P = 0.1).

Table 2 shows CpG site-specific and average percent methylation across sites of repetitive elements (LINE-1) and 3 candidate genes quantified via pyrosequencing using bisulfite converted DNA from umbilical cord blood leukocytes. Means and between-subject variance in percent methylation differed across sites from the same genic region. Supplemental Table S1 provides correlations of methylation among CpG sites for LINE-1 and 3 genes. While some pairs of sites are highly correlated with respect to their percent methylation, others are not, and yet others exhibit negative correlations.

Table 2. Percent of methylated cells at four candidate regions. Site-specific methylation and average of 4 CpG sites for repetitive elements, LINE-1, 4 sites within H19 and HSD11B2, and 3 sites for IGF2 are reported. Standardized values were adjusted for batch. LINE-1 was adjusted for batch by subtracting the 0% control value for each CpG site from the experimental plate each sample was run on. IGF2, H19, and HSD11B2 values were adjusted to a standard curve of control DNA (ranging between 0 and 100%) run for each CpG site on the corresponding experimental plate for each sample.

	CpG Site	n	Raw Data Mean ± SD	Standardized Mean ± SD
LINE-1 (%)	Average	246	83.8 ± 1.9	80.8 ± 2.3
	1	246	85.9 ± 2.1	83.4 ± 2.3
	2	246	86.1 ± 1.7	83.0 ± 3.0
	3	246	84.2 ± 2.8	81.0 ± 3.1
	4	246	79.1 ± 2.4	75.8 ± 2.5
IGF2 (%)	Average	243	59.5 ± 5.2	62.5 ± 7.3
	1	243	58.6 ± 6.0	58.0 ± 9.5
	2	243	60.4 ± 5.2	67.4 ± 6.6
	3	243	59.3 ± 5.2	62.2 ± 8.5
H19 (%)	Average	246	64.3 ± 3.4	66.3 ± 4.1
	1	246	66.7 ± 3.9	65.0 ± 4.5
	2	246	63.6 ± 3.4	62.4 ± 3.6
	3	246	63.7 ± 3.4	66.5 ± 5.3
	4	246	63.2 ± 3.5	71.4 ± 5.2
HSD11B2 (%)	Average	247	2.13 ± 1.37	4.20 ± 2.33
	1	247	2.65 ± 1.78	6.42 ± 2.61
	2	246	0.65 ± 1.64	2.40 ± 3.61
	3	236	3.92 ± 2.46	1.30 ± 4.25
	4	228	1.30 ± 1.65	6.69 ± 2.97

ANOVA revealed batch effects in both DNA methylation and lead biomarkers. Significantly different average methylation levels by experimental plate for all 4 candidate regions were found (Table S2). IGF2 batches produced the largest plate-to-plate methylation difference (8.8% difference between plates with highest and lowest averages). Tibia Pb and patella Pb levels also differed when comparing groups of subjects in each pyrosequencing batch (Table S2) with a difference as large as 9.93 μg/g average patella Pb observed between the subjects on 2 LINE-1 plates. Although differences in lead biomarkers were not related to the pyrosequencing procedure, not accounting for these differences could lead to inflated Type I error, thus complicated data analysis.

Comparisons of modeling approaches and batch adjustment

Statistical model comparisons

Using a simulation study, we compare ordinary least squares (OLS) regression (where the average among CpG sites is taken to obtain a single outcome measure per individual) and 2 marginal models for repeated measures, the General Linear Model (GLM) and generalized estimating equations (GEE) (where the actual values at each CpG site are treated as a multivariate outcome) with regard to their performance (bias and power) when estimating the association between a toxicant and DNA methylation.³³ Repeated measures models have utility in DNA methylation studies because, in addition to estimating the toxicant-DNA methylation relationship averaged across sites, they can also be used to rigorously test whether the association varies across CpG sites, and may have higher power than OLS.

Figure S1 compares several approaches to modeling DNA methylation data in the absence of batch effects using simulated data. While all approaches yield consistent estimates of the toxicant-DNA methylation association (not shown), the General Linear Model (GLM) for repeated measures that explicitly models both the site-specific variances and an unstructured correlation matrix among sites yields higher power to detect associations when data have differences in variances across sites similar to those observed in the ELEMENT data. Differences in power are more pronounced for HSD11B2 (Fig. S1B) compared to IGF2 (Fig. S1A) due to more pronounced differences in the variances (Table 2) and pairwise correlations (Table S1) among sites. In the absence of differences in variances across sites, the GLM does not lose power compared to other approaches (not shown). Since differences in variance across sites are common, and pairwise correlations among sites do not necessarily follow a particular pattern, the GLM is preferred over GEE, and much more so than OLS.

Figure 1. Simulation study findings. Bias in estimated coefficients and power to detect significant toxicant-DNA methylation associations from 5 approaches (see Table 3), when data are generated to follow a similar structure to IGF (A and B) and HSD11B2 (C and D) methylation in the ELEMENT study (see Table S3 for simulation settings). For the batch-adjustment method using control samples, the scenario with standard curves having R² = 0.99 (as we had in the ELEMENT study) is shown here. Bias and power are calculated from the average of 1000 data sets, each with 220 subjects.

Batch effect adjustment methods

Two methods to control for batch effect are compared to no batch correction. One approach is to add a random intercept for each experimental plate to the above described models. Another approach is to standardize the data using plate- and site-specific standardization curves constructed using methylation data from control samples with known methylation values. Figure 1 shows the bias and power that can be expected using these methods for batch correction (numbered according to the models described in Table 3); results when using the GLM approach in the ideal case of no batch effects and when using OLS are included as reference. Figure 1 shows that although both approaches have similar power, the approach using only random effects for plate yields biased effect estimates. The extent of the bias when the raw data are used depends on the extent to which the plate- and site-specific standardization curves deviate from the 45 degree line: if the standard reference curves are a straight line with slope other than one, biased associations will be obtained when using raw data.

Table 3. Statistical models tested. Models used to examine associations between Pb biomarker levels and DNA methylation at LINE-1, IGF2, H19, or HSD11B2. Model #4 was considered the best approach though additional modeling types that are sometimes used with candidate gene methylation data were run for comparison purposes.

Model	Batch Adjustment	DNA Methylation Outcome
1. OLS	None	Average of all CpG sites
2. GLM* observed data	None	Individual CpG sites
3. GLM* observed data + RE batch	Random intercept for batch	Individual CpG sites
4. GLM plate standardized data	Site-and plate- specific standard curve	Individual CpG sites
5. OLS plate standardized data	Site-and plate- specific standard curve	Average of all CpG sites

* General linear model for repeated measures using unstructured variance-covariance matrix

Using standard curves to adjust for batch effects will consistently give unbiased associations. However, this approach may suffer loss of power depending on the fit of the standard curves (Fig. S2). If the standard curves fit well, then the standardization process will lead to estimating unbiased associations without loss of power. When appropriate, including a random intercept for batch to Model 4 may help gain power back, while keeping associations unbiased (Fig. S2, Model 4RE).

Associations between in utero lead exposure and DNA methylation in the ELEMENT study

Primary findings

Tables 4 and S4 show the estimated associations between 3 lead exposure biomarkers and DNA methylation using the GLM for repeated measures using plate- and site- specific standardized methylation measures (i.e., model #4 in Table 3) adjusting for gestational age, sex, and maternal folate intake. We found the patella lead biomarker was associated (P = 0.05) with IGF2 hypermethylation (0.6% methylation per 10 μg/g patella Pb) and marginally (P = 0.06) associated with hypermethylation of the HSD11B2 promoter (0.1% methylation per 10 μg/g patella Pb).

Table 4. IGF2 and HSD11B2 models. Estimates and p-values for the regression of methylation at IGF2 or HSD11B2 on each Pb biomarker are noted for a series of models described in Table 3. All models adjusted for sex, gestational age, and maternal folate intake, and models 2–4 also adjust for CpG site. Models included all subjects with available data (n = 223, 237, 221 for IGF2 models with umbilical cord blood, tibia, or patella Pb, respectively; n = 227, 241, 224 for HSD11B2 models with umbilical cord blood, tibia, or patella Pb, respectively).

		Umbilical Cord Blood Pb		Tibia Pb		Patella Pb
Gene	Model	Estimate^# (SE)	P-value	Estimate* (SE)	P-value	Estimate* (SE)	P-value
IGF2	1	−0.33 (0.73)	0.65	−0.54 (0.36)	0.13	0.37 (0.23)	0.12
	2	−0.79 (0.63)	0.21	−0.53 (0.32)	0.10	0.16 (0.21)	0.45
	3	−0.48 (0.55)	0.38	−0.26 (0.28)	0.35	0.11 (0.18)	0.57
	4	−1.36 (0.92)	0.14	−0.00 (0.45)	0.99	0.61 (0.30)	0.05
	5	−1.17 (1.00)	0.24	−0.52 (0.50)	0.30	0.53 (0.34)	0.11
HSD11B2	1	−0.04 (0.19)	0.84	−0.07 (0.09)	0.46	−0.01 (0.06)	0.84
	2	−0.10 (0.17)	0.57	−0.08 (0.08)	0.34	0.01 (0.05)	0.78
	3	−0.10 (0.16)	0.56	−0.08 (0.08)	0.33	0.02 (0.05)	0.65
	4	−0.07 (0.25)	0.76	−0.06 (0.12)	0.61	0.14 (0.08)	0.06
	5	−0.18 (0.32)	0.57	−0.01 (0.16)	0.95	0.22 (0.10)	0.04

^#difference in % methylation per 1 log-unit increase in cord blood Pb

* difference in % methylation per 10 μg Pb/g bone

Site-specific effects

Differences in CpG site-specific relationships with Pb were tested with a type III test of fixed effects for the interaction between Pb and CpG site in a given genic region. Results revealed several significant interactions between Pb and CpG site: IGF2 and tibia Pb, LINE-1 and patella Pb, H19 and patella Pb. Only the interaction between patella Pb and LINE-1 sites remained significant following Bonferonni correction for multiple testing (P ≤ 0.004, Table S5). LINE-1 CpG site #2, but not other sites, exhibits hypermethylation with higher patella Pb levels, a relationship also apparent when separate linear regression models are run for each CpG site (Table S5).

Sex-specific effects

To explore sex-specific relationships, an interaction term was added between Pb and sex in Model #4. Results are reported for models with a Pb-methylation association with P < 0.1 in at least one sex following stratification (Supplemental Table S6). Three suggestive relationships between Pb biomarker levels and methylation were observed solely among females. For example, a natural-log unit higher cord blood Pb among females is associated with 3% lower methylation in IGF2. This decrease (P = 0.04) would not withstand correction for multiple hypothesis testing and is not observed among the males (Fig. 2). HSD11B2 hypermethylation was noted with higher patella Pb levels among both males and females, though effect size was greater and only marginally significant among females (Fig. 3, Tables 4 and S6).

Figure 2. Umbilical cord blood Pb and IGF2 methylation stratified by sex. Percent methylation in umbilical cord blood leukocyte DNA of the imprinted gene, IGF2, decreased with cord blood Pb solely in girls. After adjusting for gestational age and maternal folate intake, one natural-log unit increase in cord blood Pb was associated with a 3% decrease in IGF2 methylation among girls (P = 0.04) and the decrease among boys was not significant (P = 0.66). IGF2 data were first adjusted by a standard curve of known methylated or unmethylated controls from each experimental batch.

Figure 3. Maternal patella Pb and HSD11B2 methylation stratified by sex. Percent methylation upstream of the HSD11B2 promoter in blood leukocyte DNA increased with maternal patella Pb levels, indicative of in utero exposure. In a repeated measures model adjusting for gestational age and maternal folate intake, a 10 μg/g increase in maternal patella Pb of the girls was associated with a 0.23% increase in HSD11B2 methylation (near significant, P = 0.07). The increase was less notable among boys (P = 0.41). All methylation values were first standardized to known control curves from the same experimental batch.

Sensitivity analyses

Relationships between Pb and methylation remained consistent (direction, attainment of statistical significance) between continuous and categorical bone Pb models (data not shown). Furthermore, exclusion of one patella Pb outlier in the continuous Pb models did not substantially alter reported results; in IGF2 and HSD11B2 models, statistical significance increased when excluding the outlier. Because it is unknown if true methylation differences due to exposures are in an absolute scale or if the exposure effects are proportional to the variance at the CpG site, we created Z-scores for each CpG site and estimated Pb-methylation relationships using Model #4 with the Pb-by-site interaction term. Similar relationships were observed (see one example in Table S5).

Modeling ELEMENT data using other approaches

As an example of the difference in inferences that can be obtained with other common statistical approaches using real data, we also fitted the remaining models described in Table 3. In all cases, the models were adjusted for gestational age, sex, and maternal folate intake (see Tables 4 and S4 for results). Substantial changes to direction of effect and statistical significance occur when batch is accounted for via standardized methylation values. For example, the relationship between hypomethylation of LINE-1 and patella Pb apparent with Model #1 disappears (Table S4, Model 1 vs. Model 5 both using OLS). In contrast, an association between HSD11B2 hypermethylation with higher patella Pb emerges (Table 4, Model 1 vs. either Model 4 or Model 5 depicted in Fig. 3). Similarly, substantial changes to effect size and significance level were noted for associations between patella Pb and both IGF2 and HSD11B2 methylation (Table 4) when comparing to no batch adjustment (Model #2). Controlling for batch via a random intercept (Model #3) in many cases diminished the estimated association (Tables 4 and S4) as can be expected given the simulation study findings.

Discussion

This study highlights the need to carefully consider statistical modeling and batch adjustment approaches in environmental epigenetic studies with candidate region methylation data obtained using platforms such as pyrosequencing. Based on results from both real data examples and simulation studies that objectively compared modeling and batch adjustment approaches, we draw several recommendations that are particularly geared toward the environmental epigenetics field and other studies expecting small effect sizes for the variables of interest on DNA methylation. We recommend: 1) carefully checking for and correcting for batch effects; 2) utilizing models for repeated measures that incorporate the actual values at each CpG site in the interrogated region and, within those models, specifically modeling site-specific means and variances in methylation; and 3) evaluating site-specific and sex-specific effects.

In the ELEMENT study data, we observed modest evidence for an association between in utero Pb exposure and changes to umbilical cord blood leukocyte DNA methylation. IGF2 and HSD11B2 methylation was higher among those with higher maternal patella Pb levels, and several Pb-methylation relationships were observed among newborn girls but not boys. Although none of these relationships remain significant after adjusting for multiple hypothesis testing (P ≤ 0.004) except a significant interaction between LINE-1 CpG site and patella Pb driven by the second CpG site, these findings point out the importance of exploring site- and sex-specific effects to better understand potential biological mechanisms. Alterations to DNA methylation, specifically of growth-related genes, from in utero Pb exposure (for which patella Pb serves as a biomarker³⁴) have the potential to impact childhood growth and development. In ELEMENT and other cohort studies, in utero Pb exposure has been shown to impact size at birth,^35-37 growth trajectories and size throughout infancy³⁸, and childhood³⁹ in a sex-dependent fashion.³⁹ Further study of the impact of Pb exposure on the DNA methylome is warranted as previous perinatal rodent studies^11,12 and adult cohort studies^13,14,16,17 found links between Pb exposure and percent DNA methylation at various lifestages, and only 3 genes and LINE-1 repetitive elements were interrogated in this study.

Experimental batch effects are known to be problematic and are now routinely corrected for in molecular biology techniques including gene expression analysis with Real Time quantitative-PCR^28,29 and array-based platforms for interrogating DNA methylation across the epigenome.^24,26 By contrast, quality control measures vary widely for pyrosequencing and other high-throughput platforms for quantifying DNA methylation at candidate regions. Many laboratories have developed stringent validation procedures for new pyrosequencing assays that examine full standard curves from commercially available unmethylated and methylated standards or various homemade standards (e.g., plasmid control DNA,⁴⁰ multiple displacement based whole genome amplification with or without SssI-methylase treatment⁴¹). While some high-throughput studies standardize sample methylation levels to controls included in each batch,³⁰ this is not a widespread practice. In multiple batch studies, we highly recommend statistically testing for batch effects in data from each candidate region before proceeding with further data analysis to avoid reporting spurious results. Certain regions may be more prone to batch effects in bisulfite sequencing than others. Batch effects likely stem from the PCR step as biased amplification of bisulfite converted DNA can occur when the melting temperature differs by methylation status of the template.⁴² Repetitive elements such as LINE-1 represent a special case in which a consensus region of thousands of LINE-1 elements are amplified, and the proportion of specific LINE-1s amplified each time can differ leading to slight variations in quantified methylation.

Accounting for batch effects is important in environmental epigenetics studies where the magnitude of change from exposure is expected to be small (<0 .5 SD as we observed) or is similar to the batch effect size. Adjusting for batch effects will often remove the potential for bias in associations and increase power to detect associations since between-batch variation is removed. Adding batch as a random effect when evaluating associations between variables of interest (e.g., chemical exposures) and DNA methylation is one potential correction method (see Model #3, Table 3). However, in the ELEMENT data we discovered that exposure levels happened to vary by group of subjects included in the experimental batches (Table S2). Samples were placed in experimental batches (PCR plates) ordered by subject ID, a common practice in high-throughput analyses. Batches of subjects had different Pb levels which may have correlated with recruitment factors (e.g., timing, region). Randomizing subjects by exposure level on experimental plates is one practice that would eliminate this problem. Between-batch differences in Pb levels prompted us to consider and evaluate other approaches to correct for batch effects.

We found through simulation studies that standardizing to a control curve has comparable power to including random effects for batch, and in addition removes attenuation bias in the estimated effect size that remains when adjusting for batch using random effects. In the emerging environmental epigenetic literature, the practical importance of small effect sizes may be called into question since detected associations may be small. Not standardizing to batch standard-curves will yield attenuated associations when standard curves have slopes less than 1.

Platforms such as pyrosequencing and Sequenom EpiTYPER can quantify methylation at several CpG sites within the same region. Typically, percent methylation at neighboring CpG sites is highly correlated, though individual CpG sites can be differentially methylated and/or have roles in gene regulation (due to location in a transcription factor binding site, for example).^43,44 Further, between-subject variation may differ across CpG sites. Statistical analysis often involves averaging methylation across all quantified sites (see Models #1 and 5, Table 3). Along with others in the field, we recommend further evaluation of CpG site-specific effects by treating methylation measures at each site as a repeated measure using a model that explicitly models the variances and correlations across sites (Models #2–4).^1,5,45-47 This approach yields higher power to detect associations compared with the OLS approach that averages site methylation data (Model #5). Incorporating the correct variance-covariance structure in repeated measures models is known to increase the precision of effect.³³ Utilizing an unstructured covariance matrix allowed us to best approximate the true relationships between CpG sites in an interrogated region which differed widely by genic region (Table S1), and did not incur loss of power due to estimating 6–10 variance-covariance parameters. Difference in within-gene CpG site variances and correlations across assayed regions has also been observed in other studies.⁴⁸

An additional advantage of the repeated measures model is the ability to test whether specific CpG sites have distinct relationships with exposure by adding an interaction term between site and exposure to the model (Table S5), a method others have found useful in determining labile CpG sites.¹ Upon performing this analysis, a relationship between patella Pb and hypermethylation of LINE-1 CpG site #2 was identified (Table S5). Although this relationship could also be identified in individual models run for each LINE-1 CpG site, a formal assessment of whether the associations differ across sites cannot be readily computed. The repeated measures model has the advantage of estimating and testing differences in site-specific relationships within a single model.

DNA methylation patterns differ between the sexes.^5,49,50 Sex-specific relationships between toxicant exposures and DNA methylation have been observed in animals¹¹ and humans.^7,9,51 In this study, LINE-1 was hypomethylated in girls compared to boys regardless of exposure (ANOVA P = 0.04, data not shown), a pattern previously observed in the Center for the Health Assessment of Mothers and Children of Salinas (CHAMACOS) cohort⁵ and among adults.⁵⁰ Results following sex stratification suggested that DNA methylation in girls at LINE-1, IGF2, and HSD11B2 may be more sensitive to change following in utero Pb exposure, though sample size was small and significance marginal (Table S6, Figs. 2 and 3). Including sex stratification or exposure-sex interaction terms in statistical procedures is crucial when evaluating relationships between DNA methylation and exposure or health outcomes.

While many environmental epigenetic studies report associations between percent DNA methylation at specific genes and exposures, the effect size is often small in both epidemiological birth cohort studies^6-8,52 and adult cohort studies.^16,17 For example, given the estimated 0.14% difference in HSD11B2 methylation per 10 μg/g higher patella lead, we would expect those with the 75th percentile patella Pb concentration to differ by 0.27% methylation (or 0.45% methylation among girls) compared to those with the 25th percentile of patella Pb. The biological significance of such a small change on gene expression, recruitment of binding proteins, other epigenetic modifications, or downstream health outcomes is uncertain. Even so, this difference of 0.27% is equivalent to 0.2SD of HSD11B2, and an entire range increase (minimum to maximum) is expected to produce a 0.9 SD increase in HSD11B2 methylation among girls. Certain regions of the genome exhibit greater methylation variability and may be responsive to environmental influences, especially during early in utero development. If environmental factors such as Pb influenced >0.5 SD of percent methylation in these regions, termed metastable epialleles,^53-55 the change could be biologically significant.

Environmental epigenetics has emerged as a field with immense promise to identify 1) early life exposures causing persistent changes to the epigenome that may impact disease susceptibility,⁵⁶ 2) epigenetic biomarkers of past exposures that cannot be measured,^6,17 or 3) reversible epigenetic modifications to target therapeutically or nutritionally.⁵⁷ Thus, despite limitations and uncertainty in the current state of environmental epigenetics research, important links have been made between early life exposures and epigenetic modification in animal models ^1-4 and humans.^5-10 Incorporation of stringent quality control measures in the laboratory and careful statistical analysis of data following recommendations described in this article will increase confidence in results obtained from quantitative bisulfite sequencing of candidate regions. This will enable researchers to better evaluate associations between exposures and DNA methylation at key genes. Coupling methylation data with gene expression and health outcome data (statistically via meditational analyses) will further improve our understanding of the biological implications of observed changes.

Materials and methods

ELEMENT study

Sample population

In 1994 and 1995, mother-infant pairs were recruited at delivery at 3 Mexico City hospitals as the first cohort of the Early Life Exposures in Mexico to Environmental Toxicants (ELEMENT) study. ELEMENT Cohort I was originally designed to test the ability of calcium supplementation to decrease maternal to infant Pb mobilization during lactation. Exclusion criteria included living outside the metropolitan area, multiple fetuses, and maternal conditions that could interfere with in utero development (e.g., preeclampsia, gestational diabetes, seizure disorders, psychiatric, kidney or cardiac diseases). Full exclusion criteria and cohort details are described elsewhere.^36,38

Umbilical cord blood samples were collected from the neonate. Maternal and neonate anthropometry was recorded within 12 hours of delivery by trained obstetric nurses as previously detailed.^35,36 Mother-infant pairs visited the research center at one month postpartum (±5 days) for further evaluation and maternal bone Pb measurement. Of 617 mother-infant pairs who participated in the calcium supplementation trial, 248 umbilical cord blood samples were available for DNA extraction for the current investigation. Of these, high quality DNA was extracted from 247 infants, and these subjects are included in the analyses described herein.

Study procedures and strategies for Pb exposure reduction were explained to participating mothers, and mothers provided written consent. The research protocol was approved by the Human Subjects Committee of the National Institute of Public Health of Mexico, participant hospitals, and the Internal Review Board at all participating institutions including the University of Michigan.

Questionnaires

Mothers were interviewed at one-month postpartum to obtain information on demographic characteristics, health status of the mother and infant, smoking, and reproductive history. Gestational age was extracted from the medical record. A semi-quantitative food frequency questionnaire was administered to estimate maternal intake of nutrients during pregnancy. This questionnaire was previously validated in women living in Mexico City.⁵⁸

Lead (Pb) measurements

In utero Pb exposure was approximated via Pb levels in 3 biomarkers: umbilical cord blood, maternal tibia, and maternal patella using methods previously described.^36,38 In brief, blood Pb was quantified via atomic absorption spectroscopy with a Perkin-Elmer 3000 (Chelmsford, MA) at the American British Cowdray Hospital in Mexico City. Good precision and accuracy was obtained when comparing measurements with blinded replicates performed at the Wisconsin State Laboratory of Hygiene for quality control purposes.³⁶ Cortical (left mid-tibia shaft) and trabecular (left patella) bone Pb measurements were obtained from mothers at one month post-partum using a spot-source ¹⁰⁹Cd K-XRF instrument. Details on the physical principles, technical specifications, and validation of the method were previously described.^31,59,60 The Pb fluorescent signal is normalized to an indicator of bone density - the scattered γ-ray signal stemming from calcium and phosphorus in bone mineral – which results in negative values in some instances. Each K-XRF measurement also comes with an estimate of the uncertainty (e.g., measurement error which is increased by factors such as low bone density and obesity). Here, bone Pb measurements were excluded if the uncertainty estimate was >15 for patella Pb or >10 for tibia Pb.

Epigenetic analyses

DNA was isolated from umbilical cord blood leukocytes by the Harvard-Partners Center for Genetics and Genomics using PureGene kits (Gentra Systems) and stored at −80°C before transfer to the University of Michigan on dry ice.⁶¹ DNA concentration was quantified with a NanoDrop 2000 (ThermoFisher Scientific). Genomic DNA (0.5–1 μg) was bisulfite converted according to standard methods with EpiTect Bisulfite Kits (Qiagen). This treatment converts unmethylated cytosine to uracil (and ultimately thymine following PCR) while leaving methylated cytosine intact.⁶²

DNA methylation (percent of methylated cells) was quantified via pyrosequencing at global repetitive elements (LINE-1) and 3 genes important for early growth (IGF2, H19, HSD11B2).^63,64 Methylation at 4 CpG sites in LINE-1 was quantified via our previously published assay.⁴⁶ Differentially methylated regions (DMRs) of the imprinted genes, IGF2 (upstream of imprinted promoter at IGF2 exon 3) and H19 (upstream of H19 within the imprint control center), were interrogated based on validated assays.⁶⁵ PyroMark Assay Design Software version 2.0 was used to design the pyrosequencing assay for HSD11B2 to include 4 CpG sites in a CpG island upstream of the gene promoter (region described by⁶⁶). For each region, the sequence was amplified from approximately 50 ng bisulfite converted DNA by HotStartTaq Master Mix (Qiagen) and forward and reverse-biotinylated primers. Percent of methylated cells at each region was quantified by the PyroMark MD Pyrosequencer Platform (Qiagen) and computed by Pyro Q-CpG Software. The software incorporates controls to check for completed bisulfite conversion and adequate signal over background noise. All samples were run in duplicate, and duplicate reads were averaged. Each assay was validated in the lab by running 6 point standard curves made from serial dilutions of 0 and 100% methylated human control DNA (Qiagen). Every experimental plate (batch) contained at least 4 point standard curves, except LINE-1 plates which only included 0 and 100% methylated control DNA. Measures of methylated control DNA were precise across 96-well plates (e.g., <10% coefficient of variation, CV, for 100% controls). Duplicate samples had good precision for IGF2 (mean 1.8% CV), H19 (1.8% CV), and LINE-1 (1.4% CV). Methylation of HSD11B2 was consistently low (<10%) and variable (36% CV among duplicate samples).

Simulation study to evaluate modeling and batch adjustment approaches

Modeling strategies

Ordinary least squares (OLS) regression is commonly used to estimate relationships among DNA methylation and covariates; in this approach data from multiple CpG sites for a given locus are averaged to obtain a single outcome measure per individual. This approach ignores differences in mean, variance, and pair-wise correlations among CpG sites. Instead, models for multivariate data (i.e., multiple outcome measures per individual) can be used, to directly use the actual values at each CpG site.⁶⁶ Furthermore, models for repeated measures can be used to rigorously test if the toxicant-DNA methylation association varies across CpG sites. While separate OLS models can be fit to each site, a test of whether the associations differ by site cannot be readily obtained.

Models for repeated measures include random effects (RE) models or marginal models. The general linear model (GLM) for repeated measures and generalized estimating equations (GEE) are 2 marginal models that, in the case of continuous outcomes, primarily differ in how variances across sites are modeled: widely available software implementations of GEE typically assume variances are equal across sites, whereas implementations of GLM explicitly require a model for the variances as well as the correlations. A random effects model with a random intercept is equivalent to a GLM with a compound symmetry assumption, and gives the same estimated association as GEE with exchangeable correlation structure. Although robust standard error estimators are available to guard against incorrect modeling of the variance-covariance structure among repeated measures in the GEE approach, it is well known that mis-specified variance-covariance have lower power.³³ We used simulated data where toxicant-DNA methylation associations are known to compare these modeling approaches.

Batch adjustment

The most common method for batch adjustment is to consider batch effects as randomly shifting the means of the methylation values by a similar amount for all samples in a plate and to model the data including a random intercept for plate in mixed effects models. However, this adjustment method may not be appropriate because the batch effects can potentially vary by CpG site, not only in the mean but also by shifting the scale of the measure (see intercepts and slopes in Table S3). An alternative is to correct for batch effects by standardizing methylation data at each CpG site to a standard curve for that site constructed using measures for known unmethylated and methylated control DNA run on each experimental plate. In the ELEMENT data, the standard curves for H19, HSD11B2, and IGF2 were linear with slopes less than one; i.e., observed methylation for control sample = â + ŝ*known value of control sample, with ŝ < 1 being the estimated slope of the plate- and site-specific curve, and â the estimated intercept, Table S3). For the subjects in the study, the standardized values are calculated as (observed sample- â)/ŝ.

Simulated data

To objectively evaluate the modeling strategies and batch adjustment approaches, we simulated data that were similar to the ELEMENT data, but where the association between a toxicant and DNA values were known. First, data were simulated without batch effects (called “true data”) using Y_ij,true = β_0j + β₁ X_i + e_ij, where β_0j is the mean of site j (See Table 2), X_i had a normal distribution with mean 14.4 and SD = 14.8 (similar to patella bone lead), and e_ij had correlation structure similar to observed in the ELEMENT data (Table S1). The standard deviations of e_ij were set equal to those observed in ELEMENT (Table 2) or fixed to be equal for all sites (as the average of the standard deviations across the sites). The value of the association, β₁, was set to vary from 0 to 0.15; within this range, the power to detect the known associations ranges from 0.05 (Type I error) to 1.

“Observed data” with batch effects were simulated as Y_ij,obs = a_kj + s_kj * Y_ij,true, where a_kj and s_kj were random batch effects specific to site j and plate k generated from a normal distribution mean and standard deviation similar to what was observed in the ELEMENT data (see Table S3). To mimic the fact that in a real situation the intercept and slopes of the standard curve are estimated, we generated standardized data as Y_ijk,std =(Y_ij,obs − â_kj)/ŝ_kj, where â_kj = a_kj + δ_kja and ŝ_k = s_k + δ_kjs include estimation errors, δ_kja and δ_kjs; if we used the generated a_kj and s_kj we would recover exactly the “true data” Y_ij,true. Following simple linear regression formulae, the variance of the estimation errors, δ_kja and δ_kjs, were chosen such that the R² of the fitted standard curves were, on average R² = 0.999, 0.99 and 0.88. For each combination of parameters (β₁, R², variance and correlation structure of e_ij) 1000 datasets, each with 220 subjects, were generated.

Models fitted and model performance measures

We fitted the models described under “Modeling strategies” to show the performance of the various modeling approaches. Average bias (i.e., difference between estimated and known β₁ across the 1000 simulated data sets) and empirical power (the proportion of times the hypothesis H₀: β₁ = 0 was rejected in the 1000 simulated datasets) were computed. The GLM model with unstructured variance-covariance matrix was shown to have the best performance (high power and unbiased estimates). Next, we fitted the GLM to Y_ij,obs, Y_ij,std, and also the GLM with a random intercept for batch to examine impact of the adjustment method on bias and power. For comparison to what is typically done, we also fitted OLS to person-specific average (averaged across site).

Statistical analyses of ELEMENT study data

All statistical analyses were performed with SAS 9.3 (SAS Institute). Univariate statistics and distributions of variables were examined. Cord blood Pb was natural log-transformed in all analyses. Lead exposure biomarkers and methylation values for each candidate gene were compared across experimental batches (one batch = same PCR plate and same pyrosequencing date) by ANOVA (used Welch's ANOVA if variances were not homogenous). Bivariate and multivariable analyses were conducted to examine the association between biomarkers of in utero Pb exposure (cord blood Pb, maternal tibia Pb, maternal patella Pb) and DNA methylation levels at 4 candidate regions (LINE-1, IGF2, H19, and HSD11B2). Given the results of the simulation study, Model 4 (Table 3) was used as the primary modeling approach. Results were considered statistically significant at the P < 0.05 level, but also discussed in the context of multiple hypothesis testing (12 Pb-gene relationships) with a P-value ≤0.004. Multivariable analyses using different modeling strategies and methods for batch adjustments (Table 3) were also employed to compare the methods in a real data example.

Covariates

All models included child's sex, gestational age, and estimated maternal folate intake during pregnancy, which are biologically relevant predictors of methylation. Inclusion of additional covariates (indicators of maternal socioeconomic status, maternal smoking history) selected using backward elimination for specific genes did not substantially change reported results.

Site- and sex-specific effects

Model #4 was run with an interaction term between CpG site and Pb to assess whether individual CpG sites had different relationships with Pb. To examine sex-specific effects which are of increasing importance in public health research, Model #4 was also run with an interaction term between sex and Pb.⁶⁷ If the interaction term and/or main effect of Pb had P < 0.1, the model was run following sex stratification.

Sensitivity analyses

Sensitivity analyses revealed 3 highly influential data points for LINE-1, IGF2, and H19 with percent methylation outside of the expected range. These data points were removed in all reported analyses. Model #4 was also run after removing one patella Pb outlier. Further sensitivity analyses involved using quartiles of tibia and patella Pb instead of the continuous variable to remove any potential influence of the negative bone Pb measures (which were classified as the lowest quartile). Since range and variance of methylation values within genes differed across CpG sites, and it is unknown if toxicants influence methylation in the absolute scale (e.g., equal % difference across all sites), or in a scale relative to standard deviation (e.g., effect is in units of % difference/standard deviation at the CpG site), a separate analysis was performed (Model #4 with and without Pb by CpG site interaction term) using z-scores for each CpG site (i.e., data were converted to have mean 0, SD 1 for each site).

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Acknowledgments

The authors acknowledge the research staff at participating hospitals. We thank the mothers and children for participating in the study. The contents of this publication are solely the responsibility of the grantee and do not necessarily represent the official views of the US EPA or the NIH. Further, the US EPA does not endorse the purchase of any commercial products or services mentioned in the publication.

Funding

This publication was made possible by U.S. Environmental Protection Agency (US EPA) grants RD834800 and RD83543601 and National Institute for Environmental Health Sciences (NIEHS) grants P20 ES018171, P01 ES02284401, R01 ES007821, R01 ES014930, R01 ES013744, P42 ES05947, and P30 ES017885. JMG was also supported by the National Center for Advancing Translational Sciences (grant no. 2UL1TR00433).

Supplemental Material

Supplemental data for this article can be accessed on thepublisher's website.