Prenatal exposure to endocrine disrupting chemicals is associated with altered DNA methylation in cord blood

ABSTRACT

Prenatal exposure to endocrine disrupting chemicals can interfere with development, and has been associated with social-cognitive functioning and adverse health outcomes later in life. Exposure-associated changes of DNA methylation (DNAm) patterns have been suggested as a possible mediator of this relationship. This study investigated whether prenatal low-dose exposure to polychlorinated biphenyls (PCBs) and polychlorinated dibenzo-p-dioxins and dibenzofurans (PCDD/Fs) is associated with altered DNAm patterns across the genome in a Western urban-industrial population. In 142 mother-infant pairs from the Duisburg Birth Cohort Study, PCBs and PCDD/Fs levels were quantified from maternal blood during late pregnancy and associated with DNAm levels in cord blood using the Illumina EPIC beadchip. The epigenome-wide association studies (EWAS) identified 32 significantly differentially methylated positions (DMPs) and eight differentially methylated regions (DMRs) associated with six congeners of PCB and PCDD in females or males (FDRs < 0.05). DMPs and DMRs mapped to genes involved in neurodevelopment, gene regulation, and immune functioning. Weighted gene correlation network analysis (WGCNA) showed 31 co-methylated modules (FDRs < 0.05) associated with one congener of PCDF levels in females. Results of both analytical strategies indicate that prenatal exposure to PCBs and PCDD/Fs is associated with altered DNAm of genes involved in neurodevelopment, gene expression and immune functioning. DNAm and gene expression levels of several of these genes were previously associated with EDC exposure in rodent models. Follow-up studies will clarify whether these epigenetic changes might contribute to the origin for adverse mental and health outcomes.

KEYWORDS:

• Endocrine disrupting chemicals
• PCB
• PCDD/F
• EWAS
• DNA methylation
• WGCNA
• cord-blood
• Duisburg birth cohort

Background

The developing foetus is highly sensitive to perturbations in the internal and external environment [1]. Foetal exposure to a range of environmental influences, including maternal subnutrition [2,3], maternal psychosocial stress [4,5], or exposure to chemicals [6] has been associated with lifelong consequences for the developing organism [7–9]. There is a large body of literature investigating the long-term effects of prenatal exposure to endocrine disrupting chemicals (EDCs) on physiological and psychological functioning [10]. EDCs, including polychlorinated biphenyls (PCBs), polychlorinated dibenzo-p-dioxins and dibenzofurans (PCDD/Fs), are widely dispersed, environmentally persistent organic pollutants (POPs) that can act as exogenous agents interfering with endocrine processes [11]. PCDD/Fs, also known as dioxins, and PCBs are complex mixtures of congeners of polychlorinated organic compounds that share a common molecule origin and structure but vary in the number and location of chlorine atoms introducing differential functioning of each congener [12,13]. Since 1929, PCBs have been synthesized and used as synthetic oils for various industrial deployments such as in surface coating or in electrical apparatus as dielectric and coolant fluids until the late 1970s. PCDD/Fs are inadvertent by-products of synthetic processes such as smelting, chlorine bleaching of paper, or manufacturing of several herbicides and pesticides. Due to their chemical lipophilic characteristics and a generally long half-life, PCBs and PCDD/Fs persist in widely geographically dispersed environments over long periods of time and are able to bioaccumulate in the fatty tissue of humans and wildlife [11].

The exposure to EDCs has decreased in the general population over the last years [14]. Nevertheless, exposure to PCBs and PCDD/Fs below toxic level – so called low-dose exposure – is still of great concern, since exposure to even low doses can be harmful for the developing organism [15]. There is empirical evidence for adverse effects of prenatal low-dose exposure to PCBs and PCDD/Fs on cognitive [16–18] and motor development [18–20], pubertal development [21,22] as well as gender-typical play behaviour [23], attention deficits [24,25] and autistic behaviour [26,27].

The observation of such long-term influences of foetal exposures has led to the hypothesis that epigenetic mechanisms might explain how these effects retain their stability [28]. The most extensively studied epigenetic mark in humans is DNA methylation (DNAm), a chemical modification of the DNA with gene regulatory function [29]. During foetal development, DNAm is involved in the orchestration of developmental programming processes and is essential for normal cellular development and differentiation [30]. Perturbations of DNAm processes through exogenous factors such as chemicals are thought to underlie structural and functional changes in cells, tissues, and organ systems, increasing disease risk across the lifespan (developmental origin of health and disease; DOHaD) [7,31,32].

So far, it has rarely been investigated in humans whether prenatal exposure to PCBs or PCDD/Fs might be associated with altered DNAm [28,33,34]. First empirical insights into this question come from animal models that have shown global DNAm changes after prenatal exposure to PCBs or PCDD/Fs [e.g., 35,36]. Findings from animal models have been partly confirmed in humans by very few studies, reporting associations between exposure to PCBs or PCDD/Fs and global DNAm, which was estimated via the Alu and LINE1 assay [37,38].

The first epigenome-wide association study (EWAS) of prenatal exposure to various PCB congeners was realized in a cohort of Faroes individuals (N = 72) using the Infinium Methylation 450k chip [33]. These individuals were exposed to high levels of PCBs, much higher than other Western populations, due to Faroes traditional dietary habits including whale meat and blubber. Leung et al. [33] showed that higher doses of single PCB congeners were significantly associated with DNAm changes in cord blood in males and females. Among the investigated PCB congeners, mainly PCB congener 105 (PCB #105) was significantly associated with differential DNAm in sex-combined analysis. Notably, in a sex-stratified analysis, DNAm changes were significantly associated with PCB #105 exposure in females only, whereas in males, literally no single PCB congener remained significantly associated with DNAm changes.

Lastly, studies investigating other EDCs showed gene-specific or global DNAm changes after prenatal exposure to EDCs. For instance, significant associations were found between phthalates and gene-specific DNAm changes in Mexican-American newborns [39] as well as between bisphenol A (BPA) and gene-specific DNAm changes in the Michigan Mother-Infant Pairs birth cohort [40]. Also, prenatal exposure to phthalates revealed significant associations with epigenome-wide DNAm marks in placental tissues [41].

Taken together, since prenatal low-dose exposure to EDCs including PCBs and PCDD/Fs seems to be especially harmful for the developing organism [11], and has been associated with developmental sequelae, we aimed to explore whether prenatal low-dose exposure to PCBs and PCDD/Fs is associated with specific DNAm patterns in a Western urban-industrial population. In the Duisburg Birth Cohort, where maternal blood levels of PCBs and PCDD/Fs were objectively quantified during late pregnancy, DNAm in cord blood was quantified at around 850,000 sites across the genome with the Illumina Methylation EPIC beadchip. Based on previous results reporting sex-specific associations between EDC exposure and DNAm patterns, we hypothesized sex-specific effects and conducted sex stratified analyses.

Methods

Study area and participants

The Duisburg Birth Cohort study was conducted in Duisburg, Germany after one ton of dioxinated dust leaked from a recycling firm in 1999. Duisburg belongs to the Ruhr District, an industrial conurbation in Western Germany. A detailed description of the study area and study population is given in Wilhelm et al. [42]. Initially, 232 healthy mothers from Duisburg aged between 18 and 42 years who became pregnant after the leakage were recruited between 2000 and 2002. During the third trimester of pregnancy, the concentration levels of PCBs and PCDD/Fs (reflecting cumulative exposure from various potential EDC sources) were assessed from maternal blood. The mother´s children have been followed from birth until today almost annually. Healthy mother-child pairs were included if they met the following criteria: German or Turkish as a first language, born at term (weeks 38–42 of pregnancy), with an APGAR score of at least 8, of parity 1–3, without serious complications or illnesses during pregnancy or parturition and without congenital anomalies. Since 2007, regular follow-up studies have been conducted to examine child development. See additional file 1 supplementary table (ST) 1 and 2 for descriptive data on neonate, mother and pregnancy characteristics.

Ethics approval and consent to participate

The Duisburg Birth Cohort study was approved by the Ethics Committee of the Medical Faculty of the Ruhr-University Bochum (registry no. 1478); follow-up was approved by the Ethics Committees of the Faculty of Psychology of the Ruhr-University Bochum in 2016 (No. 20160126). The study was conducted in accordance with the ethical principles for medical research involving human subjects as defined by the Declaration of Helsinki. All parents gave written informed consent.

Collection and analysis of maternal blood samples

To assess the concentration levels of the congeners of PCB and PCDD/F, maternal blood samples (50 mL, N = 226) were taken between the 28th and 42nd week of pregnancy (median: 32 weeks). In seven cases, blood samples were taken during the first 38 days after birth. Blood samples were stored at −80°C until analysis. Analysis of PCBs and PCDD/Fs in maternal blood was performed at the Department of Hygiene of the Ruhr-University Bochum by standard procedure with capillary gas chromatography and high-resolution mass spectrometry described in Wittsiepe et al. [43]. Measured concentrations of single congeners of PCBs and PCDD/Fs were relativized on 0.5 limit of detection (<1 pg/g on lipid base). The measured levels of PCBs and PCDD/Fs have been classified as low-dose by comparison to reference levels of the human biomonitoring data from Germany matching sampling years and participant’s age [43].

Sixteen congeners of PCBs and PCDD/Fs with concentrations higher than the limit of detection in at least 75% of the samples were included in the data analysis (congener list in additional file 1 ST1). Concentration of single congeners of PCBs and PCDD/Fs were measured in natural logarithm transformed (nat.log.) pg/g lipid base.

DNA isolation and genome-wide methylation assessment

Genomic DNA was isolated from umbilical cord blood using the QIAamp DNA Blood Mini Kit (Qiagen, Hilden; Germany). Cord blood was taken directly after birth in the obstetrical clinic and stored at −80°C until DNA extraction and further analysis. DNA bisulphite conversion and assessment of genome-wide DNAm levels at 850,000 CpG sites using the Infinium HumanMethylationEPIC BeadChip array (Illumina, San Diego, CA, USA) were performed according to the manufacturer’s protocol at the Life&Brain Core Facility (University of Bonn, Germany).

Data analysis

All data analyses were performed with R [44] using RStudio 3.6.1 (except for WGCNA using RStudio 3.3.4). The supporting R-code of all data analyses is available at https://github.com/kattamatta/duisburgbirthcohortstudy.

DNA methylation data preprocessing

Raw Infinium EPIC data was preprocessed using the pipeline of the RnBeads package [45,46] for intensity summarization and calculation of methylation ratio (beta values). The preprocessing steps included primary quality control as well as a filtering step before and after normalization. Seven samples with low bisulphite conversion efficiency compared to control probes were removed, as well as one of two samples with identical SNP methylation patterns. Moreover, six samples that did not match a priori criteria (Caucasian ethnicity as stated by the participant’s country of origin, no missing data regarding congers and sex of neonate) were excluded. During the first filtering step, cross-reactive and SNP-enriched probes overlapping with at least three bases in the target sequence with annotated SNPs were eliminated. Low-quality probes were removed using RnBeads´ Greedycut algorithm with default settings. Normalization was conducted using the dasen method from the wateRmelon R-package [47] implemented in RnBeads. For the second filtering step, non-CpG probes (i.e., CC, CAG, CAH, CTG and CTH) and CpGs containing missing beta values were excluded. CpGs on the X chromosome for male and female neonates and Y chromosomes for male neonates were not filtered to conduct sex-stratified analyses, since recent research [13,33] has indicated that epigenetic marks due to prenatal exposure to EDCs might be sex-dependent. The processing procedure lead to a final set of 760,403 CpGs for female samples and 760,464 CpGs for male samples. Final quality control revealed three samples with distinct genome-wide DNAm patterns using multidimensional scaling plots of M values (related to beta values by logit transformation) and three samples mismatching known and predicted sex and were removed. This led to a total of 142 samples (73 females, 69 males).

Assessment of cell-type heterogeneity

To account for cell-type heterogeneity in the DNAm profiles, the reference-based approach by Houseman et al. [48] implemented in RnBeads was used with default settings to estimate cell-type compositions from the assessed DNAm profiles. There has been only one reference panel available for cord blood assessed with the EPIC Infinium array from 83 Asian neonates including CD4T, CD8T, natural killer cells, B cells, granulocytes, and monocytes (nucleated red blood cells were not available) [49]. Thus, DNAm profiles of this reference panel were used and processed together with the target DNAm profiles to estimate the relative contribution of six cell types in the target profiles. The six estimated cell-type proportions of CD4T, CD8T, natural killer cells, B cells, granulocytes, and monocytes were included as covariates in the inferential analyses.

Epigenome-Wide Association Studies (EWAS)

A series of sex-stratified EWAS was conducted to identify DNAm sites and regions significantly associated with the single PCB and PCDD/F congeners. In a first set of EWAS, differentially methylated positions (DMPs) were assessed using limma [50] by fitting multiple regression models in combination with outlier-robust empirical Bayes moderation for each retrieved CpG in M values as criterion and single PCB and PCDD/F congeners (nat.log. pg/g lipid base in maternal blood) as continuous predictor. For each congener an EWAS was conducted, since in general individual congeners are discussed to differ in the effects on health [13]. Following this argument assuming independency of congener’s effects, the single congener might also have differential effects on DNAm. Additionally, a set of possible confounders (smoking status of mother during pregnancy [yes or no], age of mother at birth [years], gestational age [weeks], weight of neonate [g]), technical batch (EPIC assay plate [plate 1 or plate 2]) and the six estimated cell type contributions to account for cell-type heterogeneity were included in each regression model simultaneously. Since PCBs and PCDD/Fs assessed during pregnancy are generally known to be associated with pregnancy parameters, such as gestational length [51] and birth weight [52], interaction terms, were included in a model, when the regarding predictor was significantly (p < 0.01, additional file 1 ST3-ST5) correlated with an above listed covariate. For each regression model separately, outlier CpGs indicating hemimethylated sites were excluded to avoid a bias regarding false-positive associations due to outliers [53].

Following EWAS that identified DMPs significantly associated with a particular congener of PCB or PCDD/F for the given sex, a second set of EWAS was implemented to identify differentially methylated regions (DMRs) associated with that same single PCB and PCDD/F congener (nat.log. pg/g lipid base from maternal blood). Thus, a separate model for each congener of PCB or PCDD/F significantly associated with a DMP in the given sex was conducted for this sex including the same confounders as in DMPs-models to identify DMRs. Each model to identify DMRs was separately implemented using dmrCate [54] with default parameters. Significant DMRs were defined by a minimum of two significant CpGs in a genomic bin smaller than 1,000 nucleotides.

To gain insight into possible over-represented gene pathways in the identified significant DMPs and DMRs, gene ontology (GO) enrichment analyses were performed for these DMPs and DMRs in a sex-stratified manner using the R-package missMethyl [55]. By taking into account all tested CpGs to correct for biases introduced due to differing numbers of probes per gene present on the EPIC array, the significant DMPs and DMRs were tested for GO term enrichment.

Weighted gene correlation network analysis (WGCNA)

To undertake a systems-level view of the DNAm data, weighted gene correlation network analysis (WGCNA) [56] was applied to the top 50% of the ranked by median absolute deviation (MAD) preprocessed CpGs in a sex-stratified manner. WGCNA offers an alternative clustering strategy to the grouping approach of DMRs by employing network analysis to the beta values of CpGs. Using the R-package WGCNA [56], scale-free correlation-based hierarchical clusters of the methylation profiles of CpGs were used to construct weighted co-methylation modules. These modules represent discrete networks of co-methylated CpGs. An unsigned network with outlier-robust bidweight midcorrelation and a soft-threshold for power adjacency transformation was chosen. For reasons of computational time, block-wise module detection of pre-clustered blocks and soft-threshold was applied on the preprocessed beta-values of the top 50% of CpGs ranked by MAD, including those on sex chromosomes for samples of each sex separately (female: 380,202 CpGs; male: 380,232 CpGs). The default procedure was used to estimate the soft-threshold power (female: ß = 5, male: ß = 2) for network construction based on scale-free topology analysis, followed by the topological overlap approach to calculate a dissimilarity between MAD filtered CpGs to reflect interconnectedness. The approach of average linkage hierarchical clustering was then used to divide CpGs with high topological overlap into modules. Modules were summarized by intramodular connectivity as well as modules EigenCpGs. The intramodular connectivity was assessed by the so-called hub-CpG, which represents the most highly connected CpG within a module. To assess whole-module behaviour, a single EigenCpG reflecting the first principal component of variance was calculated for each module separately and for distinction between EigenCpG labelled by colours. These EigenCpGs were linearly regressed onto the single PCB and PCDD/F congeners in a sex-stratified manner, including the same confounders used in the EWAS described above. GO enrichment analyses using the R-package missMethyl [55] were performed for significant modules´ CpGs, to gain insight into possible over-represented gene pathways using the same approach used in EWAS described above.

For all analyses of EWAS and WGCNA, p-values were adjusted for multiple testing by the Benjamini-Hochberg approach and the p-value cut-off was determined by a FDR less than 0.05.

Results

Descriptive data of participants, PCBs and PCDD/Fs

One hundred and forty-two mother-infant pairs were included in all data analyses; 73 (51.41%) neonates were female and 69 (48.59%) were male (additional file 1 ST1 & ST2). Levels of concentration for the individual congeners of PCBs and PCDD/Fs did not significantly differ between sexes (p > 0.05; additional file 1 ST6 & ST7). Their average birth weight was M = 3408.98 g (SD = 569.50 g). On average, the mothers were M = 31.87 years (SD = 4.67 years) old at giving birth and their BMI before pregnancy was M = 24.02 kg/m² (SD = 4.74 kg/m²). Thirty-eight (26.76%) mothers indicated that they had consumed at least one cigarette during the pregnancy. Mean gestation length was M = 39.39 weeks (SD = 1.60 weeks). Fifty (35.21%) mothers left school with University Entrance Exam (German Abitur), 13 (9.15%) with an advanced technical college certificate, 49 (34.51%) with a secondary school certificate, 25 (17.61%) with basic school qualification and five (3.52%) mothers had no or another school-leaving qualification.

Epigenome-Wide Association Studies (EWAS)

Differentially methylated positions (DMPs)

Results of the limma models showed that DNAm of 32 CpG sites were significantly associated with levels of six of the 16 PCB and PCDD/F congener as continuous predictors at FDR-corrected threshold of less than 0.05 (Figure 1, Table 1(A, B); additional file 1 ST8 & ST9). The various significant DMPs differed between female and male neonates and were found for PCB and PCDD but not PCDF congener predictors (Figure 2).

Table 1. Differentially methylated positions associated with PCB and PCDD congeners in A) female or B) male subsample.

A) Female subsample.
Congener	DMP	Chromosome: position (strand)	UCSC RefGene	log2FC (CI 95%)	FDR adj. p
PCB169	cg08931983	Chr3:97,620,931 (-)	CRYBG3	1.97 (1.08, 2.87)	0.03
PCDD67	cg22358250	Chr4:26,323,220 (+)	RBPJ	5.20 (3.94, 6.45)	0.01
	cg07612923	Chr3:117,604,196 (+)	–	6.69 (5.63, 7.75)	0.02
B) Male subsample.
Congener	DMP	Chromosome: position (strand)	UCSC RefGene	log2FC (CI 95%)	FDR adj. p
PCB126	cg02887499	Chr6: 26,158,339 (-)	H2BC5	5.74 (4.53, 6.96)	0.003
	cg04423493	Chr14:97,845,775 (+)	–	4.24 (3.56, 4.91)	0.003
	cg05889395	Chr1:231,980,899 (-)	DISC1	3.65 (2.67, 4.62)	0.003
	cg14574369	Chr12:57,483,631 (-)	NAB2	6.04 (4.19, 7.90)	0.003
PCDD66	cg02887499	Chr6:26,158,339 (-)	H2BC5	5.66 (5.04, 6.29)	6.58e-06
	cg05889395	Chr1:231,980,899 (-)	DISC1	3.19 (1.55, 4.83)	0.0006
	cg17130815	Chr7:2,558,891 (-)	LFNG	3.17 (2.46, 3.88)	0.0006
	cg23894287	Chr8:38,794,831 (-)	PLEKHA2	5.02 (4.51, 5.53)	0.0006
	cg07143052	ChrX:53449486 (-)	SMC1A; RIBC1	3.51 (2.83, 4.20)	0.002
	cg27285278	Chr5:180,650,776 (-)	TRIM41	3.23 (2.74, 3.71)	0.002
	cg15247269	Chr6:31,795,125 (+)	HSPA1B	3.75 (2.90, 4.59)	0.002
	cg07135157	Chr16:3,019,062 (-)	PAQR4	1.67 (0.58, 2.76)	0.003
	cg18883807	Chr14:96,493,999 (-)	–	2.60 (1.56, 3.63)	0.003
	cg20280727	Chr6:125,474,808 (-)	TPD52L1	3.08 (2.19, 3.97)	0.01
	cg05849676	Chr12:10,171,151 (+)	CLEC12B	−1.94 (−2.63, −1.25)	0.01
	cg24051516	Chr11:58,507,237 (+)	–	1.68 (0.86, 2.51)	0.01
	cg11550825	Chr11:62,013,092 (+)	–	2.54 (1.93, 3.14)	0.02
	cg21981270	Chr10:45,914,525 (-)	ALOX5	3.71 (3.35, 4.07)	0.02
	cg04423493	Chr14:97,845,775 (+)	–	3.49 (2.93, 4.04)	0.02
	cg16127683	Chr15:40,268,777 (-)	EIF2AK4	−4.99 (−5.78, −4.21)	0.02
	cg10295718	Chr10:86,768,400 (-)	–	−1.93 (−2.42, −1.44)	0.02
	cg15451272	Chr1:25,071,355 (-)	CLIC4	2.94 (1.92, 3.95)	0.03
	cg09173924	Chr17:73,513,338 (-)	TSEN54	3.54 (2.67, 4.41)	0.03
	cg02915720	Chr7:99,558,185 (+)	–	−2.19 (−3.19, −1.19)	0.03
	cg24735210	Chr7:38,351,468 (+)	–	2.01 (1.28, 2.75)	0.03
	cg14036884	Chr19:7,694,678 (-)	XAB2; LOC100131801	3.70 (2.88, 4.53)	0.03
	cg02846961	Chr19:812,000 (+)	PTBP1	−1.69 (−2.67, −0.72)	0.04
	cg06111627	Chr3:127,309,564 (+)	TPRA1	1.68 (1.07, 2.28)	0.04
	cg23154011	Chr1:23,340,513 (+)	C1orf234	−1.28 (−2.12, −0.44)	0.04
	cg20973769	Chr10:72,494,617 (-)	ADAMTS14	−1.41 (−2.44, −0.39)	0.04
	cg12193927	Chr5:140,571,783 (+)	PCDHB10	3.72 (3.20, 4.24)	0.048
	cg15869663	Chr3:151,986,198 (+)	MBNL1-AS1; MBNL1	1.63 (0.90, 2.36)	0.049
PCDD70	cg02887499	Chr6:26,158,339 (-)	H2BC5	5.49 (4.57, 6.42)	0.02
PCDD75	cg17130815	Chr7:2,558,891 (-)	LFNG	2.70 (1.96, 3.45)	0.04

Note: DMP = differentially methylated position, UCSC RefGene = annotation by UCSC Gene, log2FC = estimated log2-fold-change, CI 95% = Confidence interval of 95%, FDR adj. p = FDR adjusted p-value, – = no annotation.

Figure 1. Stacked manhattan plot of DMPs from six EWAS for each sex. Stacked manhattan plot of results from the six EWAS for each sex that indicates significant DMPs. -log10 transformed FDR of each CpG included in the EWAS are plotted in order of its location on the chromosome. The dashed line indicates the array-wide threshold of -log10(0.05). Each DMP is labelled by the associated congener predictor and sex, i.e., which congener was significantly associated with DNAm in the male or female subsample. DMPs that are associated with more than one congener are additionally labelled with the gene the DMP maps to (genomic location in Table 1(A, B)).

Figure 2. Significant DMPs common between different congener predictors. Euler diagram of the number of significant DMPs that were common between different congener predictors and A) males or B) females.

In the female sample, a total of three significant DMPs were identified (Table 1A; additional file 1 ST9). An increasing level of PCB169 was associated with increased DNAm of cg08931983 (FDR = 0.03). An increase of the PCDD67-level was associated with two DMPs (cg22358250, FDR = 0.01; cg07612923, FDR = 0.02).

In males, four out of 29 DMPs were significantly associated with levels of more than one congener (Table 1B; additional file 1 ST9). Cg02887499 annotated for gene H2B Clustered Histone 5 (H2BC5) was positively associated with PCB126 (FDR = 0.003), PCDD66 (FDR = 6.58×10⁻⁶) as well as PCDD70 (FDR = 0.02). In addition, cg05889395, annotated for gene Disrupted In Schizophrenia 1 (DISC1) was positively associated with PCB126 (FDR = 0.003) and PCDD66 (FDR = 0.0006). Cg04423493, not annotated for a gene, was positively associated with PCB126 (FDR = 0.003) and PCDD66 (FDR = 0.02). Cg17130815, annotated for LFNG O-Fucosylpeptide 3-Beta-N-Acetylglucosaminyltransferase (LFNG), was positively associated with PCCD75 (FDR = 0.04) and PCDD66 (FDR = 0.0006). GO analyses of the identified DMPs revealed no significant enriched pathways (FDR > 0.05; additional file 1 ST19-ST24).

Differentially methylated regions (DMRs)

Eight significantly differentially methylated regions (DMRs) were identified in sex-stratified EWAS of three different PCB and PCDD congeners (FDR-adjusted p < 0.05; mean beta-fold-changes are presented in Table 2(A, B); additional file 1 ST8 & ST10). In females, one DMR was identified. An increasing level of exposure to PCDD67 was associated with increased DNAm levels (mean beta-fold-change = 0.20, FDR = 6.52×10⁻⁹; Table 2A, additional file 1 ST10) across three CpGs located in a DMR which overlaps the promoter region (± 2,000 base pairs [bp] from the transcription start site [TSS]) of the gene coding for phosphodiesterase 6B (PDE6B; Figure 3(A, B).

Table 2. Differentially methylated regions associated with PCB and PCDD congeners in A) female or male B) subsample.

A) Female subsample.
Congener	Chromosome: start-end	width	Number of CpGs	CpG-ID (chr:pos strand)	Overlapping promoter region	meanbetafc	FDR cut-off
PCDD67	Chr4:645,256–645,495	240	3	cg27357918 (chr4:645,256 -) cg0105263 (chr4:645,440 +) cg15424112 (chr4:645,495 -)	PDE6B	0.20	6.52x10^−09
B) Male subsample.
Congener	Chromosome: start-end	width	Number of CpGs	CpG-ID (chr:pos strand)	Overlapping promoter region	meanbetafc	FDR cut-off
PCB126	Chr8:101,348,639–101,348,921	283	4	cg27190410 (chr8:101,348,639 +) cg08096750 (chr8:101,348,741 -) cg16149401 (chr8:101,348,747 -) cg12182643 (chr8:101,348,921 -)	RNF19A	−0.26	1.12x10^−15
PCDD66	Chr7:2,558,826–2,558,891	66	3	cg08947335 (chr7:2,558,826 +) cg06368153 (chr7:2,558,840 +) cg17130815 (chr7:2,558,891 -)	LFNG	0.03	8.71x10^−10
	Chr5:140,079,898–140,080,575	678	6	cg24175899 (chr5:140,079,898 +) cg02659009 (chr5:140,079,941 +) cg06225581 (chr5:140,080,094 +) cg09170127 (chr5:140,080,245 -) cg22976331 (chr5:140,080,544 -) cg03975922 (chr5:140,080,575 +)	ZMAT2	−0.12	8.71x10^−10
	Chr18:44,237,619–44,237,658	40	2	cg13561807 (chr18:44,237,619 -) cg19829856 (chr18:44,237,658 -)	LOXHD1	0.15	8.71x10^−10
	Chr15:40,268,421–40,268,777	357	5	cg02712587 (chr15:40,268,421 -) cg01081584 (chr15:40,268,610 -) cg20255370 (chr15:40,268,687 +) cg01729180 (chr15:40,268,744 -) cg16127683 (chr15:40,268,777 -)	EIF2AK4	−0.19	8.71x10^−10
	Chr19:49,244,004–49,244,126	123	4	cg05432996 (chr19:49,244,004 -) cg10347828 (chr19:49,244,040 -) cg17773655 (chr19:49,244,121 +) cg09810593 (chr19:49,244,126 +)	RASIP1, IZUMO1	0.14	8.71x10^−10
	Chr6:26,158,157–26,158,396	240	7	cg21316379 (chr6:26,158,157 -) cg08131845 (chr6:26,158,202 +) cg26376783 (chr6:26,158,315 -) cg02887499 (chr6:26,158,339 -) cg23524294 (chr6:26,158,342 +) cg26741595 (chr6:26,158,391 +) cg14108534 (chr6:26,158,396 +)	H2BC5	0.00	8.71x10^−10

Note: CpGs = Number of significantly (FDR < 0.05) differentially methylated CpGs in DMR, chr:pos = chromosome:position, overlapping promotor region in 2000 ± base pairs (bp) from transcription start site (TSS), meanbetafc = mean beta fold change within the region, FDR cut-off = significance FDR corrected p-value threshold, – = no annotation.

Figure 3. DMR associated with PCDD67 in female subsample. A) Differentially methylated region (DMR) associated with PCDD67 in female neonates. In green, the position of each assayed CpG in the DMR is shown with its -log10 transformed FDR above the FDR-corrected significance threshold (blue horizontal line). The DMR´s location related to NCBI RefSeq gene annotation (GRCh37/hg19), in a range of ± 2,000 bp from TSS is shown at the bottom of the plot. Here, it displays three PDE6B isoforms (NM_001350155.2, NM_001350154.2, & NM_001145292.1), and in addition two other PDE6B isoforms (NM_001145291.1 & NM_000283.3), that are located some 25kbp upstream of the chromosomal locus. B) For purpose of illustration only, the mean beta values of the three correlated (all p < 0.01; additional file 1 ST11) CpGs located in the DMR in chromosome 4 are linearly regressed on PCDD67 in female subsample without confounder adjustment and robustness for outliers in a scatterplot.

In males, seven DMRs were identified (Table 2B; additional file 1 ST10; additional file 2 supplementary figure (SF) 1A-SF5B). One DMR is located in the promoter region of H2B Clustered Histone 5 (H2BC5; Figure 4(A, B)). It comprises seven adjacent CpGs, including cg02887499, which was also identified as a DMP (see above). Cg02887499 [chr6:26,158,339] is located within the transcription factor-binding site motifs for YY1 (chr6:26,158,329–26,158,343), and TBP (chr6:26,158,320–26,158,339). Furthermore, this DMR includes the transcription-binding site of NFYA (chr6:26,158,003–26,158,322) and NFYB (chr6:26,158,058–26,158,316). GO analyses of the identified DMRs revealed no significant enriched pathways (FDR > 0.05; additional file 1 ST25-ST27).

Figure 4. DMR associated with PCDD66 in male subsample. A) Differentially methylated region (DMR) in chromosome 6 associated with PCDD66 in male neonates. In green, the position of each assayed CpG in the DMR is shown with its -log10 transformed FDR above the FDR-corrected significance threshold (blue horizontal line). The light red marked CpG position also represents a significant DMP (cg02887499). The DMR´s location related to NCBI RefSeq gene annotation (GRCh37/hg19) shows two H2BC5 isoforms (NM_138720 & NM_021063). B) For purpose of illustration only, the beta values of the seven uncorrelated (all p > 0.05; additional file 1 ST18) CpGs located in the DMR in chromosome 6 are linearly regressed on PCDD66 in male subsample without confounder adjustment and robustness for outliers in a scatterplot.

Weighted gene co-methylation network analysis (WGCNA)

Two weighted gene co-methylation networks were constructed separately for females and males. In female samples, 174 modules were identified from the 50% most variable CpG sites (380,202 CpGs). In males, 113 modules were identified from the 50% most variable CpG sites (380,232 CpGs).

EigenCpG-values, reflecting the first principal component of each module, were further used as indicator of the whole module and linearly regressed on the individual congeners and confounders in same order as in the EWAS-analyses in a sex-stratified manner. The first principal component in each module explained approximately half of the total variance, while the other components of a module explained minimal variance (selected scree plots in additional file 2 SF6-SF11). We did not observe any significant associations between EigenCpG-values and any PCB or PCDD/F predictor in male neonates. In females, 31 modules were identified to be significantly associated with exposure levels of PCDF130 after adjustment for multiple testing (FDR < 0.05; Table 3). For 20 of these modules, there was a negative association between the module’s EigenCpG-value and the PCDF130-congener, indicating that with increasing level of PCDF130 the weighted average DNAm profile of a module decreased.

Table 3. 31 WGCNA-identified modules significantly associated with PCDF130-congener in female subsample adjustment for confounders^a.

module description				regression
color	size	explained variance (%)	Hub-CpG	b [95% CI]	SE	p	FDR
grey60	4308	42.27	cg13984060	−.33 [−.51, −.15]	.09	.0005	.0209
saddlebrown	2694	30.39	cg05913684	.17 [.08, .26]	.05	.0004	.0209
salmon	4599	29.37	cg01733928	.24 [.12, .37]	.06	.0002	.0209
steelblue	2673	31.21	cg02124724	−.3 [−.46, −.14]	.08	.0004	.0209
blue1	28	57.54	cg10490842	−.27 [−.45, −.10]	.09	.003	.0300
brown	14,768	51.66	cg09899564	.09 [.04, .15]	.03	.002	.0300
coral4	54	52.92	cg13274612	−.22 [−.36, −.08]	.07	.003	.0300
cyan	4595	30.08	cg05098339	−.27 [−.43, −.11]	.08	.002	.0300
green4	65	47.48	cg07479988	−.27 [−.44, −.10]	.08	.002	.0300
lightcyan1	1431	51.44	cg11468953	−.07 [−.12, −.03]	.02	.003	.0300
magenta4	172	33.40	cg27420610	−.12 [−.19, −.04]	.04	.002	.0300
mediumpurple2	325	32.92	cg18169610	−.24 [−.38, −.10]	.07	.001	.0300
midnightblue	4504	33.96	cg22358291	−.17 [−.27, −.06]	.05	.002	.0300
orangered4	1963	28.54	cg23140706	.18 [.07, .29]	.06	.002	.0300
skyblue2	381	33.05	cg17356252	−.27 [−.44, −.10]	.09	.003	.0300
thistle3	201	49.14	cg21992594	−.27 [−.43, −.10]	.08	.002	.0300
yellow4	378	40.05	cg08560387	.07 [.03, .12]	.02	.001	.0300
darkgrey	3303	30.30	cg05966726	−.12 [−.19, −.04]	.04	.003	.0300
blue4	78	54.70	cg12111351	−.07 [−.11, −.02]	.02	.004	.0317
lightsteelblue	316	40.25	cg22146755	−.1 [−.17, −.04]	.03	.004	.0317
white	2997	31.42	cg10482044	−.24 [−.39, −.08]	.08	.004	.0317
orangered3	356	35.11	cg23344769	−.15 [−.26, −.05]	.05	.005	.0381
antiquewhite1	54	42.23	cg06415891	−.25 [−.43, −.08]	.09	.01	.0432
maroon	682	41.26	cg14389122	.08 [.02, .15]	.03	.01	.0468
plum2	980	33.40	cg00614413	.26 [.07, .45]	.09	.01	.0497
yellow	11,426	42.35	cg09941484	−.19 [−.32, −.05]	.07	.01	.0497
antiquewhite2	119	34.60	cg10388190	−.21 [−.37, −.06]	.08	.01	.0498
ivory	1407	35.14	cg25218840	.13 [.03, .22]	.05	.01	.0498
lavenderblush2	154	45.35	cg11905324	−.25 [−.44, −.07]	.09	.01	.0498
palevioletred1	74	33.51	cg10285122	.18 [.05, .32]	.07	.01	.0498
red	6035	48.80	cg18901378	.25 [.07, .44]	.09	.01	.0498

Note. Size = N of CpGs included in module, explained variance = variance explained by first principal component of EigenCpG analysis, chr = chromosome, pos = position, b = regression coefficient, CI = confidence interval, SE = standard error, p = p-value, FDR = false discovery rate, ^aconfounders were mother’s age at birth, gestational length, mother´s smoking status during pregnancy, neonate´s birth weight, mother´s highest education level, mother’s BMI before pregnancy, estimated fraction of (CD4T + CD8T + monocytes + natural killer cells + granulocytes + B-cells), and batch effect.

For each module significantly associated with PCDF130 in female neonates, the most connected intramodular CpG, the so-called hub-CpG, was identified (Table 3). None of these 31 hub-CpGs were identical to any DMPs or CpGs in DMRs identified in EWAS with female neonates.

Pathway analysis showed that CpGs from six (brown, light cyan, red, salmon, white, & yellow) modules of the 31 modules identified as significantly associated with PCDF130 exposure-levels in female samples were significantly enriched in gene pathways after adjustment for multiple testing (FDR < 0.05; Figure 5; additional file 1 ST28-ST33; additional file 2 SF12-SF25). Analysis of the relationship between significantly enriched GO-terms across modules identified GO-terms related to, among others, neurotransmission, synapse and neuron parts, synapse organizations, and immune response.

Figure 5. Top enriched pathways in six modules associated with PCDF130 in female subsample. Top enriched pathways in the brown, light cyan, red, salmon, white, & yellow module significantly associated with PCDF130 in female neonates. Biological pathway terms are ranked by -log10 FDR. GO aspect and GO ID are shown in parentheses. Note. GO aspects are CC = cellular component, MF = molecular function, BP = biological process.

Discussion

Prenatal exposure to endocrine disrupting chemicals (EDCs) can interfere with developmental processes, and it has been suggested that long-term effects of in utero EDC exposure might be caused by the influence of EDCs on basic gene regulatory mechanism. Altered DNAm levels might be functionally involved in these processes, or serve as biomarkers of exposure. Here, using an epigenome-wide screen, we investigated whether prenatal low-dose exposure to PCBs and PCCD/Fs was associated with altered DNAm patterns in cord blood in a Western urban-industrial population in a sex-stratified manner. Our EWAS showed that six congeners of PCB and PCDD were significantly associated with DNAm levels at 32 CpG sites in female or male neonates. In addition to these DMPs, we identified eight differentially methylated regions (DMRs), one of those in females. The implicated genes are involved in neurodevelopment processes (DISC1, LFNG, TPD52L1, PAQR4), gene regulation (SMC1A, H2BC5), and immune functioning (ALOX5, HSPA1B), and for several of these genes, there is evidence for effects of EDC exposure on gene expression and differential DNAm in rodent models.

Of note, we found one DMP annotated to DISC1 that was associated with PCB126 and PCDD66 levels in male neonates. DISC1 codes for disrupted in schizophrenia 1 and is involved in the regulation of various processes of embryonic neurogenesis. It is especially known for its role in the pathophysiology of schizophrenia and other mental health problems due to perturbation of DISC1 expression [57]. In relation to our finding, a study on mutant Disc1 (mDisc1) mice exposed to lead observed neurobehavior consistent with schizophrenia regarding locomotive activity and other mental disorders such as anxiety [58]. In human studies, it has been shown that prenatal exposure to lead or BPA is associated with increased likelihood of schizophrenia later in life [59–62].

PCDD66 and PCDD75 exposure in male neonates was also associated with DNAm levels of a CpG site and a DMR which mapped to LFNG, also involved in neurodevelopment. It codes for LFNG O-Fucosylpeptide 3-Beta-N-Acetylglucosaminyltransferase, which is involved in neuronal differentiation by acting in the Notch signalling pathway to define boundaries during embryonic development [63]. Altered expression levels of LFNG have also been reported before to be associated with ultrafine particles (UFP) exposed onto COPD-diseased human bronchial epithelial cells [64]. UFP are able to act as endocrine disruptor on the developing organism [65].

An additional DMP associated with PCDD66 in the male samples was located in the gene TPD52L1, coding for Tumour protein D52-like 1, involved in carcinogenesis and cell proliferation. PCDD66 in mixture with other PCDD/Fs, PCBs and pesticides fed to mother rats during pregnancy and nursing, led to differentially expressed Tpd52l1 in liver tissue of the offspring [66]. Furthermore, Tpd52l1 was reported to be differentially expressed in zebrafish exposed to low doses of atrazine, showing endocrine disrupting effects [e.g., 67], during early development [68].

Finally, one DMP related to neurodevelopment that was associated with PCDD66 in male neonates of our study, was annotated for the gene PAQR4. PAQR4 is coding for Progestin And AdipoQ Receptor Family Member 4 and is associated with Isolated Growth Hormone Deficiency (Type Ib), a neurodevelopmental disease that is primarily characterized by short stature and is also associated with cognitive impairment [69–71]. In male rodent neonates, which were exposed to low-dose BPA and 17β-oestradiol-3-benzoate during neonatal development, Paqr4 was differentially methylated in dorsal prostate tissue [72].

Furthermore, we found two DMPs and one DMR associated with genes involved in gene regulation. Most notably, a strong association was found between PCB126, PCDD66 or PCDD70 exposure in males with DNAm for H2BC5. In addition to one DMP in H2BC5, we identified a DMR harbouring seven CpGs in the promoter region of H2BC5. The gene codes for H2B Clustered Histone 5 and is one of the four core histones responsible for the nucleosome structure of the chromosomal fibre in eukaryotes and thus involved in modulating DNA accessibility. The DMR includes several transcription factor-binding site motifs, and cg02887499 overlaps or directly neighbours those for YY1 (Ying-Yang 1) and TBP (TATA-box binding protein).

Lastly, we found DMPs and DMRs annotated to genes involved in immune functioning. In our study, one DMP was associated with PCDD66 in male neonates, which was annotated for the gene ALOX5. ALOX5 is coding for arachidonate 5-Lipoxygenase and plays an important role in the synthesis of leukotrienes from arachidonic acid. Leukotrienes are involved in various inflammatory and allergic processes, linking ALOX5 indirectly to immune processes. In embryo-larvae of zebrafish exposed to chiral PCB149, Alox5 expression was found to be dysregulated [73]. Additionally, we identified one DMP annotated for gene HSPA1B that was associated with PCDD66 in male neonates. HSPA1B is coding for Heat Shock Protein Family A (Hsp70) Member 1B, that is involved in cellular stress regulation by maintaining the proteostasis [74] and plays a mediating role between transcription factor Foxp3 and ubiquitin ligase Sub1 in regulatory T-cells during inflammation [75]. In male offsprings of rats, fed to mixture of PCBs, PCDD/Fs (including PCDD66) and pesticides during pregnancy and nursing, Hspa1b was differentially expressed compared to male offspring of rats that were fed with corn oil [66].

Overall, DNA methylation levels associated with PCB and PCDD exposure have been previously reported as associated with exposure to other EDCs (e.g., BPA), which might suggest a set of genes with broad sensitivity to chemical exposure.

We additionally examined epigenome-wide DNAm data from a systems-level view by conducting WGCNA to identify clusters of highly connected CpGs – so-called modules. Controlling for biological and technical confounders, six modules were significantly associated with PCDF130 in female samples and significantly enriched for gene pathways of neurodevelopment and immune functioning, which converges on the results discussed above.

The following limitations need mention. First, DNAm was assessed in peripheral cord blood cells although – given the hypotheses that prenatal EDC exposure can programme neurodevelopment resulting in long-term changes in behaviour and risk for mental health problems – the brain is the primary tissue of interest. However, DNAm profiles from cord blood might be understood as proxy or marker of epigenetic patterns in neuronal tissue cells, since epigenetic modifications established during early development might become propagated soma-wide [76]. Furthermore, systematic research on tissue concordance showed that blood-based EWAS are suitable to identify DNAm patterns as biomarker for mental phenotypes manifesting in the brain [77]. Second, in our study, DNAm was measured from cord blood tissue, which compromises multiple cell types having distinct epigenetic profiles, including DNAm patterns [78]. Potential confounding effects in the epigenetic profile due to cell-type heterogeneity in cord blood were, however, accounted for [48,78,79]. Third, it is known that a large part of the variance of DNA methylation levels is accounted for by genetic variation. Whereas we accounted for SNPs overlapping probes on the array, we cannot rule out potential effects of methylation quantitative trait loci (methQTLs) as potential confounders. Fourth, the identified differentially methylated CpGs were not validated. Local deep bisulphite sequencing of the CpGs on MiSeq (Illumina, San Diego, CA, USA) would be an adequate approach to evaluate the results further. Lastly, it has to be noted that our EWAS mainly identified associations between exposure for PCB and PCDD congeners and DNAm in the male sample, whereas significant associations between exposure for PCDF congener and DNAm of co-methylated modules was limited to females. These results might indicate that there is a sex-dependent susceptibility to the different congeners of PCB and PCDD/F [13]. However, it also might be reasonable that the small sample size of each sex introduced an issue of statistical power of unidentified false-negative results.

Little is known about sex-specific effects of EDC exposure on DNAm levels. Potential mechanisms underlying the sex-dependent susceptibility and effects of the different congeners of PCB and PCDD/F might be that these chemicals are known to interact with the sex-hormonal system [28]. Different congeners of PCB and PCDD/F seem to have different binding affinities to sex-hormone nuclear receptors involved in the gene expression machinery [80,81].

A strength of our study design is the temporal proximity as well as the chronological order of measurements of levels of PCBs and PCDD/F from maternal blood during late pregnancy and DNAm in cord blood after birth, which reduces threat of reverse causation. Furthermore, EDCs levels were objectively quantified using mass spectrometry in maternal blood. Future confirmatory studies will be required to follow up on these results, and the present findings should be considered hypothesis generating rather than definitive. For instance, to further understand the regulatory effects of our identified differentially DNAm patterns associated with prenatal exposure to PCBs and PCDD/Fs, cell culture-based experiments with varying exposure levels followed by downstream analyses of gene expression are warranted. Regarding possible long-term effects, it will be of interest whether the identified differentially methylated marks in cord blood are associated with behavioural outcomes in general and in specific with developmental outcomes that have been shown to be altered by PCBs and PCDD/Fs. The longitudinal nature of our study will enable causal mediation analyses to investigate the potential role of DNAm in mediating any exposure-outcome associations [82]. Promising future studies are systematic comparisons how the measured levels of this cohort are in line with levels from other populations than the German population. The measured levels of PCBs and PCDD/Fs of the samples in the present study are representative for the German population as qualified by a German reference panel matching for sampling year and sample’s age [43]. A systematic comparison of the exposure levels between different populations would help to get insight into how representative the results of this study might be for other populations.

To conclude, a better understanding of the mechanisms involved in the exposure-outcome relationship between EDCs and mental disorders continues to be an important research aim. Our findings on differentially methylated DNAm levels of single CpG sites, regions and modules indicate that PCBs and PCDD/Fs induce epigenetic changes that are associated with neurodevelopment, gene expression and immune functioning and thereby might set origin for adverse mental and health outcomes over the life span.

Discloser of interest

The authors report no conflict of interest.

Authors’ contributions

RK, AS, & MW designed the study. Data was acquired by MKS & MW. Labor analyses were conducted by MKS, MW, & DM. KM, RK, & AS conceived of the analyses. KM & VV performed all data analyses. ST & MS reviewed the data analyses. KM, VV, & NNW drafted the manuscript. DM, ST, MS, JW, JW, JGH, & AS substantive reviewed. RK, AS, JW, & JGH supervised the research project. All authors have read and approved the manuscript.

Consent

All parents gave written informed consent about the inclusion of material pertaining to themselves and their child, about their acknowledgement that they cannot be identified via the paper, and about the full anonymization of their data. We confirm that all mandatory laboratory health and safety procedures have been complied with in the course of conducting the research reported in our paper.

Supplementary Material

Supplemental data for this article can be accessed here

Additional information

Funding

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under Grant KU 2479/5-1 and under Grant SCHO 1738/19-1 to Robert Kumsta and Axel Schölmerich; Deutsche Forschungsgemeinschaft [KU 2479/5-1]; Deutsche Forschungsgemeinschaft [SCHO 1738/19-1].