Effects of Cadmium Exposure on DNA Methylation at Imprinting Control Regions and Genome-Wide in Mothers and Newborn Children.

BACKGROUND: Imprinted genes are defined by their preferential expression from one of the two parental alleles. This unique mode of gene expression is dependent on allele-specific DNA methylation profiles established at regulatory sequences called imprinting control regions (ICRs). These loci have been used as biosensors to study how environmental exposures affect methylation and transcription. However, a critical unanswered question is whether they are more, less, or equally sensitive to environmental stressors as the rest of the genome.
OBJECTIVES: Using cadmium exposure in humans as a model, we aimed to determine the relative sensitivity of ICRs to perturbation of methylation compared to similar, nonimprinted loci in the genome.
METHODS: We assayed DNA methylation genome-wide using bisulfite sequencing of 19 newborn cord blood and 20 maternal blood samples selected on the basis of maternal blood cadmium levels. Differentially methylated regions (DMRs) associated with cadmium exposure were identified.
RESULTS: In newborn cord blood and maternal blood, 641 and 1,945 cadmium-associated DMRs were identified, respectively. DMRs were more common at the 15 maternally methylated ICRs than at similar nonimprinted loci in newborn cord blood (p=5.64×10-8) and maternal blood (p=6.22×10-14), suggesting a higher sensitivity for ICRs to cadmium. Genome-wide, Enrichr analysis indicated that the top three functional categories for genes that overlapped DMRs in maternal blood were body mass index (BMI) (p=2.0×10-5), blood pressure (p=3.8×10-5), and body weight (p=0.0014). In newborn cord blood, the top three functional categories were BMI, atrial fibrillation, and hypertension, although associations were not significant after correction for multiple testing (p=0.098). These findings suggest that epigenetic changes may contribute to the etiology of cadmium-associated diseases.
CONCLUSIONS: We analyzed cord blood and maternal blood DNA methylation profiles genome-wide at nucleotide resolution in individuals selected for high and low blood cadmium levels in the first trimester. Our findings suggest that ICRs may be hot spots for perturbation by cadmium, motivating further study of these loci to investigate potential mechanisms of cadmium action. https://doi.org/10.1289/EHP2085.

Introduction

Modulation of the epigenome has been proposed as a key mechanism through which exposure to environmental stressors during development can program altered gene expression to impact health in later life. Developmental exposure to malnourishment (Lee 2015), toxic metals (Castillo et al. 2012), and endocrine disruptors (Susiarjo et al. 2013) impacts metabolic health, and is associated with alterations to the epigenome in early life, although finding loci that maintain these epigenetic changes into adulthood has proved challenging. Identifying such loci might provide important insights into the mechanisms of this long-term programming effect.

Imprinted genes are expressed preferentially from one of their two parental alleles, dependent on allele-specific epigenetic profiles established during gametogenesis at imprinting control regions (germline ICRs) (Bartolomei and Ferguson-Smith 2011). These sites of differential methylation between the maternal and paternal alleles are protected from the wave of global demethylation that occurs immediately after fertilization, providing a parent-of-origin-specific epigenetic signature. Given that epigenetic profiles at germline ICRs are generally stable throughout life, these loci have been considered candidates for providing a long-term memory of early-life exposures. Additionally, dysregulation of imprinted gene expression can have potent effects on growth and metabolism, suggesting that changes to their epigenetic regulation could directly contribute to effects on metabolic health associated with developmental exposures. Despite the widespread use of imprinted genes as biosensors of environmentally induced epigenetic changes, it has not been determined whether ICRs are more or less sensitive to perturbation than other loci in the genome. Resolving this would provide insight into the mechanisms through which exposures interact with the epigenome.

In the current study, we address this question in humans in the context of cadmium (Cd) exposure. Cd is a toxic metal that is poorly metabolized, accumulating in the liver and kidneys. The primary routes of exposure are inhalation, typically from cigarette smoke and contaminated dust, and ingestion, particularly of cereals and vegetables that absorb Cd from contaminated soil (Järup and Åkesson 2009). Low-level chronic exposure in adulthood has been associated with metabolic disorders, including nonalcoholic fatty liver disease and nonalcoholic steatohepatitis (Hyder et al. 2013), cardiovascular disease (Agarwal et al. 2011), and reduced bone mineral density, causing an increased incidence of bone fractures and osteoporosis (Järup and Åkesson 2009). Postnatal exposure to Cd has also been associated with adverse neurodevelopmental outcomes (Ciesielski et al. 2012) and psychiatric disorders (Orisakwe 2014).

In utero exposure to Cd causes intrauterine growth restriction (IUGR) (Xu et al. 2016) and reduced birth weight in rodents (Ba et al. 2016), and can program an altered metabolic state in later life, including increased systemic fat accumulation (Ba et al. 2016) and the mis-regulation of genes involved in hepatic lipid metabolism (Castillo et al. 2012). Mechanisms linking developmental Cd exposure to pathophysiological outcomes are not well defined, but studies of cultured cells (Gadhia et al. 2015; Xiao et al. 2015a), animal models (Castillo et al. 2012), and human cohorts (Kippler et al. 2013; Sanders et al. 2014) suggest a role for epigenetics. In utero Cd exposure in humans has been associated with differences in methylation at genes involved in hepatic lipid metabolism and bone mineralization, among other functions (Kippler et al. 2013; Sanders et al. 2014), supporting the idea that epigenetic modifications could contribute to the programming of adult disease (Heijmans et al. 2009). However, previous studies have used array-based approaches to assay methylation, limiting the analysis to a small proportion of the genome, potentially missing sites of differential methylation that could contribute to disease pathogenesis. Additionally, Cd-associated epigenetic changes have been studied at only a few ICRs, despite evidence suggesting that exposure to other toxic metals, such as lead, may perturb methylation at multiple ICRs, including at the MEST, IGF2, DLK1/MEG3, and PLAGL1 loci (Li et al. 2016; Nye et al. 2016).

To estimate associations between Cd exposure and DNA methylation in an unbiased and global manner, we performed whole-genome bisulfite sequencing (WGBS) on newborn cord blood DNA from individuals exposed to low and high levels of Cd in utero, selected on the basis of maternal blood Cd levels in the first trimester of pregnancy. We also performed WGBS on the paired maternal blood samples. We compared ICRs to similar, nonimprinted sequences in the genome to determine whether they appear to be more, less, or equally sensitive to Cd-associated methylation changes.

Materials and Methods

Recruitment of Study Participants

Study participants were recruited between January 2009 and June 2011 as part of the ongoing Newborn Epigenetics STudy, a prospective cohort study of women and their children. Details of participant identification, enrollment, and data and specimen collection procedures have been described elsewhere (Hoyo et al. 2014; Liu et al. 2012). Briefly, English- or Spanish-speaking pregnant women aged 18 y or older were recruited from prenatal clinics at Duke and Durham Regional Hospitals, North Carolina, United States. From these, women who planned to relinquish custody of the child or had established HIV infection were excluded, as were those who received obstetric care at facilities other than Duke or Durham Regional Hospitals, and those who were illiterate and could not reasonably provide informed consent. By December of 2011, 1,700 (67%) of the 2,548 eligible women had been enrolled. Of the 1,700 enrolled, 396 were excluded, including 115 due to fetal demise and 281 due to refusal of further participation or inability to locate the participant for subsequent follow-up, leaving 1,304 (77% of those enrolled or 51% of those approached). The median gestational age at enrollment was 12 wk (interquartile range: 7–14 wk). The protocol was approved by Duke University, Durham Regional Hospital, and North Carolina State University. All final participants in the study provided informed consent.

At enrollment, of peripheral blood was obtained in ethylenediaminetetraacetic acid (EDTA) vacutainers from all women, of which was stored as whole blood in aliquots at . For the first 310 participants, one aliquot was used to measure Cd using established solution-based inductively coupled plasma mass spectrometry (ICP-MS) methods that have been previously described (Darrah et al. 2009; McLaughlin et al. 2011; Vidal et al. 2015). Briefly, the frozen samples were equilibrated at room temperature, homogenized, and then pipetted into a trace-metal-clean test tube and verified gravimetrically using a calibrated mass balance. These samples were then spiked with internal standards comprising known quantities (10 and , respectively) of indium and bismuth (SCP Science) used to correct for instrumental drift. Cd concentrations were measured using a DRC II (dynamic reaction cell; PerkinElmer) axial field ICP-MS at the University of Massachusetts, Boston. The method detection limit (MDL) was calculated according to the two-step approach using the method at 99% confidence interval (). To facilitate comparisons with prior studies, trace metal concentrations were converted from ng/g to based on blood density of . The MDL yielded a value of . The limit of detection (LOD) was calculated as the sensitivity of the methods estimated against an established calibration curve and was . Two samples were below the LOD. The limit of quantification estimated according to Long and Winefordner (1983) was .

In the current study, we analyzed specimens from 20 pregnant women and their offspring, 10 with the highest blood Cd concentrations in the first trimester of pregnancy and 10 who were selected from the majority of women with blood concentrations close to zero (Table 1). Among the 310 women in whom Cd concentrations were measured, those who fell in the top quartile were defined as exposed. To ensure that DNA methylation responsive to the entire spectrum of exposure was captured in the pooled samples, we divided exposed women into quartiles and selected at least two women per quartile, for a total of 10, and classified these samples in this study as the high Cd–exposed group. These were matched as best as was possible to 10 low Cd–exposed women (Table 1), whose Cd concentrations fell below the median for the 310 women measured. We attempted to minimize confounders including maternal body mass index (BMI), age, smoking status, and other factors, which have been shown to either be directly associated with DNA methylation and thus provide a potential explanation for associations observed, or to modify associations between Cd or other toxic metals and DNA methylation. We and others have previously shown associations between DNA methylation and cigarette smoking (Joubert et al. 2016, 2012; Murphy et al. 2012a), maternal obesity (Sharp et al. 2017), diabetes (Quilter et al. 2014), gestational age at delivery (Liu et al. 2013), and sex (Murphy et al. 2012a).

Table 1

Characteristics of study participants.

Characteristic	Low Cd–exposed newborns and mothers; n (%) or median (min–max)	High Cd–exposed mothers; n (%) or median (min–max)	High Cd–exposed newborns; n (%) or median (min–max)
Total n	10	10	9
Maternal age at delivery (years)
≤25	4 (40)	3 (30)	3 (33)
≥26	6 (60)	7 (70)	6 (67)
Prepregnancy BMI
≥30	3 (30)	2 (20)	2 (22)
25≤29	2 (20)	3 (30)	2 (22)
<25	5 (50)	5 (50)	5 (56)
Maternal race/ethnicity
Black non-Hispanic	4 (40)	3 (30)	3 (33)
White or Hispanic	6 (60)	5 (50)	4 (44)
Other	0 (0)	2 (20)	2 (22)
Smoking status
Non-smoker	10 (100)	9 (90)	8 (89)
Smoker	0	0 (0)	0 (0)
Missing	0	1 (10)	1 (11)
Maternal diabetes
Gestational diabetes	0 (0)	1 (10)	1 (11)
Normal	10 (100)	9 (90)	8 (89)
Gestational age at delivery (weeks)
≥37	9 (90)	9 (90)	8 (89)
<37	1 (10)	1 (10)	1 (11)
Infant sex
Male	8 (80)	7 (70)	7 (78)
Female	2 (20)	3 (30)	2 (22)
Birth weight (g)
<2,500	1 (10)	1 (10)	1 (11)
≥2,500	9 (90)	9 (90)	8 (89)
Maternal blood Cd μg/dL, median (minimum–maximum)	0.012 (0.004–0.023)	0.165 (0.090–0.338)	0.161 (0.090–0.338)

Note: Values are the same for the 10 mothers and newborns in the low Cd–exposed group. High Cd–exposed mothers and newborns are shown separately since one newborn was excluded from the analysis due to insufficient DNA. BMI, body mass index.

Offspring of the 20 women were included except one whose cord blood DNA was degraded, leaving 19 samples available for analysis.

DNA Isolation, Bisulfite Conversion, and Whole-Genome Bisulfite Sequencing

Genomic DNA was isolated from buffy coat cells (mixed leukocytes) separated from blood collected in EDTA vacutainers and stored at until extraction. Maternal and newborn cord blood samples were processed using PureGene reagents (Qiagen, Inc.). DNA was bisulfite converted using EpiTect Bisulfite kits (Qiagen, Inc.) according to manufacturer’s protocols. The bisulfite conversion rate was across all samples. Bisulfite-converted maternal DNAs were pooled and converted to sequencing libraries using TruSeq DNA Methylation kits (Illumina, Inc.) according to manufacturer’s protocols. Due to resource limitations at the time, we pooled maternal DNA samples for sequencing. Newborn cord blood samples were sequenced at a later time when we had sufficient resources to sequence them individually. Sequencing libraries from bisulfite converted newborn cord blood DNAs were prepared using TruSeq DNA Methylation kits and index tagged using TruSeq DNA Methylation Index PCR Primers (Illumina, Inc.). The index tagged libraries were quantitated using the ddPCR™ Library Quantification Kit for Illumina TruSeq (Bio-Rad, Inc.) according to manufacturer’s protocols, and a QX200 Droplet Digital PCR system (Bio Rad, Inc.) at the North Carolina State University Genomic Sciences Laboratory (NCSU GSL). Normalized libraries were pooled, with four or five libraries per pool. Libraries were sequenced on an Illumina HiSeq 2500 instrument (Illumina, Inc.) at the NCSU GSL, producing 100-base pair (bp) single-end reads.

Whole-Genome Bisulfite Sequencing Analysis

Bisulfite sequencing reads were aligned to the most recent reference human genome (hg38) using BSMAP (version 2.9; GNU GPL) (Xi and Li 2009), which automatically trims adapters and low-quality sequence reads. Model Based Analysis of Bisulfite Sequencing (MOABS, version 1.2.2; GNU GPL) (Sun et al. 2014) was used to identify differentially methylated CpG dinucleotides and to call differentially methylated regions (DMRs). MOABS uses a beta binomial hierarchical model to capture sampling and biological variation, and it generates a credible methylation difference (CDIF) metric. This metric does not require the high coverage necessary for other approaches, which mainly implement Fisher’s exact test and typically require coverage. MOABS supports three models to call differentially methylated CpGs and DMRs: Fisher’s exact test, hidden Markov model, and CDIF. We used CDIF in this study, which is a single metric combining biological and statistical significance. DMRs were defined by a minimum of three methylated CpG dinucleotides within a maximum interval of 300 bp, each with an absolute methylation difference of at least 10% in the same direction (i.e., lower or higher methylation) between low and high Cd–exposed groups.

Overlap with Genomic Features

The University of California, Santa Cruz (UCSC) Table Browser function was used to download the coordinates of genomic features from the hg38 genome build (Karolchik et al. 2004; Kent et al. 2002). The following datasets were used: Gencode_v22 for gene, exon, and intron coordinates (Harrow et al. 2012); FANTOM5 transcriptional start sites (TSSs) (Lizio et al. 2015); and DNase I hypersensitivity peak clusters from ENCODE (95 cell types) filtered for clusters scoring 200 or above (Thurman et al. 2012). CpG island shores were defined as the regions flanking CpG islands (Irizarry et al. 2009). Intergenic regions were defined as those outside of the other categories (, , exon, intron, CpG island, CpG shore, DNase I hypersensitive site) in Figure 1A. High to intermediate CpG content promoters (HICPs) were obtained from Weber et al. (2007). High CpG content promoters contain a 500-bp region with a CpG ratio above 0.75 and GC content . Low CpG content promoters, which were not used in our study, are defined as not containing a 500-bp region with a CpG ratio above 0.48. Intermediate CpG content promoters are between these categories, and include many CpG islands smaller than 500 bp with a moderate CpG richness and/or a GC content below 55%. Coordinates of HICPs were converted from hg17 to hg19 and then to hg38 using the UCSC Genome Browser LiftOver function. The list of HICPs and their coordinates in hg38 are presented in Table S1. Regions of monoallelic methylation identified in human blood were obtained from Schalkwyk et al. (2010); specifically, 21 single-nucleotide polymorphisms (SNPs) showing strong evidence of an association with allele-specific methylation as demonstrated by an average change in relative allele score of were selected [see Table 1 of Schalkwyk et al. (2010) for a complete list]. Regions of monoallelic methylation were defined as the regions flanking these SNPs. SNRPN, an imprinted gene, was excluded from this analysis. Overlaps between these regions and Cd-associated DMRs were identified using the intersection tool in Galaxy (version 17.01; Galaxy Project). The same tool was used to identify genes overlapping DMRs. CTCF, RAD21, and RNA polymerase II (RNA pol II) (Pol2-4H8 antibody, Abcam) chromatin immunoprecipitation (ChIP)-seq data on H1-hESCs were obtained from ENCODE/HAIB (hg19) and the coordinates of the peaks converted to hg38 using the UCSC Genome Browser LiftOver tool.

Figure 1.

Cd-associated differentially methylated regions (DMRs). (A) Distribution of the 641 DMRs identified in newborn cord blood, the 1,945 DMRs identified in maternal blood, and the 24 DMRs identified in newborn cord blood that overlap DMRs in maternal blood (overlapping), relative to genomic features. Data are expressed as the percentage of DMRs that overlap each category. Categories are not mutually exclusive. DMRs were defined by a minimum of three methylated CpG dinucleotides within a maximum interval of 300 bp, each with an absolute difference of at least 10% in the same direction (i.e., lower or higher methylation) between low and high Cd–exposed groups. The complete list of genes overlapping differentially methylated regions is presented in Table S4. (B) Numbers of hypomethylated and hypermethylated DMRs in newborn cord blood and maternal blood. (C) Intersection of genes overlapping DMRs in newborn cord blood and maternal blood. Ninety-eight genes overlap DMRs in both datasets. Note: TSS, transcriptional start site.

Identification of Human Germline Imprinting Control Regions

To identify likely germline ICRs in humans, i.e., those with differential methylation established in gametogenesis rather than after fertilization in somatic cell types, we searched for evidence in the literature of imprinting and allele-specific methylation at sequences orthologous to the 22 murine germline ICRs. Three of the murine ICRs showed no or inconclusive evidence of imprinting in humans, and we excluded them from the list of likely human ICRs. The Gpr1/Zdbf2 ICR in mice was originally described as a paternally methylated ICR (Kobayashi et al. 2009), but this pattern of methylation has since been shown to be dependent on methylation first established in the oocyte at an alternative promoter of Zdbf2 (Duffié et al. 2014; Kobayashi et al. 2012). Given that an equivalent alternative promoter of human ZDBF2 has not yet been described and the nature of the imprint at this locus remains unclear, we excluded this locus from our analysis. The Rasgrf1 ICR is paternally methylated in mice, but its status as an imprinted locus in humans is debated (Pitamber et al. 2012; Yuen et al. 2011). Thus, this locus was excluded from our analysis. We therefore identified 15 likely human germline ICRs that are maternally methylated and two that are paternally methylated (Table 2). Given that the boundaries of ICRs in humans have not been well defined, we used the coordinates of the overlapping CpG islands as proxies for the positions of the 15 maternally methylated ICRs. The two paternally methylated ICRs (H19/IGF2 and DLK1/MEG3) have relatively low CpG content (observed/expected CpG ratio of 0.4 compared to 0.56 for maternally methylated ICRs) and do not meet the criteria of CpG islands (Schulz et al. 2010); thus, for paternally methylated ICRs, we used coordinates determined experimentally by others (Geuns et al. 2007; Kerjean et al. 2000). Table 2 also includes three loci that may be germline ICRs in humans but do not have imprinted mouse orthologs. These were analyzed separately in our study.

Table 2

Imprinting control regions in mice and humans.

Imprinted murine locus	Imprinted human locus	Reference for imprinting in human	Human ICR coordinates (hg38)	Notes
Maternally methylated ICRs
Peg13	PEG13	(Court et al. 2014)	chr8:140097739-140100885	—
Kcnq1ot1	KvDMR	(Smilinich et al. 1999)	chr11:2699181-2700857	—
Igf2r	IGF2R	(Xu et al. 1993)	chr6:160005233-160006470	—
Gnasxl	GNASXL	(Hayward et al. 1998)	chr20:58853970-58856184	—
Nnat	NNAT	(Evans et al. 2001)	chr20:37520202-37521734	—
Gnas Ex1A	GNAS Ex1A	(Hayward et al. 1998)	chr20:58888598-58892684	—
Mcts2	MCTS2	(Wood et al. 2007)	chr20:31547274-31547489	—
Zac1/Plagl1	PLAGL1	(Kamiya et al. 2000)	chr6:144007780-144008710	—
U2af1-rs1	—	—	—	No imprinting of human ortholog
Inpp5f_v2	INPP5F_V2	(Wood et al. 2007)	chr10:119818018-119818873	—
Impact	—	—	—	No imprinting of human ortholog
Grb10	GRB10	(Blagitko et al. 2000)	chr7:50782056-50783174	—
Nap1l5	NAP1L5	(Wood et al. 2007)	chr4:88697614-88697903	—
Cdh15	—	—	—	No imprinting of human ortholog
Mest	MEST	(Kobayashi et al. 1997)	chr7:130490899-130493270	—
Zim2/Peg3	ZIM2/PEG3	(Murphy et al. 2001)	chr19:56839916-56840916	—
Snurf/Snrpn	SNURF/SNRPN	(Sutcliffe et al. 1994)	chr15:24954889-24955907	—
Peg10	PEG10	(Ono et al. 2001)	chr7:94655547-94657215	—
Gpr1/Zdbf2	GPR1/ZDBF2	(Kobayashi et al. 2009)	—	Nature of ICR in humans unclear; exclude from analysis
Paternally methylated ICRs
H19/Igf2	H19/IGF2	(Giannoukakis et al. 1993)	chr11:1999841-2000164	—
Dlk1/Gtl2	DLK1/MEG3	(Kobayashi et al. 2000)	chr14:100810900-100811156	—
Rasgrf1	—	—	—	Nature of ICR and imprinting status in humans unclear; exclude from analysis
Human-specific ICRs
—	RB1	(Kanber et al. 2009)	chr13:48318500-48319721	—
—	DIRAS3	(Yu et al. 1999)	chr1:68050554-68050790	—
—	L3MBTL1	(Li et al. 2004)	chr20:43514571-43514951	—

Note: To identify likely germline imprinting control regions (ICRs) in humans, we searched for evidence in the literature of imprinting and allele-specific methylation at sequences orthologous to the 22 murine ICRs. —, no information; GPR1/ZDBF2 and RASGRF1 were excluded from the analysis, as described in the Methods.

Overlap with CpG Dinucleotides on the Infinium MethylationEPIC BeadChip (Illumina, Inc.)

Coordinates for the CpG dinucleotides represented on this array (hg19) were converted to hg38 using the UCSC Genome Browser LiftOver tool and intersected with the Cd-associated DMRs identified by WGBS.

Functional Enrichment Analysis

Enrichr was used to intersect genes that overlapped DMRs with the Database of Genotypes and Phenotypes (dbGaP) and with the Mouse Genome Informatics (MGI) Mammalian Phenotype 2017 database to identify enriched functional categories (Chen et al. 2013; Kuleshov et al. 2016). The Database for Annotation, Visualization and Integrated Discovery (DAVID; version 6.8) was also used to identify enriched disease categories using the same set of input genes (Huang et al. 2009a, 2009b). Specifically, GAD (Genetic Association Database) Diseases and GAD Disease Classes were examined with DAVID. Ingenuity Pathway Analysis (IPA) was also used to identify common upstream regulators using the same input genes (Krämer et al. 2014).

Results

Study Participants

Table 1 summarizes the characteristics of the 20 pregnant women and their offspring. The high and low Cd–exposed groups have similar distributions with respect to maternal age, race/ethnicity, infant sex, birth weight, gestational age at delivery, prepregnancy obesity, and smoking during gestation. In this population, only one woman had gestational diabetes, a factor that may also influence DNA methylation. The median blood Cd levels in first-trimester maternal blood were (range: ) in the high Cd–exposed group compared with (range: ) in the low Cd–exposed group. One female infant born to a non-Hispanic white woman who was aged was excluded due to the limited quantity of available DNA. For this reason, the characteristics of the infants in the high Cd–exposed group are also presented separately in Table 1.

Whole-Genome Bisulfite Sequencing

For newborn cord blood samples, 460 and 535 million reads were generated for low- and high-Cd exposure groups, respectively, of which 61.5% and 61.1% were uniquely aligned to the reference genome (Figure S1A), representing 13 times and 15 times coverage. For maternal blood samples, 201 million and 188 million reads were generated for low- and high-Cd exposure groups, respectively, of which 56.2% and 59.5% of reads were uniquely aligned to the reference genome (Figure S1B), representing 5 times and 6 times coverage. The percentage of uniquely aligned reads is consistent with data reported in a separate study in which four alignment tools for genome-wide bisulfite sequencing data were compared and yielded 55.5% uniquely mapped reads, on average (Chatterjee et al. 2012).

Distribution of Differentially Methylated Regions

A methylation difference of at least 10% between high Cd–exposed and low Cd–exposed groups was evident for 169,691 and 276,452 CpG dinucleotides in newborn cord blood and maternal blood, respectively (Figures S2A and S2B). In newborn cord blood, 641 Cd-associated DMRs were identified, defined as regions containing at least 3 CpG dinucleotides within a 300-bp window, each showing a methylation change in the same direction (i.e., hypomethylation or hypermethylation). Table S2 shows the complete list of DMRs identified in newborn cord blood.

The DMR size ranged from 6–431 bp (Figure S2C). Figure 1A shows the distribution of the 641 DMRs relative to genomic features. Of note, 66% of the DMRs overlapped with gene bodies or with the flanking regions of genes, 24% were within of an annotated TSS, and a similar proportion overlapped CpG islands. Forty-one percent of DMRs overlapped DNase I hypersensitive sites. Only 15% of DMRs were intergenic. Together, these data suggest that Cd-associated DMRs are localized at genes and regulatory elements.

In maternal first-trimester blood samples, 1,945 Cd-associated DMRs were identified (see Table S3 for a complete list). These ranged from 6–555 bp in length (Figure S2D). Of the 1,945 DMRs, 68% overlapped with gene bodies or the flanking regions of genes, 29% were within of TSSs, 27% overlapped CpG islands, and 36% were at DNase I hypersensitive sites (Figure 1A).

Applying the same criteria for calling DMRs, DNA from first-trimester maternal blood contained more than three times the number of DMRs identified in newborn cord blood (1,945 DMRs in maternal blood, 641 in newborn cord blood). However, the proportions of DMRs situated near to genes, TSSs, CpG islands, and other genomic features were similar for maternal and newborn cord blood samples (Figure 1A).

Four percent of DMRs in newborn cord blood (24/641) overlapped with DMRs in maternal blood (Table S1), identifying loci that may be susceptible to Cd irrespective of the timing of exposure. These sites were enriched at genes and regulatory elements compared to all newborn and maternal DMRs, particularly at TSSs, CpG islands, CpG island shores, and DNase I hypersensitive sites (Figure 1A), although this should be interpreted with caution given the small number of overlapping DMRs.

Cd exposure was associated with both hypomethylated and hypermethylated loci. However, while there was a nearly equal proportion of hypomethylated and hypermethylated DMRs in newborn cord blood (51% hypomethylated, 49% hypermethylated), maternal blood DMRs were mostly hypermethylated (64%; Figure 1B).

Methylation at Imprinting Control Regions

We next used our genome-wide data to investigate the relationship between Cd exposure and methylation at ICRs.

First, we identified 17 likely human germline ICRs, based on evidence from the literature and orthologous sequences in the mouse that function in this capacity (Table 2), as described in the Methods section. This includes 15 maternally methylated and two paternally methylated ICRs. Of the 15 maternally methylated germline ICRs, four (IGF2R, KvDMR, SNURF/SNRPN, GNASXL) overlapped Cd-associated DMRs in newborn cord blood, nine (SNURF/SNRPN, GNASXL, PEG3, MCTS2, NNAT, GNAS, PLAGL1, GRB10, PEG13) overlapped DMRs in maternal blood, and four did not overlap with DMRs in either newborn cord blood or maternal blood (Figure 2A). Some ICRs were hypomethylated and others hypermethylated, while the GNASXL ICR in newborn cord blood and the PLAGL1 ICR in maternal blood each contained regions of both hypo- and hypermethylation. Neither of the two paternally methylated ICRs overlapped DMRs.

Figure 2.

Cd-associated differentially methylated regions (DMRs) at imprinting control regions (ICRs). (A) Hypomethylation and hypermethylation at ICRs in newborn cord blood and maternal blood. The GNASXL ICR in newborn cord blood and the PLAGL1 ICR in maternal blood overlap with both hypomethylated and hypermethylated DMRs, indicated by the dual shading. Human-specific ICRs: These have no known imprinted orthologs in mice, and their status as germline or somatic ICRs has not been resolved. Refer to the Methods for how these ICRs were identified, and to Table 2 for details of the ICRs. (B) Cd-associated hypermethylation at the KvDMR ICR. The height of the vertical black lines in the boxed portion of the figure represents percent methylation at individual CpG dinucleotides (0–100%) across all samples in each group. The hypomethylated DMR identified by whole-genome bisulfite sequencing (WGBS) is near the 3‘end of the ICR (not highlighted). The red box highlights the region of differential methylation between high and low Cd–exposed groups in both newborn cord blood and maternal blood at the RNA pol II binding site of the KCNQ1OT1 transcriptional start site (TSS), identified by visual inspection. (C) Cd-associated hypermethylation at the GRB10 ICR in maternal blood. The region of differential methylation identified by WGBS, indicated by a red box, overlaps CTCF and RAD21 binding sites, and is coincident with the TSS of the paternally expressed transcript of GRB10 (GRB10-pat). Note: GRB10-mat, maternally expressed transcript of GRB10.

To determine whether maternally methylated ICRs are significantly enriched for Cd-associated DMRs, we first compared these ICRs to the total number of CpG islands in the genome that overlap at least one DMR. For both newborn cord blood and maternal blood, the enrichment was highly significant ( for newborn cord blood; for maternal blood; Fisher’s exact test, two-sided). However, maternally methylated ICRs, which are all located at promoters, have a very high CpG density; because the CpG density of CpG islands is variable, we therefore compared ICRs to a more appropriate set of control sequences than all CpG islands. We used a subset of 11,168 human promoter sequences that are classified as HICPs, which have comparable CpG content to maternally methylated ICRs (Schulz et al. 2010; Weber et al. 2007). In newborn cord blood, only 0.3% (30/11,168) of HICP sequences overlapped DMRs, whereas 27% (4/15) of maternally methylated ICRs overlapped (, Fisher’s exact test, two-sided). In maternal blood, 1.4% (153/11,168) of HICP sequences overlapped DMRs, whereas 60% (9/15) of maternally methylated ICRs overlapped ().

Our list of 17 human germline ICRs is likely to be a conservative estimate of the total number, since the human genome contains imprinted loci that are absent in mouse. This includes the RB1, DIRAS3, and L3MBTL1 loci, for which regions of allele-specific methylation have been identified, but their status as germline or somatic ICRs has not been resolved (Kanber et al. 2009; Li et al. 2004; Niemczyk et al. 2013). In maternal blood, all three of these ICRs overlapped directly with DMRs, while there was no overlap in newborn cord blood (Figure 2A).

In addition to ICRs, other regions of the genome display monoallelic DNA methylation, but the methylated allele is not determined by parental origin. To test whether the observed enrichment of Cd-associated DMRs is a property of monoallelically methylated regions in general or of ICRs specifically, we looked at the overlap between Cd-associated DMRs and monoallelically methylated, but not imprinted, loci (Schalkwyk et al. 2010). We identified no overlap between these loci and Cd-associated DMRs in newborn cord blood or maternal blood.

Consistent with an enrichment of DMRs at imprinted loci, the most significant associations of genes overlapping DMRs with phenotypes in the MGI Mammalian Phenotype 2017 database are maternal imprinting and abnormal imprinting for newborn cord blood and maternal blood, respectively (Table S2). Furthermore, IPA identified DNMT3L as one of the top five upstream regulators of genes overlapping DMRs in newborn cord blood ( as calculated by IPA using the right-tailed Fisher’s exact test) (Table 3). DNMT3L is a cofactor of the de novo methyltransferases DNMT3A and DNMT3B, and is required for the establishment of methylation in oocytes at maternally methylated ICRs (Bourc’his et al. 2001). Additionally, in maternal blood, the methylation-sensitive insulator protein CTCF, which maintains allele-specific expression at multiple imprinted loci (Prickett et al. 2013), and its binding partner RAD21, were identified as being two of the top five upstream regulators by IPA ( for CTCF; for RAD21). Other upstream regulators identified with this approach are not directly implicated in the regulation of imprinting, but several have been associated with gene regulation more generally (Table 3).

Table 3

Enriched upstream regulators of differentially methylated genes, from Ingenuity Pathway Analysis (IPA).

Top upstream regulators	p-Value of overlap	Function	Reference
Newborn cord blood
miR-210	5.01×10^−4	miRNA, posttranscriptional regulator	(Bavelloni et al. 2017)
NEAT1	0.00129	Long noncoding RNA, possible transcriptional regulator	(Yu et al. 2017)
DNMT3L	0.00213	Required for establishment of maternal genomic imprints	(Bourc’his et al. 2001)
CREBBP	0.00224	Transcriptional coactivator and acetyltransferase	(McManus and Hendzel 2001)
FSH	0.00262	Gonadotropin	(Rose et al. 2000)
Maternal blood
Actinonin	1.09×10^−6	Antibacterial agent	(Chen et al. 2000)
LRPPRC	3.93×10^−6	Transcriptional regulator	(Siira et al. 2017)
RAD21	9.45×10^−6	Regulates chromatin conformation; colocalizes with CTCF	(Nativio et al. 2009)
DAP3	1.20×10^−5	Regulator of mitochondrial translation	(Xiao et al. 2015b)
CTCF	2.89×10^−5	Regulates chromatin conformation and gene expression; controls allele-specific expression of imprinted genes	(Franco et al. 2014)

Note: Genes overlapping DMRs were input into IPA and top upstream regulators identified. p-Values are calculated using the right-tailed Fisher’s exact test.

Predicted Impacts of Cd-Associated DNA Methylation Changes on Imprinting Regulation

We next used our WGBS data to make predictions about how Cd-associated changes to DNA methylation at ICRs might impact the regulation of imprinting. We aligned our data against experimentally determined binding sites for CTCF and RAD21—two regulators of imprinted gene expression—and RNA pol II. For analysis, we focused on two ICRs (KvDMR and GRB10) that regulate genes involved in intrauterine growth and glucose homeostasis, two subclinical indicators of altered metabolic function known to be affected by Cd exposure (Bell et al. 1990; Xu et al. 2016).

At the KvDMR ICR, we identified a hypomethylated DMR in newborn cord blood (Figure 2A), but this did not align with binding sites for CTCF, RAD21, or RNA pol II (data not shown). However, visual examination of the methylation state across the entire ICR revealed a region of Cd-associated hypermethylation in both newborn cord blood and maternal blood that overlapped the binding site of RNA pol II at the KCNQ1OT1 promoter (Figure 2B). Transcription of this long, noncoding RNA is important for the regulation of a cluster of imprinted genes with roles in intrauterine growth and glucose metabolism, as determined by in vivo manipulation of its expression (Mohammad et al. 2010; Redrup et al. 2009). We would predict that hypermethylation at the promoter would reduce KCNQ1OT1 transcription by inhibiting RNA pol II binding, potentially impacting imprinted expression across the locus.

At the GRB10 ICR, the Cd-associated DMR identified in maternal blood is at the promoter of an alternative transcript of GRB10 (Figure 2C). This region overlaps with the binding sites of CTCF and RAD21, which have been predicted to regulate imprinting at the locus (Plasschaert and Bartolomei 2015; Prickett et al. 2013). We would predict that Cd-associated hypermethylation would abrogate CTCF binding, which is methylation sensitive, leading to a relaxation of imprinting and an altered dosage of GRB10. Further work, potentially using animal models, will be required to empirically test these hypotheses.

Differential Methylation at the RNA5S Cluster

Changes to DNA methylation at ribosomal DNA (rDNA) have been associated with environmental insult during development, specifically protein restriction in mice (Holland et al. 2016). In our study, multiple Cd-associated DMRs were identified at the RNA5S cluster on chromosome 1 (Figure S3). In each of newborn cord blood and maternal blood, seven DMRs were identified within the cluster, all of which were hypomethylated in the high Cd–exposed group (Tables S2 and S3).

Functional Analysis of Differentially Methylated Genes Genome-Wide

To understand how methylation changes genome-wide might contribute to Cd-associated disease, we identified all annotated genes that directly overlapped DMRs, including imprinted and biallelically expressed genes. We identified 444 genes in newborn cord blood and 1,340 genes in maternal blood, of which 98 were common between the two (Figure 1C; Table S4).

These differentially methylated genes were analyzed for common functional annotations using Enrichr. Specifically, genes were compared to the dbGaP. The five categories most significantly associated with differentially methylated genes in maternal blood, i.e., those with the smallest adjusted p-value after applying the Benjamini-Hochberg method for multiple testing correction were BMI (), blood pressure (), body weight (), albumins (), and hematocrit (0.0201) (Table 4). In newborn cord blood, BMI, atrial fibrillation, hypertension, heart function tests, and sleep were the five most significantly associated functional categories, although these were not significant after multiple testing correction (Table 4).

Table 4

The five dbGaP (Database of Genotypes and Phenotypes) categories most significantly associated with differentially methylated genes in newborn cord blood, maternal blood, and common to both, ranked by adjusted p-value (adj p-value).

dbGaP category	p-Value	Adj p-value
	Newborn cord blood
Body mass index	8.3×10^−4	0.0988
Atrial fibrillation	0.0026	0.0988
Hypertension	0.0029	0.0988
Heart function testsa	0.0041	0.1094
Sleep	0.0073	0.1611
	Maternal blood
Body mass index	9.0×10^−8	2.0×10^−5
Blood pressure	3.5×10^−7	3.8×10^−5
Body weight	1.9×10^−5	0.0014
Albumins	9.2×10^−5	0.0051
Hematocrit	0.0006	0.0201
	Common to newborn cord blood and maternal blood
Heart failure	0.0135	0.1758
Platelet count	0.0288	0.1758
Forced expiratory volume	0.0319	0.1758
Monocyte chemoattractant protein 1	0.0329	0.1758
Respiratory function testsa	0.0382	0.1758

Note: Differentially methylated genes were analyzed using Enrichr, as described in the Methods. dbGaP, Database of Genotypes and Phenotypes.

The heart function tests and respiratory functions tests categories include genes that have been associated with affecting the outcomes of measurements for normal heart and respiratory function, respectively.

Consistent with the results from Enrichr, analysis of differentially methylated genes in maternal blood and newborn cord blood using DAVID identified associations with various metabolic and cardiovascular functions, among others (see Tables S3 and S4 for full results).

For the 98 differentially methylated genes common between newborn cord blood and maternal blood, the top functional annotation using dbGaP was heart failure. Other functions most significantly associated with these genes are presented in Table 4. However, none of these associations were statistically significant after applying a multiple testing correction, which may reflect the small number of genes analyzed.

Discussion

We undertook a genome-wide analysis of DNA methylation to identify CpG dinucleotides that are associated with exposure to Cd in adult women and their newborn offspring. Most regions of differential methylation were at functional elements of the genome, suggesting they may affect genome regulation. This included an enrichment at genes involved in cardiovascular and metabolic functions, suggesting that epigenetic changes may form part of the mechanism through which Cd causes disease and that disease susceptibility may be programmed in early life following in utero exposure. Our findings also suggest that ICRs may be particularly sensitive to Cd exposure.

ICRs have been proposed and are widely used as biomarkers for epigenetic modulation by the environment (Hoyo et al. 2009). Whether ICRs are more or less sensitive to perturbation than other loci has not been resolved. Since the differential methylation of ICRs established in the gametes must survive mitosis and differentiation in the embryo, retaining a memory of allelic origin, it could be argued that their methylation state is robust and likely to be protected from environmental insult by multiple layers of regulation. On the other hand, the unique methylation dynamics of ICRs in early development, which retain their methylation state while most of the genome undergoes methylation erasure (Cowley and Oakey 2012), may make them uniquely susceptible to perturbation during this time window. Our data suggest that maternally methylated ICRs may be more susceptible than similar sequences in the genome to DNA methylation changes associated with in utero Cd exposure. Although we focused predominantly on 15 maternally methylated loci that are likely to be true germline ICRs, three additional ICRs (RB1, DIRAS3, L3MBTL1), which are also methylated on the maternally inherited allele but whose status as germline or somatic ICRs has not been fully resolved, also overlapped with DMRs in maternal blood, adding further support to the idea that DNA methylation at ICRs is sensitive to Cd. Our findings also suggest that the possible increased sensitivity is not a property of monoallelically methylated regions in general, but is specific to ICRs.

In a mouse model of in utero nutrient restriction, imprinted genes were neither more nor less affected than other genes at the level of transcription (Radford et al. 2012). A subset of imprinted genes including Peg3 showed altered expression, but this was not accompanied by DNA methylation changes at the Peg3 ICR. There are several potential explanations for the different findings of Radford et al. and our own study, and the observations are not necessarily inconsistent. First, ICRs in mice and humans may respond differently to exposure, although mechanisms of imprinting and DNA methylation dynamics in development are highly conserved (Bartolomei and Ferguson-Smith 2011; Smith et al. 2014). Second, the differences could reflect the nature of the exposures. However, like maternal nutrient restriction, Cd likely causes IUGR by impacting placental function (Xu et al. 2016) rather than by acting directly on the fetus, since Cd is inefficiently transported across the placenta (Lauwerys et al. 1978). Alternatively, the differences could reflect the timings of exposure. We selected study participants based on maternal blood Cd measured in the first trimester, consistent with fetal exposure during very early development. Radford et al. (2012) induced nutrient restriction in the mouse during the third week of gestation, after the postfertilization genome-wide methylation erasure from which ICRs are protected. This time window may represent a critical period during which ICRs are susceptible to epigenetic modulation by the environment. Further work will be required to determine how environmental stressors may interact with ICRs and their protection machinery at this stage.

Many imprinted genes regulate growth and metabolism. Enrichment of DMRs at ICRs suggests that perturbation of imprinted gene expression might contribute to growth restriction and metabolic disease associated with Cd exposure. To predict how DMRs might affect imprinting regulation, we aligned our data against CTCF, RAD21, and RNA pol II ChIP-seq data. We focused on two ICRs: the KvDMR and GRB10 ICRs. Methylation at the KvDMR controls the imprinted transcription of a long noncoding RNA, KCNQ1OT1, which is required for the imprinted expression of several growth regulators, including PHLDA2 and CDKN1C, whose expression levels are negatively correlated with birth weight (Apostolidou et al. 2007; Piyasena et al. 2015). The hypermethylation at the KCNQ1OT1 promoter and RNA pol II binding site observed in this study would be expected to reduce transcription of the noncoding RNA, potentially relaxing the imprinting state across the locus. This hypothesis is supported by a mouse model of Kcnq1ot1 deletion, which shows loss of Cdkn1c and Phlda2 imprinting (Mohammad et al. 2010). Altered dosage of PHLDA2 and CDKN1C could contribute to the growth restriction associated with Cd exposure. The GRB10 ICR was differentially methylated only in maternal blood. Mouse studies have shown that during development, Grb10 is a potent growth suppressor (Charalambous et al. 2003; Cowley et al. 2014), while in adulthood, it negatively regulates insulin signaling, controlling glucose homeostasis (Smith et al. 2007). We predict that the hypermethylation observed in maternal blood at the CTCF/RAD21 binding site at GRB10 might cause loss of imprinted expression, potentially contributing to the Cd-associated effects on glucose metabolism observed in experimental models. These observations provide a proof of principle that our WGBS data can be used to make predictions on how Cd-associated methylation changes may affect the regulation of imprinting, potentially contributing to Cd-associated growth and metabolic deficiencies. Further studies will be required to empirically test these hypotheses and determine the contribution of methylation changes at ICRs to Cd-associated disease.

We identified hypomethylation across a cluster of genes on chromosome 1 encoding the 5S rRNA subunit. Studies in mice, Drosophila, and yeast have identified rDNA as a target of environmental stress, causing altered rDNA copy number (Aldrich and Maggert 2015; Jack et al. 2015; Shiao et al. 2011). Recently, mouse models of nutrient restriction were found to display differential methylation at 45S rDNA in a manner dependent on underlying genetic variation (Holland et al. 2016). To our knowledge, our study is the first to suggest that rDNA may be a target for epigenetic modulation by an environmental stressor in humans. Our preliminary findings, together with experimental data from mice, suggest that epigenetic changes at rDNA may represent a general response to multiple environmental insults, with the potential to impact cellular metabolism (Holland 2017).

Analyses of all differentially methylated genes using Enrichr and DAVID suggested an enrichment for functions associated with metabolic and cardiovascular processes. Although these associations were not significant for newborn cord blood, they suggest that developmental exposure to Cd might program DNA methylation changes during development that persist into adulthood, potentially contributing to disease susceptibility in later life. While newborn cord blood DMRs were hypomethylated and hypermethylated in approximately equal proportions, 64% of DMRs in maternal blood were hypermethylated. This may reflect differences between direct and indirect Cd action on the epigenome; since Cd is inefficiently transported across the placenta (Lauwerys et al. 1978), effects observed in newborn cord blood may be a consequence of altered placental function, rather than the direct effects observed following adult exposure.

Previously, DNA methylation changes associated with in utero Cd exposure have been measured using array-based approaches, capturing a small proportion of the genome. We performed WGBS to assay methylation across the entire genome, providing the opportunity to identify DMRs outside of regions profiled on arrays such as the Infinium MethylationEPIC BeadChip, which interrogates CpG dinucleotides at CpG islands and other functional regions. Of the DMRs we identified in newborn cord blood and maternal blood by WGBS, only 18% and 21%, respectively, overlapped with at least one CpG dinucleotide on this array. Thus, many of the DMRs we identified by WGBS would have been missed if an array-based approach had been used. In most cases, our DMRs do overlap with or are within 250 bp of a CpG island or other feature represented on the Infinium MethylationEPIC BeadChip; however, only a subset of CpG dinucleotides at those features are measured on the array, frequently missing our relatively small DMRs. A small, localized change in DNA methylation may be sufficient to have a functional effect, such as through interference with transcription factor binding. Thus, while arrays provide a standardized platform for assaying DNA methylation that allows for more reliable comparisons to be made between samples processed in different laboratories than WGBS, for which methods are variable, our global approach highlights the potential for array-based methods to miss a large proportion of methylation changes associated with an exposure.

We identified DMRs in DNA isolated from mixed leukocytes and did not adjust for potential variation in blood cell type composition, since ICR methylation is established in gametogenesis and is thought to be consistently maintained through development in all cell lineages. This consistency of methylation at ICRs has been demonstrated for multiple somatic tissues, including for major newborn cord blood fractions (Murphy et al. 2012b; Pervjakova et al. 2016; Woodfine et al. 2011).

Our findings should be interpreted in the context of the small sample size and limitations inherent in human data. Although the sample size precluded adjustment by known confounding factors, we matched as best as was possible sex, maternal smoking status, birth weight, and other factors to minimize potential confounding. However, the Centers for Disease Control and Prevention (CDC) have demonstrated that the average American is exposed to chemicals in addition to nonchemical stressors such as anxiety and stress, all of which may also alter the epigenome. We cannot exclude the possibility that the reported methylation differences may be related to other exposures. As methods to measure multiple exposures become more accessible, the potential effects of other exposures may be clarified. We also cannot exclude a possible role for genetic variation in modulating epigenetic susceptibility. An improved understanding of gene-by-environment interactions in this context, and the study of larger cohorts, will be important for making further progress in the field. Furthermore, the higher read depth we obtained from sequencing individual newborn cord blood samples than from sequencing the pool of maternal samples means that greater confidence may be assigned to the findings from cord blood. These limitations notwithstanding, our data suggest that exposure to Cd may be associated with epigenetic perturbations at genes that map to disease pathways previously linked to this ubiquitous toxic metal.

Conclusions

Our findings suggest a possible role for DNA methylation in contributing to the cardiovascular and metabolic outcomes that have been associated with Cd exposure. Our data also suggest that ICRs may be more sensitive to Cd exposure than similar, nonimprinted loci in the genome, motivating further study of these regulatory sequences to understand their importance in the response to environmental stress.

Click here for additional data file.

Click here for additional data file.