Copper associates with differential methylation in placentae from two US birth cohorts

ABSTRACT

Copper is an essential trace nutrient and an enzymatic cofactor necessary for diverse physiological and biological processes. Copper metabolism is uniquely controlled in the placenta and changes to copper metabolism have been linked with adverse birth outcomes. We investigated associations between patterns of DNA methylation (DNAm; measured at >485 k CpG sites) and copper concentration measured from placentae in two independent mother-infant cohorts: the New Hampshire Birth Cohort Study (NHBCS, n = 306) and the Rhode Island Child Health Study (RICHS, n = 141). We identified nine copper-associated differentially methylated regions (DMRs; adjusted P < 0.05) and 15 suggestive CpGs (raw P < 1e-5). One of the most robust variably methylated CpGs associated with the expression of the antioxidant, GSTP1. Our most robust DMR negatively associates with the expression of the zinc-finger gene, ZNF197 (FDR = 4.5e-11). Genes co-expressed with ZNF197, a transcription factor, are enriched for genes that associate with birth weight in RICHS (OR = 2.9, P = 2.6e-6, N = 194), genes that are near a ZNF197 consensus binding motif (OR = 1.34, P = 0.01, N = 194), and for those classified in GO biological processes growth hormone secretion (P = 3.4e-4), multicellular organism growth (P = 3.8e-4), and molecular functions related to lipid biosynthesis (P = 1.9e-4). Further, putative transcriptional targets for ZNF197 include genes involved in copper metabolism and placentation. Our results suggest that copper metabolism is tied to DNAm in the placenta and that copper-associated patterns in DNAm may mediate normal placentation and foetal development.

KEYWORDS:

• Copper
• copper metabolism
• DNA methylation
• placenta
• placental copper
• placental epigenetics
• placental methylation

Introduction

Copper is an essential trace nutrient, utilized in many cuproenzymes and vital to a diverse array of cellular functions that are crucial to life. Copper deficiency, while uncommon in the United States, can occur when dietary sources of copper (yeast bread, potatoes, tomatoes, ready-to-eat cereal, and beef) are scarce, or when intestinal absorption is impaired. Copper toxicity is extremely rare, but can result when copper dust is inhaled (e.g., when mining or grinding copper-containing metal), or when acidic water or food is exposed to copper and ingested [1]. High copper concentrations are toxic, necessitating strict copper ion homoeostasis. Menkes and Wilson diseases, which result from mutations in the copper metabolism genes ATP7A and ATP7B, respectively, demonstrate the essentiality of copper and copper metabolism in humans [2].

This balance is particularly critical during in utero development when copper ion homoeostasis is required for proper foetal development and placental function; variation in copper concentration in placenta [3–6], maternal serum [7–11], and cord blood [7,9] has been associated with adverse birth outcomes, including low birth weight. As one of the only human tissues that expresses both ATP7A and ATP7B, the placenta uniquely controls copper metabolism and foetal access to copper. ATP7A controls copper transfer from the placenta to foetal the circulation, while ATP7B controls the return of excess placental copper to maternal circulation [12]. Likely in coordination with changing hormone levels, placental and cord blood copper levels appear to increase in the third trimester to accommodate the storage of copper in the foetal liver [13,14]. Although the mechanism by which copper is transferred from maternal to foetal circulation is not yet fully elucidated, maternal serum copper is higher than umbilical cord serum copper, suggesting that the placenta acts as a partial barrier [9,13]. Unlike the heavy metal cadmium, copper does not appear to accumulate in the placenta under maternal excess [5,15,16]. Loss of function in copper metabolism genes can result in copper accumulation in the placenta and miscarriage [12,17,18].

Though the specific biological processes that associate with variations in copper concentration have not been thoroughly studied in the placenta, studies of other tissues provide evidence that DNA methylation is involved. Animal models of Wilson disease suggest that liver damage associated with copper accumulation may be partially mediated by global hypomethylation [19–21]. Through inhibition of S-Adenosylhomocysteine Hydrolase (SAHH), a key enzyme that regulates the amount of s-adenosylmethionine (SAM) available for methylation reactions, copper accumulation reduces a cell’s methylation potential [20,22]. In mice, global hypomethylation led to a decreased ability for liver cells to respond to oxidative stress [21]. In this model, copper chelation or methyl donor supplementation restored global methylation levels and improved liver function [19,21]. It remains unclear how variability in placental copper concentrations impacts placental functions and to what extent changes in DNA methylation patterns are involved.

Although copper metabolism is tightly controlled by necessity, it can be influenced by endogenous factors beyond the overt alterations observed in Mendelian disease. Recent research in placental cells has found that insulin decreases the action of ATP7B and increases the action of ATP7A [23]. Accordingly, studies of type I and II diabetes have found altered copper metabolism in cases vs. controls and rat models of diabetes have higher tissue concentrations of copper [24]. In Cohen diabetic rats, placental hypomethylation is ameliorated by copper supplementation [6]. Taken together, these findings support a scenario in which copper metabolism in the placenta could become uncoupled by environmental factors, leading to changes in the epigenetic landscape of the placenta. To more broadly consider the impact of variation in placental copper concentration on DNA methylation in a human population, we assess site-specific and regional changes to DNA methylation that associate with placental copper concentration and how these changes associate with the birth weight in two independent cohorts using a meta-analytic framework.

Results

This study analysed data from 306 mother–infant pairs from the New Hampshire Birth Cohort (NHBCS) and 141 mother–infant pairs from the Rhode Island Child Health Study (RICHS), with the demographics of the participants for this study shown in Table 1. A few notable differences exist between the two cohorts. RICHS contains more small for gestational age (SGA) and large for gestational age (LGA) newborns than NHBCS by design (Table 1). RICHS samples have slightly higher median copper concentration (880.1 ng/g vs 847.1 ng/g in NHBCS, Table 1). Lastly, the RICHS cohort is slightly more racially diverse than the NHBCS cohort (Table 1).

Table 1. Participant demographics.

	NHBCS ^α	RICHS ^α
Maternal Characteristicsβ
Educational attainment
●High school or less	59.5% (182)	48.2% (68)
●More than high school	40.5% (124)	51.8% (73)
Age (yrs)	31.7 (19–44.8)	30 (18–40)
BMI(kg/m2)	16.6–48.2 (24.9)	16–45.4 (25.6)
Pregnancy Smoking status
●Any	12.7% (39)	13.5% (19)
●None	87.3% (267)	86.5% (122)
Infant Characteristicsβ
Race γ
●Not white	0.3% (1)	34% (48)
●White	99.7% (305)	63.8% (90)
Gestational Age (wks)	39.6 (31.6–42.1)	39 (37–41)
Birthweight (g)	3452.5 (1680–4912)	3600 (2160–5465)
Birthweight Category
●SGA	4.9% (15)	17.7% (25)
●AGA	85.9% (263)	45.4% (64)
●LGA	8.5% (26)	36.9% (52)
Sex
●Male	52.6% (161)	51.8% (73)
●Female	47.4% (145)	48.2% (68)
Placental copper (ng/g)	847.1 (602.7–2428.5)	880.1 (623.1–2643.5)
Below experiment-wide median copper concentration	53% (162)	43% (61)

^α NHBCS and RICHS have 306 and 141 samples, respectively.

^β Where data are continuous, median and range are given; per cent and count are given for count data.

^γ Three samples were missing from NHBCS for this measure.

To analyse the associations between placental copper concentration and genome-wide epigenetic landscape, we utilized an epigenome-wide association study, or EWAS, design. Within each cohort and for each CpG individually, methylation M-values for all placental samples were regressed on corresponding dichotomized placental copper concentrations, while also adjusting for infant sex, maternal pregnancy smoking status, highest achieved maternal education level, and the estimated proportions of putative placental cell types. The genomic inflation factors (median (T-statistic)²/0.4549) for RICHS (1.2) and NHBCS (0.94) EWAS, respectively, did not indicate substantial statistical confounding. The RICHS EWAS produced 21 CpGs with suggestive P-values (p < 10⁻⁵) and two CpGs that were significant after Bonferroni correction (p < 10⁻⁷, Figure 1a). The NHBCS EWAS produced 9 CpGs with suggestive P-values (Figure 1b). These two EWAS were combined via meta-analysis, which produced 15 CpGs with suggestive P-values for association with placental copper status (Figure 1c).

Figure 1. Epigenome-wide association study results for placental copper status.

Volcano plots display the p-value and effect estimate (log fold change) for the a. RICHS, b. the NHBCS cohort-level associations, as well as for c. the meta-analysis. d. The Manhattan plot of meta-analysis results indicates -log₁₀(p-value) for individual CpG associations with placental copper status, as well as locations on the chromosomes (indicated on x-axis). Red vertical lines indicate differentially methylated regions found by CombP. Blue dots represent CpGs within the DMR.

In both the cohort-level EWASs and the meta-analysis, suggestive and significant hits were spread throughout the genome (Figure 1d). Copper status was more often associated with lower methylation (hypomethylation) than higher methylation (hypermethylation) among top results for all three analyses (Figure 1a-c; Tables S1-S3). Among CpG sites passing a liberal threshold of significance in the meta-analysis (p < 0.001), we find that CpG sites that trend negatively with placental copper are enriched among putative enhancers (OR = 2.15, P = 0.003, N = 813) and amongst regions that flank active transcription start sites (OR = 3.65, P = 1.3 × 10⁻⁴, N = 813, Figure 2b). Accordingly, though not statistically significant, we find those CpG sites are also enriched among CpG shores (1.5 kb out from CpG island regions; OR = 1.4, P = 0.08, N = 813, Figure 2a).

Figure 2. Distribution of results among genomic compartments.

a. Proportion of CpG sites in CpG islands (CGI), CGI shores and CGI shelves. Open sea CpGs are not near any CGI feature. b. Proportion of CpG sites within ChromHMM states. For both figures, compartment distributions are shown for all CpGs on the 450k array, all significant CpG sites (Sig), positively associated CpG sites (Pos), or negatively associated CpG sites (Neg). TSS and Enh are abbreviations of transcription start site and enhancer, respectively.

One of the most significant differentially methylated positions (DMP) in the meta-analysis is found in a CpG island shore and is annotated in foetal placenta cells as being in an enhancer region [25]. Cg25866895 resides in a 5ʹ CpG island shore, upstream of the GSTP1 promoter. It is also in a DNaseI hypersensitivity site, indicating that the CpG is in accessible chromatin [26]. For both cohort-level EWASs, methylation at cg25866895 is negatively associated with copper status, meaning that higher copper associates with lower CpG methylation at the locus (Table S3). Methylation at this position is negatively associated with the expression of GSTP1 (P = 0.005) among 194 placental samples with matched CpG methylation and transcript expression data in the RICHS cohort.

When unaccounted for in the statistical model, population admixture is a confounder of association studies. RICHS is slightly more ethnically diverse than NHBCS (Table 1). Ethnicity was not included in the cohort-level EWASs in this study because of missing race/ethnicity data in both cohorts. However, the inflation factors for the cohort-level results do not suggest population stratification. Additionally, CpGs with significant inter-study heterogeneity (Cochran’s Q-test p-values < 0.05) constituted only 5.5% of tested CpGs, were not enriched among the most significant results (Fisher exact test, Cochran’s Q-test p-value < 0.05 vs. meta p-value < 0.0001) and were removed from downstream analysis. Finally, in a sensitivity analysis conducted on the top significant CpGs in each cohort-level EWAS, the associations of placental copper status and CpG methylation were robust to the addition of race in the regression model (Table S4).

We employed the Comb-p method to find local genomic regions enriched in low P-values from the meta-analyses, considering that these differentially methylated regions (DMRs) might be more biologically relevant than individually significant CpGs. Using the meta-analysis P-values, we found nine significant DMRs (regional P-value < 0.05, Figure 1d, Table 2). Zinc-finger genes are overrepresented in the significant DMRs (4 of 9, Table 2). The most significant DMR is on chromosome 3 at the 5ʹ region of the gene, ZNF197. For the most significant CpG in the DMR at ZNF197, the average difference in per cent methylation between high and low copper levels at ZNF197 is 6.3% (RICHS – 8%, NHBCS – 4.6%, Table S5).

Table 2. Copper-associated differentially methylated regions.

Chr	Start	End	Regional P-value	refGene Symbol	refGene feature^α	In CGI
3	44666362	44667011	5.2E-14	ZNF197*	TSS, 5ʹ UTR	Y
19	52531438	52531957	3.3E-09	ZNF614*	TSS, 5ʹ UTR	Y
11	66624256	66624854	9.4E-09	LRFN4	TSS, 5ʹ UTR	Y
3	44726891	44727579	4.8E-08	ZKSCAN7-AS1	NC intron/exon	Y
6	31783272	31783546	1.1E-07	HSPA1B	TSS, CDS	Y
12	130502003	130502587	6.6E-07	LINC02418	Intergenic	N
8	52722023	52722540	1.5E-06	PXDNL	Intergenic	N
19	37825307	37825680	1.8E-06	HKR1	TSS, 5ʹ UTR	Y
7	94953653	94953811	1.3E-02	PON1	5ʹ UTR, CDS	Y

* Most significant CpG in DMR associates with the transcription of gene.

^α TSS – Transcription start site; UTR – Untranslated region; NC – noncoding; CDS – coding sequence.

Within the RICHS cohort, we assessed 194 placental samples with matched CpG methylation and transcript expression data for functional consequences to differential methylation levels among the significant DMRs. For two of the nine DMRs, methylation level at the most significant CpG within the DMR was associated with lower transcription of the nearest gene (ZNF197 FDR = 6.6e-12; ZNF614 FDR = 2.9e-07, Table 2). Methylation level at the ZNF197 DMR is positively associated with copper status (regional adjusted P < 5.2e-14, Figure 3a,c) and negatively associated with transcription of ZNF197 (FDR = 6.6e-12, Table 2) [26,27].

Figure 3. Copper-associated DMR at ZNF197.

Mean methylation proportion for CpGs within the ZNF197 DMR, between high and low placental copper groups, indicates that higher copper results in higher methylation across the DMR in a. NHBCS and c. RICHS. b. DMR location relative to the ZNF197 start site, and chromatin information for foetal placenta from the ROADMAP epigenomics project [26] (DNaseI hypersensitivity region, RNA-seq signal and ChromHMM annotation) and the UCSC genome browser [27] (CpG island) are also displayed.

The DMR at ZNF197 overlaps the promoter and 5ʹ untranslated region of the ZNF197 gene. Accordingly, chromatin annotation by ChromHMM in foetal placenta cells marks the region as an active promoter or the region flanking an active promoter. The DMR also overlaps a CpG island and DNaseI hypersensitivity site (Figure 3b). Members of the DMR are consistently more methylated in the higher copper status samples, across both cohorts (Figure 3a,c, Table S5).

ZNF197 is known to have a DNA binding transcription factor activity. Specifically, ZNF197 is reported to act as a transcriptional repressor of the master regulator, hypoxia-inducible factor 1 alpha, HIF1α [28]. Given its potential to affect the transcriptional landscape, we sought potential transcriptional targets and cascades of ZNF197. First, we found regions in the genome that matched the consensus binding motif for ZNF197 [29]. We also assessed how the transcript levels of other genes correlate with ZNF197 transcript abundance (co-expression). One hundred and thirty-eight transcripts were significantly co-expressed with ZNF197 (Figure 4). Co-expressed genes are enriched for genes that associate with birth weight in RICHS (OR = 2.9, P = 2.6e-6, N = 194) and genes that are near a ZNF197 consensus binding motif (OR = 1.34, P = 0.01, N = 194). ZNF197 co-expressed genes are enriched for GO biological processes growth hormone secretion (P = 3.4e-4), multicellular organism growth (P = 3.8e-4), and molecular functions related to lipid biosynthesis (P = 1.9e-4; Table S6).

Figure 4. ZNF197 Genome-Wide Co-expression.

Volcano plots display the p-value and effect estimate (log fold change). One hundred and thirty-eight genes are significantly co-expressed with ZNF197 (larger points in the volcano plot). Large grey dots indicate that the association had a false discovery rate (by Benjamini Hochberg method) less than 0.05 and the large black dots indicate a Bonferroni-adjusted p-value less than 0.05.

We used both ZNF197 co-expression and motif analysis to generate a list of potential targets of ZNF197. We narrowed putative binding regions to only those that are found within open chromatin in placenta (ROADMAP Epigenomics, DNaseI hypersensitivity) [26] and were immediately upstream of the transcription start site of a gene annotated in Ensembl (v37) [30]. We narrowed the co-expression results to just those transcripts with P-values less than 0.01 (equivalent to an FDR<0.18). The intersection of genes that were analysed for co-expression with ZNF197 and genes with an accessible ZNF197 binding motif was ranked by P-value, and a combined, optimized super-list of genes were generated by a Cross-Entropy Monte Carlo algorithm [31] (Table 3).

Table 3. Super-list of genes co-expressed with and containing a binding motif for ZNF197.

Gene symbol	Ensembl ID	Motif P-value	Co-expression P-value
GLI2	ENSG00000074047	7.32E-08	1.20E-04
BSG	ENSG00000172270	1.41E-07	4.18E-04
LINC00968	ENSG00000246430	3.07E-07	6.98E-05
ADAMTS9-AS2	ENSG00000241684	4.02E-07	1.02E-05
L3MBTL4	ENSG00000154655	8.68E-10	9.58E-03

Higher placental copper concentration in RICHS is associated with almost 3 times the odds of SGA than AGA while accounting for LGA (multinomial linear regression; N = 141; OR = 2.98; P-value = 0.02). This association was not observed in NHBCS (N = 343; OR = 0.95, P-value = 0.93). While the expression of ZNF197 alone does not associate with birth weight status, the transcripts that are significantly co-expressed with ZNF197 are enriched (OR = 2.5, P = 1.3 × 10⁻¹⁰) among transcripts that are significantly associated with birth weight in a transcriptome-wide association study of birth weight in RICHS (N = 194, FDR<0.05).

Discussion

We combined data from two mother-infant cohorts to evaluate associations between placental copper concentration and genome-wide DNA methylation in placenta. The methods used for quantifying DNA methylation and copper concentration were identical in both cohorts, and the placental concentrations of copper were similar to those reported in other cohorts [32]. On a cohort-level, we did not find any consistent evidence of an association between copper status and methylation status at any single CpG. However, from the meta-analysis, we found evidence of nine larger regions of the genome that were differentially methylated with copper concentration.

In both individual cohorts, and in the meta-analysis – hypomethylation is more common than hypermethylation. This finding is in accordance with previous findings that link increased copper in liver with decreased methylation through the disruption of methionine metabolism [20]. The CpGs that were less methylated in the higher copper placentas were enriched among putative enhancers, regions that flank active transcription start sites and CpG island shores. These regions of the genome represent areas of active DNA methylation turnover. Whereas hypermethylation of normally unmethylated CpG island promoters is a well-described mechanism for transcriptional inhibition, variation in the methylation status of such regions is low and generally a poor predictor of gene expression. Conversely, CpG island shores and enhancer elements exhibit higher variation in DNA methylation, and recent studies stress the importance of altered methylation at enhancer elements as an important contributor to variation in gene expression [33].

GSTP1, one of our most significant DMPs, encodes the dominant isoform of the important detoxifying enzyme, glutathione-S-transferase in the placenta, and is responsible for inactivating reactive oxygen species [34]. In a study of methylation changes in liver associated with Wilson disease, GSTP1 was hypomethylated in both mouse models and Wilson disease patients [21]. Interestingly, while other differentially methylated gene promoters were ameliorated by choline treatment, GSTP1 was not, leading the authors to conclude that differential methylation of GSTP1 was a feature of ATP7B deficiency and copper accumulation, but not the result of the hypomethylation mediated by the copper-induced inhibition of SAHH [21]. It has been suggested by others that upregulated GST expression is a biomarker of cell’s first response to oxidative stress [35].

Our most robust DMR was at the promoter of zinc-finger gene and transcription factor, ZNF197. The DMR overlaps a DNaseI hypersensitivity region, a CpG island and the transcription start site of ZNF197. Copper status is positively correlated with DNAm in the DMR; meaning that higher copper concentration associates with more DNAm. DNAm at the most significant CpG in the DMR was negatively associated with transcript abundance. As the DMR is in a promotor-region, this result is in accordance with the canonical relationship of promoter DNAm and transcriptional repression. ZNF197 is a member of the SCAN family of zinc-finger protein, which is substantially expanded in eutherians (placental mammals) [36,37]. Additionally, ZNF197 specifically appears to be under positive selection among Homininae (includes genera Homo, Pan and Gorilla), and likely plays a role in placentation [38]. Finally, in a hepatic cell line (HepG2), ZNF197 is differentially expressed with cupric chloride treatment [39].

ZNF197 was significantly co-expressed with 138 other gene transcripts, throughout the genome. These genes are enriched for GO biological process terms for growth and growth hormone secretion, which suggests that the ZNF197 transcriptional network could be involved in foetal growth signalling. Further, transcripts that were co-expressed with ZNF197 were enriched among transcripts that are significantly associated with birth weight. We observe a negative association between placental copper and SGA birth in RICHS, but not in NHBCS. RICHS is better powered to detect this association, given the oversampling of extremes, but likely the effect estimate is biased, and would be attenuated in a more general population. Reports of both negative [5] and positive [3,4] associations between placental copper and birth weight exist in extant research studies. This disagreement is not altogether unexpected, since both low and high concentrations of copper are detrimental and the prevalence of copper deficiency varies by the study population. In populations where copper deficiency is common [40,41], one might expect to find a positive correlation between copper and birth weight.

The motif analysis of ZNF197 suggests that it may be involved in copper metabolism, by targeting CCDC93, which is involved in trafficking ATP7A [42]. However, by sharing information between ZNF197 motif-binding and co-expression, we find that targets of ZNF197 may be more broadly involved in placentation, placental development, metabolism and foetal growth. GLI2 is a regulator of the HedgeHog pathway, is required for embryonic development and maintaining placental function [43]. GLI2 is also differentially methylated in placentae from early-onset preeclampsia [44]. BSG, a plasma membrane protein in the immunoglobulin superfamily, is important to implantation [45]. Though it has not been previously described in placenta, in cancer studies, LINC00968 is associated with the activation of the WNT signalling pathway [46]. WNT, a key patterning gene, has been implicated in placental invasion and development through DNA methylation studies [47,48].

Copper concentration in placental tissue and maternal serum has been linked to adverse birth outcomes, including low birth weight [3–5,9], diabetes [49–51] and preeclampsia [11,52]. All of these adverse outcomes are also linked to oxidative stress [6,10,11]. It is intriguing to consider our findings from that perspective. One of our most significant DMPs, GSTP1, is linked to oxidative stress [35,53] and our most significant DMR, ZNF197, has been linked to Hypoxia Inducible Factor 1a (HIF1A) [28], which is stabilized by copper [54]. While our top results, an antioxidant and a transcription factor linked to hypoxia response, may seem paradoxical – hypoxia leads to oxidative stress through the generation of superoxide anions by mitochondria [55]. However, expression of ZNF197 and HIF1A were not significantly associated in this study, and we did not find evidence of differential expression of any other markers of oxidative stress, including the cuproenzyme Cu, Zn-superoxide dismutase (SOD1). Yet, these enzymes are more prominently controlled in the post-translational level, and so not identifying transcript changes may be indicative of their normal physiologic regulation.

To our knowledge, this is the only study yet to examine the relationships between human placental copper and DNAm. Despite this, many of our individual findings are consistent with copper research in other models and tissues. In rodent livers and placentae, higher copper correlates with global hypomethylation [19–21]. Antioxidant enzymes within glutathione metabolism (specifically GSTP1) associate with copper in liver among humans, rodents and canines. Although it has not been studied extensively, ZNF197 is differentially expressed in cultured liver cells treated with copper salts. This consistency across multiple studies, with different technologies, models, and tissues, is notable. It lends support that our results are robust and relevant within the broader topic of placental health and early development. However, our findings should be interpreted within the context of this study’s limitations.

This is an observational study in which placental copper concentrations and DNAm levels were both measured in placenta at term. Thus, we cannot conclude that our results are representative of copper and DNAm associations throughout development. We adjusted for likely confounders in our study, but cannot rule out the possibility that unmeasured or residual confounding remains in our analysis. Two potential confounders are self-reported maternal smoking and ethnicity. We adjusted for any reported smoking during pregnancy in our primary analysis, but self-reported smoking lacks the sensitivity of the biomarker, cotinine, which could also detect exposure to secondhand smoke [56]. Similarly, while sensitivity analysis suggests that self-reported ethnicity had little to no effect in our observations, a genetic estimate of heritage may be more informative. To limit statistical confounding due to differing cellular composition between individuals, our models are adjusted for putative cellular heterogeneity. However, a paucity of information about placental cell composition and methylation patterns requires that we use a reference-free approach [57]. Very few of our EWAS results reached any genome-wide significance threshold. Our DMR and enrichment analyses utilized methylation results with raw p-values < 1 × 10⁻⁴. Thus, there are likely some false-positive findings among these hits. Power is a common problem in epigenetic epidemiology given the high dimensionality of the data and the small effect sizes [58]; however, the accordance with our findings with others is encouraging and sufficient motivation for further study of the role of copper metabolism and methylation in pregnancy and placental function. Further, both cohorts utilized in this study consisted predominantly of healthy white mothers from the New England region of the United States, limiting the generalizability of our findings. Lastly, the placental copper concentrations in this study are similar to those previously reported and are likely within the physiologically normal range. As copper deficiency and toxicity are rare in the US, our results are more representative of general population patterns that result from normal variation in copper metabolism.

Conclusion

Placental and foetal development rely on the balance of copper within the placenta. Here we have shown that DNAm patterns within the placenta associate with placental copper status. Specifically, the methylation level and associated expression level of the antioxidant GSTP1 and the transcription factor gene ZNF197 associate with copper level. The expression of ZNF197 is associated with regulators of placentation, genes that associate with birth weight and are enriched in biological processes including growth hormone secretion and multicellular organism growth. The present study provides evidence that placental copper and copper metabolism play a role in placentation, in part through the differential methylation of GSTP1 and ZNF197.

Methods

Cohorts

The New Hampshire Birth Cohort Study (NHBCS) was initiated in 2009 and is an ongoing study comprised of a cohort of mother-infant pairs. Pregnant women between 18 and 45 years of age were recruited from the study’s participating prenatal care clinics in New Hampshire, USA. Women were included in the cohort if their primary source of drinking water was from an unregulated residential well, they had resided in the same household since their last menstrual period and had no plans to move before delivery. All participants provided written informed consent in accordance with the requirements of the Institutional Review Board (IRB) of Dartmouth College. In this study, NHBCS participants were singleton pregnancies recruited between February 2012 and September 2013. Data provided from this cohort are placental copper concentration and site-specific measures of DNA methylation (n = 306). Questionnaires and medical record abstraction were utilized to collect sociodemographic, lifestyle, and anthropometric data (Table 1).

The Rhode Island Child Health Study (RICHS) enrolled mother-infant pairs from the Women & Infants Hospital in Providence, Rhode Island, between September 2010 and February 2013. Mothers were at least 18 y of age, had no life-threatening conditions, and delivered singletons free of congenital/chromosomal abnormalities at or after 37 weeks of gestation. Infants who were born small for gestational age (≤10th birth weight percentile) or large for gestational age (≥90th birth weight percentile) were oversampled. Infants who were adequate for gestational age (between the 10th and 90th birth weight percentiles), and matched for gestational age and maternal age were coincidentally enrolled. All participants provided written informed consent and all protocols were approved by the IRBs at the Women & Infants Hospital of Rhode Island and Emory University, respectively. Data provided by this study include placental copper concentrations, site-specific measures of placental DNA methylation (Illumina 450 k beadarray), and placental transcript expression (RNA-seq). Interviewer-administered questionnaires were utilized to collect sociodemographic and lifestyle data. Structured medical record review was used to collect anthropometric and medical history data (Table 1).

Data collection

Anthropomorphic measures Both of the studies have parallel demographic and anthropomorphic measures for mothers and babies, as well as placental genome-wide CpG methylation and metal concentration data. Birth weight (g), head circumference (cm), and birth length (cm) were abstracted from medical records. z-Scores were calculated for each, standardized by gestational time and infant sex, via Fenton growth curves [59]. Infants with standardized BW percentiles below the 10th percentile are classified as small for gestational age (SGA); those above the 90th percentile as large for gestational age (LGA); and those between the 10th and 90th percentile as adequate for gestational age (AGA). Mother’s BMI was calculated from height and weight collected from medical records (Table 1).

Study covariates Self-reported highest achieved maternal education level was coded as either high school or less or more than high school. Maternal smoking was defined as any self-reported smoking during any point of pregnancy versus reported no smoking during pregnancy. Foetal sex was obtained from medical records. Foetal race was self-reported and coded as white or not white.

Sample collection The same placental sampling protocol was used by both cohorts. Placenta samples were biopsied from the foetal side adjacent to the cord insertion site within 2 h of delivery. Maternal decidua was removed. Samples were placed in RNAlater (Life Technologies) and then frozen at −80°C. Both RNA and DNA were extracted (Norgen Biotek) and quantified using a Qubit Fluorometer (Life Technologies) and were subsequently stored at −80°C.

Measurements

CpG methylation Both cohorts measured DNA methylation via Illumina Infinium Human Methylation 450 K BeadArray (Illumina) at the University of Minnesota Genomics Centre (NHBCS n = 306, RICHS n = 141). Bisulphite modification was performed with the EZ Methylation kit (Zymo Research, Irvine, CA). Samples were randomized across multiple batches and batch variables were recorded to allow for batch-effect corrections. Data were assembled using BeadStudio (Illumina). The raw methylation array data are available via the NCBI Gene Expression Omnibus (GEO) for NHBCS and RICHS via accession numbers GSE71678 and GSE75248, respectively.

Quality control, normalization and filtering of methylation data were all performed in the R statistical programming language. We excluded CpG probes that: had poor detection p-values (p-value > 0.001). No samples had to be excluded due to poor quality (as indicated by having more than 5% failing CpG probes). Background correction, dye bias, and functional normalization were performed via the minfi Bioconductor package [60–62]. Methylation data were exported from minfi as β-values, representing the proportion of methylated alleles for each CpG site. Standardization across probe types was performed via beta mixture quantile (BMIQ) normalization [63]. β-values were adjusted for batch effects using the combat function in the R package, sva [64,65]. β-values were logit-transformed (M-values) to better approximate a normal distribution [66]. Prior to downstream analysis, CpG probes were filtered to remove probes that were located on the X- or Y-chromosomes, had a single nucleotide polymorphism (SNP; with minor allele frequency >1%) within 10 bp of the target CpG [67], or were non-specific [68]. After QC and filtering, 409,882 CpG probes from 306 placental samples for NHBCS and 398,332 CpG probes from 141 placental samples for RICHS remained for analysis.

Using the 450 K CpG methylation data, we estimated the proportions of constituent putative cell types in each placental sample with non-negative matrix factorization via RefFreeEWAS package in R [57]. Using the 10,000 CpGs with the greatest variability in methylation, we identified four putative constituent cell mixtures in the RICHS and NHBCS cohorts, respectively, and then used the full set of CpGs for each cohort (NHBCS p = 408,367 and RICHS p = 397,040) to estimate the relative proportions of each putative cell type per placental sample.

Gene expression Transcriptome-wide sequencing was performed on 197 placental samples from RICHS, including the 141 samples with 450 K methylation array data and measured placental copper. We used an RNeasy Mini Kit (Qiagen) to isolate total RNA, which was then stored at −80°C until analysis. We quantified RNA using a Nanodrop Spectrophoto-metre (Thermo Scientific), assessed RNA integrity via Agilent Bioanalyzer (Agilent), removed ribosomal RNA with a Ribo-Zero Kit [69], converted to cDNA using random hexamers (Thermo Scientific), and performed transcriptome-wide RNA sequencing via the HiSeq 2500 platform (Illumina) [70]. Raw reads have been deposited at the NCBI sequence read archive (SRP095910).

QC was performed in FastQC. Reads passing QC were mapped to the human reference genome (h19) using the Spliced Transcripts Alignment to a Reference (STAR) aligner. Further QC was performed with the Bioconductor package, edgeR [71]. Low-expressed transcripts were excluded and the remaining transcripts normalized via the trimmed mean of M-values (TMM) method [72]. Robust estimates of dispersion were calculated [73]. Annotation for the transcripts was provided by the biomaRt package for the ensembl GCh37 human reference genome [30,74]. Transcripts that had no annotation were removed. Normalized transcript data were exported from edgeR as counts per million (CPM) for surrogate variable analysis.

To adjust for unmeasured confounding effects, such as cell-type heterogeneity, we generated surrogate variables from the normalized transcript reads (as CPM) via the svaseq function in the sva package [65,75]. The number of surrogate variables was estimated with the num.sv function using the Buja and Eyuboglu (‘be’) method and 100 permutations. In the svaseq function, the iteratively re-weighted least-squares method was used to empirically estimate of control transcripts. The model (design) argument used for num.sv and svaseq differed by the application. To compare methylation and expression, target CpG M-values alone were used for the model. For co-expression analyses, only log₂(CPM) of the target transcript was used for the model. The recommended number of surrogate variables were added to downstream models.

Copper quantification For both RICHS and NHBCS placental samples, copper concentrations were quantified at the Dartmouth Trace Elements Analysis Core using inductively coupled plasma mass spectrometry (ICP-MS). Details of the processing are described elsewhere [76]. None of the samples in either cohort had non-detectable copper concentrations. We dichotomized copper concentrations by the study-specific median concentration (857.2 μg/g) for the EWAS analysis.

Statistical analysis

Robust linear modelling of methylation data For each CpG, we used robust linear models to regress methylation M-values on copper status for all samples. Each model included covariates for infant’s sex, maternal smoking status, maternal education level and the four cell-type covariates determined by refFreeEWAS. This analysis was conducted in R. The robust linear models were run with the rlm function from the MASS package [77]. Estimates, standard errors, test statistics, and P-values were generated with the coeftest function in the lmtest package [78]. The covariance matrix of the estimated coefficients (vcov argument) were supplied by the vcovHC function (type = ‘HC0’) in the sandwich package [79,80].

METAL meta-analysis We used inverse variance weighted fixed-effects models to meta-analyse the cohort-specific copper-associations via the METAL software package. We filtered CpGs with significant interstudy heterogeneity (Cochran’s Q-test p-values < 0.05) [81]. In R, the p.adjust function in the stats package was used to calculate false discovery rates [82].

Comb-p for DMRs To find contiguous genomic regions with differential methylation (differentially methylated region, DMR), we used the Comb-P method [83]. Comb-P identifies genomic regions enriched for low p-values (meta P < 1e-4), corrects for auto-correlation with neighbouring CpGs (within 1000 bp) using the Stouffer-Liptak method, then generates a regional P-values adjusted for multiple testing, using the Sidak correction. DMRs with an adjusted regional P-value less than 0.05 were considered significant.

CpG Annotation Genomic intervals for CpG islands (CGI), chromHMM states [25,26], and DNaseI hypersensitivity were downloaded from the UCSC Genome Browser (hg19/GRCh37) [27]. CpG sites were annotated to a CGI if they were within a CGI interval, annotated to the CGI shores if they were within the 2 kb region flanking a CGI, and annotated to a CGI shelf if they were within the 2 kb flanking CGI shores. Open sea CpGs are not within CGI, shores or shelves. CpGs were also annotated to a chromatin-based functional genomic category (e.g. promoter, enhancer, transcribed region) using ChromHMM, established by Ernst and Kellis [25] and applied to foetal placenta cells [26]. DNaseI hypersensitivity intervals were used to confirm open chromatin around significant DMRs.

edgeR-robust QL F-tests of sequence transcripts To model RNA-seq data, we used the R Bioconductor package, edgeR. Normalized and filtered data were fit to a negative binomial-generalized linear model, with robust quasi-likelihood (QL) methods to account for gene-specific variability (glmQLFit, robust). Hypothesis testing via edgeR quasi-likelihood F-tests (glmQLFtest) compared methylation M-values from the most significant CpG within each DMR to the log(CPM) of transcripts within 1 MB of the DMR, while also accounting for the recommended number of surrogate variables [73]. Transcripts with an FDR < 0.05 were considered significantly differentially expressed.

edgeR was used as described above to assess the co-expression of genes. For this application, log(CPM) for ZNF197 (ENSG00000186448) was tested against the log(CPM) of all other transcripts passing QC and filtering, while also accounting for the recommended number of surrogate variables. Transcripts with an FDR < 0.05 were considered significantly differentially co-expressed with ZNF197.

edgeR was used as described above to model the associations between birth weight percentile and RNA transcript abundance. Here, each infant’s birth weight percentile [59] was tested against the log(CPM) of all other transcripts passing QC and filtering, while also accounting for race, sex, known batch and the recommended number of surrogate variables. Transcripts with an FDR < 0.05 were considered significantly differentially co-expressed with ZNF197.

MEME-FIMO motif analysis The consensus binding motif for ZNF197 [29] was compared the Ensembl human reference genome (GRCh38) in the Find Individual Motif Occurrences (FIMO) tool from MEME suite 5.0.2 (http://meme-suite.org; [84]). The results were re-annotated to the GRCh37 Ensembl human reference genome and then filtered to only those occurrences that coincided with open chromatin as indicated by the DNase-seq of a placenta sample from a 113d foetus (Roadmap Epigenomics Project; ENCSR035QHH [26]). Motifs were also filtered to those that matched orientation with and were within 1500 bp of a gene TSS (in Ensembl GRCh37). Motifs were ranked by the P-values assigned by FIMO.

Rank aggregation Ensembl IDs and (minimum, if multiple IDs were present) p-values from the ZNF197 co-expression analysis and motif analysis, respectively, were filtered to the intersecting Ensembl IDs and ranked by P-value. The RankAggreg function from the RankAggreg R package [31] was used to create an optimized rank of the top 25 ensembl IDs, using the Cross-Entropy Monte Carlo algorithm, (P-value) weighted-Spearman distance, and a rho value of 0.01.

Enrichment tests Fisher’s exact tests were used to test the enrichment of significant ZNF197 motif occurrences among genes that are significantly co-expressed with ZNF197. Fisher’s exact tests were also used to test the enrichment of genes that are co-expressed with ZNF197 among the genes that are associated with birth weight percentile. Finally, among significant CpGs from the meta-analysis, fisher’s exact tests were used to assess enrichment of positively or negatively associated CpG sites within each of the ChromHMM chromatin states, as well as within CGI, shores and shelves.

Gene Ontology analysis Gene ontology analysis was conducted via the goana function in the limma package [85]. GO analysis was used to test for the enrichment of gene ontology terms among the genes that are co-expressed with ZNF197.

Multinomial logistic regression Multinomial logistic regression was used to determine the odds of SGA or LGA, with reference to AGA status with placental copper concentration. Birth weight group (defined above) was regressed on dichotomized placental copper concentration, with the function multinom function in the R package nnet [77].

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the National Institutes of Health under Grants NIH-NIEHS R24 ES028507, NIH-NIEHS R01 ES022223, NIH-NIEHS R21 ES028226, NIH-NIGMS P20 GM104416, and NIH-NIEHS P01 ES022832; and the U.S. Environmental Protection Agency (EPA) under U.S. EPA grant RD83544201. Its contents are solely the responsibility of the grantee and do not necessarily represent the official views of the U.S. EPA. Further, the U.S. EPA does not endorse the purchase of any commercial products or services mentioned in the presentation.