Polychlorinated biphenyl exposure and DNA methylation in the Anniston Community Health Survey.

Anniston, Alabama was home to a major polychlorinated biphenyl (PCB) production facility from 1929 until 1971. The Anniston Community Health Survey I and II (ACHS-I 2005-2007, ACHS-II 2013-2014) were conducted to explore the effects of PCB exposures. In this report we examined associations between PCB exposure and DNA methylation in whole blood using EPIC arrays (ACHS-I, n = 518; ACHS-II, n = 299). For both cohorts, 35 PCBs were measured in serum. We modelled methylation versus PCB wet-weight concentrations for: the sum of 35 PCBs, mono-ortho substituted PCBs, di-ortho substituted PCBs, tri/tetra-ortho substituted PCBs, oestrogenic PCBs, and antiestrogenic PCBs. Using robust multivariable linear regression, we adjusted for age, race, sex, smoking, total lipids, and six blood cell-type percentages. We carried out a two-stage analysis; discovery in ACHS-I followed by replication in ACHS-II. In ACHS-I, we identified 28 associations (17 unique CpGs) at p ≤ 6.70E-08 and 369 associations (286 unique CpGs) at FDR p ≤ 5.00E-02. A large proportion of the genes have been observed to interact with PCBs or dioxins in model studies. Among the 28 genome-wide significant CpG/PCB associations, 14 displayed replicated directional effects in ACHS-II; however, only one in ACHS-II was statistically significant at p ≤ 1.70E-04. While we identified many novel CpGs significantly associated with PCB exposures in ACHS-I, the differential methylation was modest and the effect was attenuated seven years later in ACHS-II, suggesting a lack of persistence of the associations between PCB exposures and altered DNA methylation in blood cells.

Introduction

From 1929 until 1971, Anniston, Alabama was home to the major polychlorinated biphenyl (PCB) manufacturing facility in the United States. Releases of PCBs from the Anniston production facility resulted in substantial environmental contamination. Since PCBs are relatively resistant to biological and chemical decomposition and are highly lipophilic, they easily bioaccumulate in the food chain. In 2005, the Anniston Community Health Survey was established to assess the health effects of general population PCB exposure (ACHS-I) with a follow-up study conducted in 2014 (ACHS-II) [1,2]. Previous ACHS-I studies reported associations between PCB exposure and race [1], hypertension [3], blood pressure [4], leukocyte telomere length [5], liver disease [6], metabolic syndrome [7], serum lipid levels [8], and diabetes [9]. Several previous population studies have found associations between PCB or dioxin exposures and altered blood DNA methylation levels in Dutch men [10], Greenlandic Inuits [11], Koreans [12], Faroe Islanders [13] and Japanese [14].

Some PCBs are ligands for the aryl hydrocarbon receptor (AHR), and it is hypothesized that most, but not all, of the PCB-associated adverse health outcomes are related to persistent activation of AHR [15]. In some tissues, the aryl hydrocarbon receptor repressor (AHRR) is upregulated following AHR activation, and AHRR protein acts as an AHR negative regulator [16]. Numerous recent epigenome-wide association studies (EWAS) of tobacco smoke exposure have demonstrated that in smokers the AHRR CpG cg05575921 is strongly hypomethylated. Also, among smokers, AHRR gene expression is upregulated, particularly in CD14+ monocytes [17-19]. We and others have hypothesized that this altered state could be caused by the presence of polyaromatic hydrocarbons (PAH) in tobacco smoke that activate AHR and lead to AHRR hypomethylation and persistent upregulation of AHRR mRNA. PCB exposures, as AHR ligands, could result in similar effects. So, in addition to examining epigenome-wide associations with PCB exposures, we also assessed whether PCBs were associated with the methylation of the AHRR smoking biomarker.

Here we carry out an EWAS examining the associations between DNA methylation levels using the Illumina EPIC array and the sum of 35 measured PCBs (Σ35PCBs), as well as various exposure groupings: mono-ortho substituted PCBs, di-ortho substituted PCBs, tri/tetra-ortho substituted PCBs, and PCBs with putative oestrogenic or antiestrogenic activity [9,20-22].

Methods

Study population

ACHS-I was conducted between 2005 and 2007, and details of the study design have been reported [1,23]. In ACHS-I, we recruited 765 participants using random-stratified sampling, and residents in the area closest to the PCB manufacturing facility (west Anniston) were over-sampled (two-thirds of eligible participants). The residents were interviewed and had blood samples taken for serum level measurement of 35 PCB congeners, serum lipids, and DNA isolation. We measured DNA methylation in a subset (n = 518) of the ACHS-I cohort members using DNA from stored frozen blood clot samples [5].

ACHS-II was conducted between in 2014 as a follow-up study [2]. For ACHS-II, surviving participants from ACHS-I were contacted and 359 eligible individuals were enrolled. The same 35 PCBs measured in ACHS-I were determined again in ACHS-II. DNA methylation was measured in a subset (n = 299) of the ACHS-II cohort. Both studies were approved by the Institutional Review Boards at the Center of Disease Control and the University of Alabama at Birmingham. See Table 1 for study demographics.

Table 1.

Characteristics of the Anniston Community Health Survey, phases I and II.

Characteristic	African-Americans(n = 208)	ACHS-I Whites (n = 310)Mean ± SE	Total (n = 518)	African-Americans (n = 155)	ACHS-II Whites (n = 144)Mean ± SE	Total (n = 299)
Age	53.3 ± 1.1	55.6 ± 1.0	54.6 ± 0.7	61.2 ± 0.9	64.5 ± 1.2	62.8 ± 0.8
Body Mass Index (BMI)	31.8 ± 0.6	30.5 ± 0.4	31.1 ± 0.3	32.3 ± 0.6	30.9 ± 0.7	31.6 ± 0.5
Years of residence	26.0 ± 1.3	32.0 ± 1.1	29.6 ± 0.9	-	-	-
Residential distance from plant (km)	2.5 ± 0.1	4.3 ± 0.1	3.6 ± 0.1	-	-	-
		n (%)			n (%)
Female	148 (71.2)	221 (71.3)	369 (71.2)	120 (77.4)	104 (72.2)	224 (74.9)
Age Group (years)
<40	39 (18.8)	61 (19.7)	100 (19.3)	6 (3.9)	10 (6.9)	16 (5.4)
40-59	99 (47.6)	111 (35.8)	210 (40.5)	64 (41.3)	41 (28.5)	105 (35.1)
≥ 60	70 (33.7)	138 (44.5)	208 (40.2)	85 (54.8)	93 (64.6)	178 (59.5)
BMI Class (kg/m2)
<25	40 (19.2)	69 (22.3)	109 (21.0)	28 (18.1)	35 (24.3)	63 (21.1)
25-29	43 (20.7)	93 (30.0)	136 (26.3)	41 (26.5)	43 (29.9)	84 (28.1)
≥ 30	123 (59.1)	148 (47.7)	271 (52.3)	86 (55.5)	66 (45.8)	152 (50.8)
Current smoker	75 (36.1)	101 (32.6)	176 (34.0)	36 (23.2)	30 (20.8)	66 (22.1)
Residence in west Anniston	180 (86.5)	249 (80.3)	429 (82.8)	145 (93.5)	115 (79.9)	260 (87.0)
Occupational PCB exposure	53 (25.5)	74 (23.9)	127 (24.5)	33 (21.3)	40 (27.8)	73 (24.4)

Serum PCBs and lipid measurements

The CDC’s National Center for Environmental Health laboratory analysed PCBs in serum [1]. In both ACHS-I and ACHS-II, 35 PCBs congeners were measured: 28, 44, 49, 52, 66, 74, 87, 99, 101, 105, 110, 118, 128, 138–158, 146, 149, 151, 153, 156, 157, 167, 170, 172, 177, 178, 180, 183, 187, 189, 194, 195, 196–203, 199, 206, and 209. Wet-weight substituted measurements were grouped into seven categories: Σ35PCBs (the sum of all 35 PCBs), mono-ortho substituted PCBs (28, 66, 74, 105, 118, 156, 157, 167, and 189); di-ortho substituted PCBs (99, 138–158, 146, 153, 170, 172, 180, and 194); and tri/tetra-ortho substituted PCBs (177, 178, 183, 187, 195, 196–203, 199, 206, and 209), oestrogenic PCBs, E1 (PCBs 66, 74 and 99) [9,20] and E2 (PCBs 99 and 153) [24], and PCBs with antiestrogenic activity, AE (PCBs 66, 74, 105, 118, 156, and 167) [21] (Table 2). PCB values below the limit of detection were imputed using the limit of detection divided by square root of 2. All PCBs were used in calculating the Σ35PCBs. For the ortho-substituted and oestrogenic/antiestrogenic groups, we excluded any PCB for which ≥ 60% individuals had measures below the limit of detection (PCBs 44, 49, 52, 87, 101, 110, 128, 149, and 151).

Table 2.

PCBs measured in ACHS-I and ACHS-II.

Σ35PCBs		Mono-ortho substituted PCBsa	Di-ortho substituted PCBsa	Tri/tetra-ortho substituted PCBsa	Oestrogenic PCBs Group 1b	Oestrogenic PCBs Group 2c	Antiestrogenic PCBsd
PCB28	PCB156	PCB28	PCB99	PCB177	PCB66	PCB99	PCB66
PCB44	PCB157	PCB66	PCB138-158	PCB178	PCB74	PCB153	PCB74
PCB49	PCB167	PCB74	PCB146	PCB183	PCB99		PCB105
PCB52	PCB170	PCB105	PCB153	PCB187			PCB118
PCB66	PCB172	PCB118	PCB170	PCB195			PCB156
PCB74	PCB177	PCB156	PCB172	PCB196-203			PCB167
PCB87	PCB178	PCB157	PCB180	PCB199
PCB99	PCB180	PCB167	PCB194	PCB206
PCB101	PCB183	PCB189		PCB209
PCB105	PCB187
PCB110	PCB189
PCB118	PCB194
PCB128	PCB195
PCB138-158	PCB196-203
PCB146	PCB199
PCB149	PCB206
PCB151	PCB209
PCB153

PCB values (wet-weight substituted, ng/g) below the limit of detection were imputed by dividing the limit of detection for the assay by the square root of 2.

aThe following PCBs were excluded because ≥ 60% of participants had levels below the limit of detection: PCB44, PCB49, PCB52, PCB87, PCB101 PCB110, PCB128, PCB149,

and PCB151.

bOestrogenic PCBs Group 1 also includes PCBs 44, 49, 110, and 128 (Silverstone et al., 2012; DeCastro et al. 2006).

cOestrogenic PCBs Group 2 also includes PCBs 52, 101, and 110 (Warner et al., 2012).

dAntiestrogenic PCBs includes PCBs 66, 74, 105, 118, 156, 167 (Warner et al., 2012 combined Wolff and Cooke groupings).

The Clinical Chemistry Laboratory at the Jacksonville Medical Center (Jacksonville, AL) measured serum lipid levels for total cholesterol, HDL cholesterol, LDL cholesterol, and triglycerides [1]. Total lipids were calculated using the method as described by Bernet et al. [25]

DNA methylation measurement

For ACHS-I, DNA was extracted from stored blood clot samples using the PureGene (QIAGEN, 158,389) protocol and Clotspin Baskets (QIAGEN, 158,932). For ACHS-II, whole blood samples (300 µL per extraction) were aliquoted into deep well 96-well plates and genomic DNA was isolated using the Agencourt Genfind v2 Solid Phase Reversible Immobilization (SPRI) paramagnetic bead-based technology (Beckman Coulter, A41497). Prior to bisulphite-conversion, we measured DNA concentrations using a QUBIT dsDNA BR assay kit (Invitrogen, Q32850). The National Cancer Institute’s Cancer Genomics Research Laboratory performed DNA bisulphite conversion using the EZ-96 DNA Methylation MagPrep kit (Zymo Research, D5040) and then ran the bisulphite-converted DNA on Illumina EPIC methylation arrays (Illumina, WG-317-1001).

Statistical analysis

Bioconductor’s ChAMP package v2.12.0 [26-28] was used to normalize and batch correct methylation data. We excluded any samples which failed array QC standards. We excluded from analysis all CpG probes on the X and Y chromosomes, probes containing a SNP with a minor allele frequency ≥ 1% in the CpG site, and probes failing QC standards. We also removed an additional 43,254 probes reported to hybridize to one or more non-target sites in the genome [29]. There were 706,625 CpG probes remaining after exclusions.

Figure 1 outlines our analysis strategy. In brief, associations between PCB exposures and methylation in ACHS-I were analysed using robust multivariable linear regression (M-estimation) with the rlm function in Modern Applied Statistics with S (MASS v.7.3–51.1) [30]. We log10-transformed PCB measurements, consistent with previous ACHS-I analyses. All models were adjusted for age, race, sex, smoking status (yes/no), log10-transformed total lipids, bisulphite-conversion batch, and white blood cell percentages. Other covariates e.g., body mass index and alcohol use, were tested in our models, but did not contribute to model fit. The six white blood cell types (CD4 and CD8 T-cells, B-cells, monocytes, natural killer cells, and granulocytes) were estimated using the Houseman et al. method [31] based on the Reinius reference panel [32]. For CpGs significant at Bonferroni threshold p ≤ 6.70E-08, differential methylation (ΔM) was calculated between the highest and lowest PCB exposure quartiles for CpG sites in ACHS-I. Benjamini-Hochberg false discovery rate (FDR) significant CpGs at p ≤ 5.00E-02 (n = 283) in ACHS-I were then analysed in ACHS-II.

Figure 1.

Overview of data analyses strategy. Robust linear regression models were adjusted for age, race, sex, smoking status, total serum lipids, and percentages for six different white blood cell types. Seven different PCBs groupings were used in regression: 1) Σ35PCBs, 2) mono-ortho substituted PCBs, 3) di-ortho substituted PCBs, 4) tri/tetra-ortho substituted PCBs, 5) oestrogenic PCBs group 1, 6) oestrogenic PCBs group 2, 7) antiestrogenic PCBs. Stratified and DMR analyses were not conducted in ACHS-II.

Potential effect measure modification by race and sex for ACHS-I CpGs with an absolute ΔM ≥ 1.00% was evaluated in ACHS-I. We identified effect modification by including an interaction term in our regression models. For CpGs with nominally significant interaction terms (p ≤ 5.00E-02), we generated stratified ΔM results. We conducted interaction analyses using SAS v9.4 PROC ROBUSTREG (SAS Institute Inc.). The ACHS-II had a smaller sample size, which resulted in our having lower statistical power to detect true effect measure modification. As a result, we did not conduct stratification analyses for ACHS-II CpGs.

Differential methylated regions (DMRs) in ACHS-I were identified using Bioconductor’s DMRcate v.1.18.0 [33], comparing the highest versus lowest tertiles for each exposure. For each DMR we calculated FDR p-values. The model was adjusted for age, race, sex, smoking status (yes/no), log10-transformed total lipids, bisulphite-conversion batch, and white blood cell percentages. Again, because of the smaller sample size in ACHS-II, we did not perform DMR analysis in ACHS-II.

CpG/DMR functional annotations and associations in the comparative toxicogenomics database

To find functional annotations of FDR CpGs and DMRs, we used The Genomic Regions Enrichment of Annotations Tool (GREAT) method [34]. We also retrieved curated chemical–gene association data from the 2019 Comparative Toxicogenomics Database (CTD) [35], MDI Biological Laboratory, Salisbury Cove, Maine, and NC State University, Raleigh, North Carolina. World Wide Web (URL: http://ctdbase.org). [May 2019].

Results for demographics comparisons used a two-sided Welch’s t-test. For testing involving PCB distributions by demographics and PCB exposure quartiles, we used Kruskal-Wallis tests and Wilcoxon rank-sum pair-wise tests, respectively. Statistical analyses were conducted in R [36], SAS v9.4 and JMP 13.0.0 (SAS Institute Inc.).

Results

Demographics

Demographics for individuals available for DNA methylation analysis in ACHS-I and ACHS-II are found in Table 1. Participants were mostly female (ACHS-I female, 71.2% and ACHS-II female, 74.9%) and were heavier (ACHS-I BMI, 31.1 ± 0.3; ACHS-II BMI, 31.6 ± 0.5) than the U.S. general population as measured in the National Health and Nutrition Examination Survey 2015–2016 [37]. African-American participants were significantly more likely to live closer to the PCB production facility (ACHS-I African-Americans’ residential distance from plant, 2.5 km ± 0.1 and whites’, 4.3 km ± 0.1; t-test p = 2.51E-27).

Serum PCB levels by age and race

Serum PCB levels in the ACHS-I were previously reported to be 1.5- to 3.0-fold higher than the general population [1]. Figure 2 shows the distribution of the Σ35PCBs by age group and race in both ACHS-I and ACHS-II. In ACHS-I, mean PCB levels were significantly different among the three age groups: 19–39 v. 40–59, p = 4.87E-32; 19–39 v. ≥ 60, p = 1.27E-40; 40–59 v. ≥ 60, p = 4.50E-15 (Figure 2a) and between African-Americans and whites, p = 5.14E-15 (Figure 2b). These distributions were also similar in ACHS-II: age groups; 19–39 v. 40–59, p = 4.17E-10; 19–39 v. ≥ 60, p = 3.90E-11; 40–59 v. ≥ 60, p = 6.60E-08 and race, p = 1.87E-11 (Figure 2c,d, respectively). We also compared the correlation of Σ35PCBs between ACHS-I and ACHS-II for the 245 subjects with exposure data in both cohorts (Figure S1, r2 = 0.86, p = 1.30E-107). The average measured Σ35PCBs in ACHS-II was reduced approximately 10%, but differences by age and race varied considerably as they did in ACHS-I.

Figure 2.

Sum of 35 PCB levels by age group and race for ACHS-I (n = 518) and ACHS-II (n = 299) participants. PCB values are graphed on a log10 scale. Box/violin plot labels are non-transformed mean PCB level (ng/g wet-weight). (a) ACHS-I age group distribution. (b) ACHS-I race distribution. (c) ACHS-II age group distribution. (d) ACHS-I race distribution. Wilcoxon rank-sum test p-values are presented.

Cell type analysis in ACHS-I

As described in the Method section, DNA methylation profiles were used to estimate distributions of six blood cell types, and the percentages were used as covariates in the linear regression models. Because PCB exposures have been associated with cell type shifts and other immune system effects [38], we hypothesized that exposure might alter methylation-based estimated cell type composition. We tested if any PCB exposures were associated with shifts in estimated immune cell composition in ACHS-I and observed no significant results.

Smoking and methylation in ACHS-I

We examined if tobacco smoke exposure in the ACHS-I population was associated with altered DNA methylation at the AHRR CpG cg05575921, a site which has been confirmed to be a biomarker of tobacco smoke exposure in numerous studies [17-19,39-42]. The multivariable robust linear regression model that tested smoking versus methylation at AHRR cg05575921 was highly significant (smoking regression coefficient p = 7.27E-269, r2 = 0.58; adjusted for age, sex, race, lipids, cell types and batch). The differential methylation of AHRR cg05575921 between non-smokers and smokers was −18.3% (Wilcoxon rank-sum pair-wise test, p = 1.62E-64). There was no association between the smoking biomarker AHRR cg05575921 and PCB levels in non-smokers (p = 9.44E-01).

PCBs and DNA methylation in ACHS-I and ACHS-II

In ACHS-I, we found 28 significant genome-wide associations, representing 17 unique genes, (Bonferroni threshold p ≤ 6.70E-08) between CpG methylation and PCB exposure, across the seven PCB exposure groups: Σ35PCBs, mono-ortho substituted PCBs, di-ortho substituted PCBs, tri/tetra-ortho substituted PCBs, two groups of PCBs with oestrogenic activity, and PCBs with antiestrogenic activity (Table 3) and 369 associations (286 unique CpGs) at the FDR p ≤ 5.00E-02 threshold (Table S1). For the 28 Bonferroni significant associations in ACHS-I, absolute ΔM ranged from 0.06% to 5.02%, with 15 (54%) of the associations being hypomethylation and ten associations having an absolute ΔM ≥ 1.00%. Ten of the 28 (36%) associations were most significant in the tri/tetra-ortho-substituted PCB exposure group (PCBs 177, 178, 183, 187, 195, 196–203, 199, 206, and 209). The Venn diagram in Figure 3 illustrates the distribution of FDR significant CpGs (n = 280) in ACHS-I across the four ortho-substituted exposure groups. The majority of CpGs (n = 229, 82%) were significant exclusively in the tri/tetra-ortho substituted PCB exposure group. Figure S2 shows the Venn diagram of the 15 FDR significant CpGs for the oestrogenic/antiestrogenic groups (E1, E2, AE). There was only one CpG, cg04991747, in polypeptide N-acetylgalactosaminyltransferase 2 (GALNT2), that was FDR significant among all seven PCB exposure groups.

Table 3.

Top ACHS-I PCB-associated CpGs, selected by Bonferroni (BF) threshold p ≤ 6.70E-08.

ProbeID	GeneID	Coordinate	ACHS-I PCB BF Significant Exposuresa	ACHS-I PCB Exposure	ACHS-I % Differential Methylation (Mean ± SD)b	ACHS-Ip-valuec	ACHS-IPCB Regression Coefficent	ACHS-II % Differential Methylation (Mean ± SD)b	ACHS-IIp-valuec	ACHS-IIPCB Regression Coefficent	Replication in ACHS-II	PCB/Dioxin Associationd
cg14251777	GIMAP8	chr7:150,147,655	|S|-|D|-|-|-|-|	Σ35PCBs	−0.50 ± 0.16	4.13E-08	−0.0061	−0.03 ± 0.20	7.47E-01	−0.0004	Yes	PCB28, PCB52, PCB101, PCB138, PCB153, PCB180
cg12803754	OTULIN	chr5:14,676,460	|S|-|-|T|-|-|-|	Σ35PCBs	−1.16 ± 0.38	3.42E-08	−0.0135	−0.17 ± 0.48	4.21E-01	−0.0029	Yes
cg00475490	PRSS23	chr11:86,517,110	|S|-|D|T|-|-|-|	Σ35PCBs	−4.89 ± 0.60	2.47E-08	−0.0224	−3.12 ± 0.73	2.79E-03	−0.0181	Yes	TCDD, PCB28, PCB52, PCB101, PCB138, PCB153, PCB180
cg00941989	SH3TC1	chr4:8,232,817	|S|-|-|T|-|-|-|	Σ35PCBs	0.14 ± 0.32	4.89E-08	0.0103	−0.72 ± 0.41	8.38E-01	−0.0007	No	TCDD
cg25153882	CTTNBP2	chr7:117,499,375	|-|M|-|-|-|-|AE|	Mono-ortho PCBs	−0.29 ± 0.28	1.41E-08	0.0086	−0.58 ± 0.30	8.97E-02	−0.0039	No
cg04991747	GALNT2	chr1:230,207,211	|-|-|-|D|E1|E2|AE|	Di-ortho PCBs	1.12 ± 0.25	3.40E-08	0.0110	0.14 ± 0.32	9.68E-01	−0.0001	No	TCDD
cg14251777	GIMAP8	chr7:150,147,655	|S|-|D|-|-|-|-|	Di-ortho PCBs	−0.51 ± 0.16	1.65E-08	−0.0059	0.02 ± 0.20	6.69E-01	−0.0006	Yes	PCB28, PCB52, PCB101, PCB138, PCB153, PCB180
cg00475490	PRSS23	chr11:86,517,110	|S|-|D|T|-|-|-|	Di-ortho PCBs	−5.02 ± 0.60	6.10E-09	−0.0224	−3.09 ± 0.70	7.00E-03	−0.0160	Yes	TCDD, PCB28, PCB52, PCB101, PCB138, PCB153, PCB180
cg18861197	PTK2B	chr8:27,209,336	|-|-|D|T|-|-|-|	Di-ortho PCBs	−0.83 ± 0.30	6.41E-08	−0.0108	0.07 ± 0.42	4.85E-01	−0.0019	Yes	TCDD, PCB28, PCB52, PCB101, PCB126, PCB138, PCB153, PCB180
cg21566642	AC068134.2	chr2:233,284,661	|-|-|-|T|-|-|-|	Tri/tetra-ortho PCBs	−2.00 ± 1.06	6.13E-08	−0.0382	−1.64 ± 1.43	8.98E-01	−0.0015	Yes
cg25506215	CA9	chr9:35,670,577	|-|-|-|T|-|-|-|	Tri/tetra-ortho PCBs	0.55 ± 0.27	2.16E-08	0.0112	−0.80 ± 0.39	5.30E-01	−0.0020	No	TCDD
cg06005913	CEP97	chr3:101,451,428	|-|-|-|T|-|-|-|	Tri/tetra-ortho PCBs	−0.06 ± 0.24	1.89E-08	0.0059	−0.26 ± 0.23	6.73E-01	0.0007	Yes	TCDD
cg18974766	CPNE5	chr6:36,775,640	|-|-|-|T|-|-|-|	Tri/tetra-ortho PCBs	0.12 ± 0.12	1.35E-08	0.0053	−0.17 ± 0.15	7.24E-02	−0.0024	No
cg26266985	LOC101928185	chr2:133,617,227	|-|-|-|T|-|-|-|	Tri/tetra-ortho PCBs	0.38 ± 0.25	4.20E-08	0.0111	−0.54 ± 0.34	9.24E-01	−0.0003	No	TCDD
cg12803754	OTULIN	chr5:14,676,460	|S|-|-|T|-|-|-|	Tri/tetra-ortho PCBs	−0.97 ± 0.38	6.45E-09	−0.0128	0.01 ± 0.49	2.60E-01	−0.0041	Yes
cg23235442	PFDN4	chr20:52,824,228	|-|-|-|T|-|-|-|	Tri/tetra-ortho PCBs	−0.28 ± 0.08	5.38E-08	−0.0034	0.02 ± 0.09	1.60E-01	0.0013	No	TCDD, Chlorodiphenyl (54% Chlorine)
cg00475490	PRSS23	chr11:86,517,110	|S|-|D|T|-|-|-|	Tri/tetra-ortho PCBs	−4.93 ± 0.60	8.23E-09	−0.0217	−3.00 ± 0.72	1.33E-04	−0.0207	Yes	TCDD, PCB28, PCB52, PCB101, PCB138, PCB153, PCB180
cg18861197	PTK2B	chr8:27,209,336	|-|-|D|T|-|-|-|	Tri/tetra-ortho PCBs	−0.89 ± 0.30	1.93E-08	−0.0106	0.18 ± 0.44	2.67E-01	−0.0029	Yes	TCDD, PCB28, PCB52, PCB101, PCB126, PCB138, PCB153, PCB180
cg00941989	SH3TC1	chr4:8,232,817	|S|-|-|T|-|-|-|	Tri/tetra-ortho PCBs	0.28 ± 0.32	3.23E-08	0.0099	−0.47 ± 0.39	5.74E-01	0.0018	Yes	TCDD
cg04991747	GALNT2	chr1:230,207,211	|-|-|-|D|E1|E2|AE|	Oestrogenic PCBs Group 1	1.13 ± 0.26	1.84E-08	0.0107	0.17 ± 0.31	9.28E-01	0.0002	Yes	TCDD
cg24512644	ST13	chr22:41,248,902	|-|-|-|-|E1|-|AE|	Oestrogenic PCBs Group 1	0.43 ± 0.37	5.30E-08	0.0138	0.16 ± 0.48	1.86E-02	−0.0089	No	TCDD
cg08435853	CCDC166	chr8:144,790,729	|-|-|-|-|-|E2|-|	Oestrogenic PCBs Group 2	2.33 ± 0.91	3.98E-08	0.0402	2.10 ± 1.42	4.22E-01	0.0097	Yes	TCDD
cg04991747	GALNT2	chr1:230,207,211	|-|-|-|D|E1|E2|AE|	Oestrogenic PCBs Group 2	1.27 ± 0.26	7.05E-09	0.0092	0.35 ± 0.31	8.80E-01	−0.0004	No	TCDD
cg11821245	LDHC	chr11:18,433,683	|-|-|-|-|-|E2|-|	Oestrogenic PCBs Group 2	2.81 ± 2.61	4.35E-08	0.0423	2.36 ± 3.26	9.25E-01	−0.0015	No	TCDD
cg25153882	CTTNBP2	chr7:117,499,375	|-|M|-|-|-|-|AE|	Antiestrogenic PCBs	−0.28 ± 0.28	7.39E-09	0.0084	−0.55 ± 0.30	2.71E-01	−0.0026	No
cg04991747	GALNT2	chr1:230,207,211	|-|-|-|D|E1|E2|AE|	Antiestrogenic PCBs	0.96 ± 0.27	6.17E-08	0.0105	0.30 ± 0.32	9.67E-01	−0.0001	No	TCDD
cg24512644	ST13	chr22:41,248,902	|-|-|-|-|E1|-|AE|	Antiestrogenic PCBs	0.30 ± 0.35	4.63E-08	0.0130	−0.06 ± 0.49	1.13E-02	−0.0092	No	TCDD
cg09639964	TIMM22	chr17:909,300	|-|-|-|-|-|-|AE|	Antiestrogenic PCBs	−0.50 ± 0.79	6.41E-08	0.0275	−0.14 ± 0.99	5.25E-01	−0.0037	No	TCDD, 2378-TCDF, PCB28, PCB52, PCB101, PCB138, PCB153, PCB180

CpGs are listed in order of PCB exposure category and GeneID. PCBs with ≥ 60% participants below the limit of detection were excluded from the ortho-substituted and oestrogenic PCB groups.

aBF significance = p ≤ 6.70E-08. PCB exposure groups: S = Sum of all PCBs (n = 35), M = Mono-ortho substituted PCBs (n = 9), D = Di-ortho substituted PCBs (n = 8),

T = Tri/tetra-ortho substituted PCBs (n = 9), E1 = Oestrogenic Group 1 (PCBs 66, 74, and 99); E2 = Oestrogenic Group 2 (PCBs 99 and 153), AE = Anti-oestrogenic Group (PCBs 66, 74, 105, 118, 156, 167).

bDifferential methylation was calculated by subtracting methylation in lowest quartile from the highest.

cRobust linear regression models were adjusted for age, race, sex, smoking status, total serum lipids, bisulphite-conversion batch, and estimated percentages of CD4+ and CD8 + T-cells, CD19 + B-cells, monocytes, granulocytes, and natural killer cells.

dGenes associated with PCB or dioxin exposure based on the 2019 Comparative Toxicogenomics Database.

Figure 3.

Venn diagram of 280 significantly differentially methylated CpGs (BH FDR p ≤ 5.00E-02) in ACHS-I for four PCB groups: sum of 35 PCBs, mono-ortho substituted PCBs, di-ortho substituted PCBs, and tri/tetra-ortho substituted PCBs. Robust linear regression models were adjusted for age, race, sex, smoking status, total serum lipids, bisulphite-conversion batch, and estimated percentages of CD4+ and CD8 + T-cells, CD19 + B-cells, monocytes, granulocytes, and natural killer cells.

The CpG in ACHS-I with the largest absolute ΔM was in the serine protease 23 gene (PRSS23). The PRSS23 CpG cg00475490 was differentially methylated in the Σ35PCBs (−4.89 ± 0.60%, p = 2.47E-08), di-ortho substituted PCBs (−5.02 ± 0.60%, p = 6.10E-09), and tri/tetra-ortho substituted PCBs (−4.93 ± 0.60%, p = 8.23E-09). Figure 4 shows the region of PRSS23 where cg00475490 is located, which is an actively transcribed enhancer region in natural killer cells, as defined by ENCODE/Roadmap projects [43,44]. PRSS23 cg00475490 is 79-bp downstream of a MAF bZIP transcription factor (MAF)/nuclear factor (erythroid-derived 2)-like 2 (NFE2L2, NRF2) transcription binding site.

Figure 4.

PRSS23 genome browser view.

a. Tracks defined from top:Roadmap Epigenome Project Whole Genome Bisulphite Sequencing shows methylation level of each CpG in region.Refseq genes display three mRNA isoforms.CpG island is marked in green.Roadmap Chromatin HMM model for 6 cell types indicates only Natural Killer cells express PRSS23 in blood.ChromHMM code: Red = open chromatin; yellow = active enhancer; Green = active transcription; purple = poised promoter; grey = repressed, closed chromatin.NK Cell H3K27ac ChIP-seqENCODE Transcription Factor ChIP-seqb. Cg004755490 indicated in enhancer region defined by H3K27ac and H3K4me1 histone ChIP-seq in NK cells as determined by Roadmap project. ENCODE transcription factor ChIP-seq for small MAF proteins (binding partners of NRF2).c. Small MAF/NRF2 binding motif located within ChIP-seq peak.

To test the reproducibility of the ACHS-I results, and to see if the observed effects persisted in 2014, we did a replication analysis on the 369 FDR significant (p ≤ 5.00E-02) ACHS-I CpG associations in the ACHS-II. Thirty-three were nominally significant (p ≤ 5.00E-02) in ACHS-II. Of the 28 associations at p ≤ 6.70E-08 in ACHS-I, 14 (50%) were replicated in ACHS-II (based on the directionality of the PCB regression coefficient); however, only one association was statistically significant at Bonferroni p ≤ 1.70E-04 (the Bonferroni threshold for ACHS-II): PRSS23 cg00475490 versus tri/tetra-ortho substituted PCBs, p = 1.33E-04 (Table 3). Of the associations with directional replication both the regression coefficient and the significance were reduced. Absolute ΔM in ACHS-II ranged from 0.01% to 3.12%, with four CpG sites with absolute ΔM ≥ 1.00%.

Quartile analysis in ACHS-I

For the four ACHS-I CpGs significant at Bonferroni p ≤ 6.70E-08 and with an absolute ΔM ≥ 1.00%, we explored the methylation trend across quartiles of PCB exposure (Figure 5). These CpGs were cg21566642 in Homo sapiens BAC clone RP11-762N20 (AC068134.2), cg04991747 in GALNT2, cg12803754 in OTU deubiquitinase with linear linkage specificity (OTULIN), and cg00475490 in PRSS23. For AC068134.2 cg21566642 in the tri/tetra-ortho substituted PCBs, the only significant difference in methylation was between Q1 and Q3, Kruskal-Wallis p = 1.00E-01. For GALNT2 cg04991747, the methylation in the lowest quartile of oestrogenic PCB group 2 was significantly different than the three other quartiles, Kruskal-Wallis p = 2.04E-06. For OTULIN cg00475490, methylation in Q1 of tri/tetra-ortho substituted PCBs differed from Q3 and Q4, but not Q2, Kruskal-Wallis p = 1.98E-03. For PRSS23 cg00475490, Q1 methylation in the di-ortho substituted PCBs was significantly higher than the other three quartiles, Kruskal-Wallis p = 1.52E-17.

Figure 5.

PCB quartile distributions for CpGs with ≥ 1.0% absolute differential methylation. Kruskal-Wallis test and Wilcoxon rank-sum test p-values are presented.

Race and sex stratified analysis in ACHS-I

For ACHS-I Bonferroni significant CpGs with an absolute ΔM ≥ 1.00%, we tested for race and sex specific effects by running multivariable robust linear regression models with an interaction term (either sex or race by exposure) for each of the seven PCB exposure groups and the CpG. For the exposure*CpG interaction terms with p ≤ 1.00E-02, we generated stratified ΔM results. We observed potential effect measure modification by sex for AC068134.2 cg21566642 and tri/tetra-ortho substituted PCBs, with a much larger effect in males (ΔM males −6.11 ± 1.97; ΔM females −0.76 ± 1.22, p for interaction = 8.95E-03). GALNT2 cg04991747 and the di-ortho substituted PCBs displayed a modest difference between males and females (ΔM males 0.87 ± 0.50; ΔM females 1.15 ± 0.30, p for interaction = 2.14E-02).

Look up analysis

Leung et al. identified 241 CpG sites that were associated with PCB105 levels in a Faroe Island study [13]. We looked up the significance levels of these CpGs sites within our Σ35PCBs results and observed 23 CpGs in our analysis were nominally significant (p ≤ 5.00E-02) but none of these were significant after FDR correction (Table S2).

DMR analysis in ACHS-I

To determine if there were groups of CpGs (i.e. DMRs) that differed by PCB exposure, we grouped methylation data into PCB exposure tertiles and used the DMRcate analysis method. We identified 185 DMR/PCB exposure associations at FDR p ≤ 5.00E-02 (Table S3). Of the 115 unique FDR DMRs, only one DMR contained an ACHS-I Bonferroni significant individual CpG, GIMAP8 cg14251777. The 56 DMR/PCB associations with absolute ΔM ≥ 1.00% are shown in Table 4, representing 37 unique DMRs.

Table 4.

ACHS-I DMRs selected by FDR p ≤ 5.00E-02 and absolute differential methylation ≥ 1.00%.

Gene ID	DMR Coordinates	Number of CpGs	PCB Exposure	FDRp-valuea	Maximun ΔMb	Mean ΔMc	PCB/Dioxin Associationsd	DMR CpGs
C8orf31	chr8:144,120,106-144,120,681	9	Σ35PCBs	2.28E-06	1.99%	1.10%		cg07770222, cg22915785, cg08795786, cg00002224, cg22829325, cg08540654, cg03058664, cg21015266, cg03029255
FBLN2	chr3:13,590,221-13,590,503	10	Σ35PCBs	2.14E-04	4.43%	1.69%	TCDD	cg02027717, cg09543427, cg19392656, cg16604516, cg17561417, cg18603228, cg01896761, cg04195512, cg18406197, cg20605980
HLA-DPA1, HLA-DPB1	chr6:33,048,286-33,048,592	9	Σ35PCBs	2.89E-04	−4.50%	−3.47%	TCDD	cg09234582, cg25045942, cg13349035, cg02692313, cg10850215, cg03229061, cg17588455, cg26645432, cg06437840, cg14870156, cg19990651, cg01132696, cg14801692
LINC01833	chr2:45,159,894-45,160,389	6	Σ35PCBs	1.09E-04	4.87%	3.57%		cg06391932, cg02272993, cg23325963, cg08319079, cg24981593, cg01362243
LINC02036	chr3:193,921,489-193,922,370	6	Σ35PCBs	4.43E-06	−3.50%	−2.49%		cg24935598, cg10781513, cg27427514, cg21399510, cg01821429, cg00604425
LOC728392, NLRP1	chr17:5,403,516-5,404,337	8	Σ35PCBs	1.48E-05	−3.72%	−1.01%		cg00343839, cg16558432, cg01433610, cg04241623, cg08554115, cg12848098, cg17666974, cg22298430, cg27230784
SHANK1	chr19:51,171,712-51,172,144	5	Σ35PCBs	2.04E-05	2.88%	1.04%	PCB28, PCB52, PCB101, PCB126, PCB130, PCB153	cg01427575, cg00774550, cg10286380, cg06936631, cg18411150
SOX18	chr20:62,680,232-62,680,565	6	Σ35PCBs	1.60E-05	3.35%	1.38%	TCDD, Chlorodiphenyl (54% Cholrine)	cg06602546, cg04881720, cg08855301, cg23509896, cg12630336, cg18194354
TNFAIP8	chr5:118,603,742-118,603,938	5	Σ35PCBs	2.22E-04	−2.19%	−1.33%	TCDD	cg01810267, cg18689486, cg07376834, cg03723497, cg21130861
DDO	chr6:110,736,772-110,737,053	5	Mono-ortho PCBs	1.99E-03	−14.02%	−7.41%	TCDD	cg02872426, cg06413398, cg00804078, cg07164639, cg14956327
LINC02036	chr3:193,921,489-193,922,370	6	Mono-ortho PCBs	8.58E-06	−3.58%	−2.57%		cg24935598, cg10781513, cg27427514, cg21399510, cg01821429, cg00604425
LINGO3	chr19:2,291,373-2,291,872	6	Mono-ortho PCBs	3.34E-04	4.11%	2.74%		cg01294327, cg08565469, cg21869609, cg00835193, cg10068045, cg09370594
LMNB2, TIMM13,TMPRSS9	chr19:2,424,006-2,426,333	8	Mono-ortho PCBs	1.58E-05	3.42%	2.14%	TCDD	cg16170346, cg13496660, cg08569517, cg07418971, cg19497523, cg19906122, cg10341482, cg07593415
METTL22	chr16:8,735,575-8,736,433	5	Mono-ortho PCBs	3.92E-04	2.76%	1.15%		cg14384920, cg14513841, cg04858335, cg26821542, cg09085220
PSMA8	chr18:23,713,595-23,714,084	11	Mono-ortho PCBs	1.53E-06	4.15%	2.24%		cg11858305, cg24247743, cg25983544, cg01630687, cg06377543, cg03162994, cg00262344, cg15865827, cg22027766, cg01070760, cg21248196
BRF1, BTBD6	chr14:105,714,713-105,715,035	6	Di-ortho PCBs	6.90E-05	−3.55%	−1.05%	TCDD	cg11537379, cg09049585, cg15474367, cg16561543, cg02201215, cg01237565, cg14696445
C8orf31	chr8:144,120,106-144,120,681	9	Di-ortho PCBs	6.35E-06	2.07%	1.21%		cg07770222, cg22915785, cg08795786, cg00002224, cg22829325, cg08540654, cg03058664, cg21015266, cg03029255
EPHX2	chr8:27,348,165-27,348,453	5	Di-ortho PCBs	1.12E-04	−2.60%	−1.04%	TCDD, PCB81, PCB126, PCB153	cg15420071, cg25399743, cg14010279, cg21608792, cg22099723
GALR1	chr18:74,962,133-74,962,369	5	Di-ortho PCBs	1.34E-04	2.73%	2.01%	PCB28, PCB52, PCB101, PCB126, PCB130, PCB153	cg03502002, cg01178451, cg21938435, cg19056926, cg04534765
LINC02036	chr3:193,921,489-193,922,693	8	Di-ortho PCBs	9.15E-08	−3.40%	−2.27%		cg24935598, cg10781513, cg27427514, cg21399510, cg01821429, cg00604425, cg10148968, cg25548594, cg23003225, cg25862063
MAJIN	chr11:64,739,253-64,739,374	5	Di-ortho PCBs	5.63E-05	2.06%	1.01%		cg15681626, cg02700626, cg27639942, cg09278980, cg24408436
NRP2	chr2:206,628,415-206,628,773	9	Di-ortho PCBs	5.80E-05	−6.55%	−5.25%	TCDD, PCB126	cg01565529, cg04088940, cg22308949, cg05348875, cg14157435, cg10126788, cg20351668, cg25715429, cg10807027
NWD1	chr19:16,830,287-16,830,859	9	Di-ortho PCBs	7.50E-07	−9.20%	−5.51%		cg08260395, cg20398091, cg13208063, cg15845821, cg01649837, cg20249566, cg19784428, cg19344626, cg21035374
PHKG1	chr7:56,160,687-56,161,020	8	Di-ortho PCBs	1.12E-04	−4.10%	−2.65%	TCDD	cg19759064, cg15666071, cg26422861, cg17603988, cg08370546, cg21940640, cg13278478, cg15307593, cg09872392
SHANK1	chr19:51,171,712-51,172,144	5	Di-ortho PCBs	1.79E-05	2.81%	1.03%	PCB28, PCB52, PCB101, PCB126, PCB130, PCB153	cg01427575, cg00774550, cg10286380, cg06936631, cg18411150
TRIM6, TTRIM34	chr11:5,617,703-5,618,023	7	Di-ortho PCBs	4.86E-05	−2.30%	−1.49%	TCDD	cg19872996, cg14275946, cg14304349, cg00104484, cg00827519, cg06653507, cg22133704, cg15137954
ADORA2A-AS1, UPB1	chr22:24,890,690-24,890,936	8	Tri/tetra-ortho PCBs	1.81E-05	5.13%	3.51%	TCDD, PCB126	cg02109793, cg03364108, cg15552843, cg03503516, cg12492885, cg02469909, cg12550055, cg11210890
GALR1	chr18:74,961,724-74,962,794	19	Tri/tetra-ortho PCBs	7.15E-08	3.91%	1.41%	PCB28, PCB52, PCB101, PCB126, PCB130, PCB153	cg00854242, cg26721264, cg15146859, cg15343119, cg10486998, cg03175653, cg03032214, cg10122698, cg03659519, cg20872937, cg17911318, cg03502002, cg01178451, cg21938435, cg19056926, cg04534765, cg06360427, cg02611744, cg10390058
IFT172, KRTCAP3, NRBP1	chr2:27,664,918-27,665,711	10	Tri/tetra-ortho PCBs	5.48E-06	5.27%	3.31%	TCDD, PCBs	cg26350635, cg18428193, cg04845466, cg24768116, cg12000995, cg12648201, cg21248554, cg17158414, cg02592271, cg11618577, cg20102877, cg26034919
LAD1, TNNI1	chr1:201,368,390-201,369,031	11	Tri/tetra-ortho PCBs	4.19E-05	3.31%	1.16%	TCDD	cg19713196, cg02838178, cg10095226, cg14952599, cg07561547, cg21964551, cg27425612, cg23291161, cg25947945, cg23715582, cg26539023
LOC105369632, SPSB2, TPI1	chr12:6,982,384-6,982,548	6	Tri/tetra-ortho PCBs	9.89E-05	−2.13%	−1.25%	TCDD, PCB81, PCB153	cg13449363, cg14553260, cg22985390, cg22855052, cg06176834, cg26309352, cg00613384
LOC107985544, MORC2, MORC2-AS1	chr22:31,318,103-31,318,546	9	Tri/tetra-ortho PCBs	4.58E-06	2.31%	1.75%		cg05450667, cg17468100, cg08837037, cg15425811, cg07459594, cg16210088, cg20548231, cg17940740, cg13896476
MEIOC	chr17:42,733,527-42,733,994	10	Tri/tetra-ortho PCBs	3.72E-05	3.62%	1.77%		cg08496953, cg22605919, cg12892303, cg17269633, cg15000379, cg02932314, cg12127472, cg12781700, cg22190438, cg21095561
NWD1	chr19:16,830,287-16,830,859	9	Tri/tetra-ortho PCBs	3.52E-05	−9.32%	−6.13%		cg08260395, cg20398091, cg13208063, cg15845821, cg01649837, cg20249566, cg19784428, cg19344626, cg21035374
RBM46	chr4:155,702,172-155,703,138	13	Tri/tetra-ortho PCBs	4.53E-06	3.24%	1.73%		cg10920378, cg02992596, cg10316270, cg06318837, cg26413432, cg22103164, cg08282428, cg23483495, cg01827861, cg22496683, cg01289218, cg20745987, cg19340941, cg05914674, cg00137234
RP11-66B24.2, RP11-66B24.7	chr15:101,389,272-101,390,350	13	Tri/tetra-ortho PCBs	7.76E-09	4.37%	1.66%		cg17221377, cg25878441, cg16548362, cg07882398, cg23117796, cg07035436, cg05000474, cg18304498, cg09785344, cg05500125, cg10405604, cg04392029, cg13494481
RPL37, SNORD72	chr5:40,835,533-40,836,182	7	Tri/tetra-ortho PCBs	3.12E-05	−2.51%	−1.24%	TCDD, PCB126	cg11110213, cg01087697, cg27450992, cg04986097, cg05478818, cg04002187, cg17351974, cg20777740
SEMA3B, SEMA3B-AS1	chr3:50,304,930-50,305,083	5	Tri/tetra-ortho PCBs	2.97E-05	−2.03%	−1.57%	TCDD	cg12999941, cg11039694, cg09034667, cg06192381, cg16216311
SNORD36	chr13:23,309,689-23,310,225	7	Tri/tetra-ortho PCBs	1.44E-05	3.82%	2.66%		cg26361286, cg15973954, cg01863042, cg08083251, cg03042692, cg05215994, cg20395040
B4GALNT2	chr17:47,209,546-47,210,166	10	Oestrogenic PCBs Group 1	2.64E-06	3.60%	1.57%		cg11383850, cg19687959, cg24238356, cg25446098, cg09367266, cg20233029, cg01688264, cg01147550, cg01429662, cg22259831
DPPA5	chr6:74,063,982-74,064,434	8	Oestrogenic PCBs Group 1	5.98E-05	4.02%	1.76%		cg01652244, cg15414745, cg23414861, cg18052665, cg09071762, cg08808852, cg02931604, cg14276083
NCL	chr2:232,348,334-232,348,794	5	Oestrogenic PCBs Group 1	7.56E-05	3.97%	2.63%		cg03727500, cg20788612, cg05868531, cg15371801, cg11559198
PSMA8	chr18:23,713,595-23,714,084	11	Oestrogenic PCBs Group 1	4.74E-07	4.61%	2.67%		cg11858305, cg24247743, cg25983544, cg01630687, cg06377543, cg03162994, cg00262344, cg15865827, cg22027766, cg01070760, cg21248196
B4GALNT2	chr17:47,209,546-47,210,166	10	Oestrogenic PCBs Group 1	2.64E-06	3.60%	1.57%		cg11383850, cg19687959, cg24238356, cg25446098, cg09367266, cg20233029, cg01688264, cg01147550, cg01429662, cg22259831
DPPA5	chr6:74,063,982-74,064,434	8	Oestrogenic PCBs Group 1	5.98E-05	4.02%	1.76%		cg01652244, cg15414745, cg23414861, cg18052665, cg09071762, cg08808852, cg02931604, cg14276083
NCL	chr2:232,348,334-232,348,794	5	Oestrogenic PCBs Group 1	7.56E-05	3.97%	2.63%		cg03727500, cg20788612, cg05868531, cg15371801, cg11559198
PSMA8	chr18:23,713,595-23,714,084	11	Oestrogenic PCBs Group 1	4.74E-07	4.61%	2.67%		cg11858305, cg24247743, cg25983544, cg01630687, cg06377543, cg03162994, cg00262344, cg15865827, cg22027766, cg01070760, cg21248196
B4GALNT2	chr17:47,209,546-47,210,368	14	Oestrogenic PCBs Group 2	4.01E-06	3.25%	1.36%		cg11383850, cg19687959, cg24238356, cg25446098, cg09367266, cg20233029, cg01688264, cg01147550, cg01429662, cg22259831, cg16206138, cg19133221, cg18208707, cg03167683
B4GALNT2	chr17:47,209,546-47,210,368	14	Oestrogenic PCBs Group 2	4.01E-06	3.25%	1.36%		cg11383850, cg19687959, cg24238356, cg25446098, cg09367266, cg20233029, cg01688264, cg01147550, cg01429662, cg22259831, cg16206138, cg19133221, cg18208707, cg03167683
ZNF875	chr19:37,825,009-37,826,008	12	Oestrogenic PCBs Group 2	3.57E-07	9.83%	7.41%		cg12147799, cg10237978, cg26734888, cg14166009, cg13687570, cg05280698, cg24834889, cg12948621, cg08565796, cg23756236, cg12024906, cg23448505
LINGO3	chr19:2,291,373-2,291,872	6	Antiestrogenic PCBs	5.13E-05	4.22%	2.89%		cg01294327, cg08565469, cg21869609, cg00835193, cg10068045, cg09370594
LINGO3	chr19:2,291,373-2,291,872	6	Antiestrogenic PCBs	5.13E-05	4.22%	2.89%		cg01294327, cg08565469, cg21869609, cg00835193, cg10068045, cg09370594
PSMA8	chr18:23,713,595-23,714,084	11	Antiestrogenic PCBs	6.96E-08	4.18%	2.24%		cg11858305, cg24247743, cg25983544, cg01630687, cg06377543, cg03162994, cg00262344, cg15865827, cg22027766, cg01070760, cg21248196
PSMA8	chr18:23,713,595-23,714,084	11	Antiestrogenic PCBs	6.96E-08	4.18%	2.24%		cg11858305, cg24247743, cg25983544, cg01630687, cg06377543, cg03162994, cg00262344, cg15865827, cg22027766, cg01070760, cg21248196
RBM46	chr4:155,702,069-155,703,070	13	Antiestrogenic PCBs	1.19E-05	3.49%	1.77%		cg02201858, cg10920378, cg02992596, cg10316270, cg06318837, cg26413432, cg22103164, cg08282428, cg23483495, cg01827861, cg22496683, cg01289218, cg20745987, cg19340941, cg05914674
RBM46	chr4:155,702,069-155,703,070	13	Antiestrogenic PCBs	1.19E-05	3.49%	1.77%		cg02201858, cg10920378, cg02992596, cg10316270, cg06318837, cg26413432, cg22103164, cg08282428, cg23483495, cg01827861, cg22496683, cg01289218, cg20745987, cg19340941, cg05914674

DMRs are listed in order of PCB exposure and GeneID. DMRs were identified using DMRcate v.1.18.0, comparing the highest v. lowest tertiles for each exposure.

aFDR p-value for the entire DMR. Model was adjusted for age, race, sex, smoking status, total serum lipids, bisulphite-conversion batch, and estimated percentages of CD4+ and CD8 + T-cells, CD19 + B-cells, monocytes, granulocytes, and natural killer cells.

bThe differential methylation value (ΔM) for the DMR CpG with the largest absolute differential value. Differential methylation was calculated by subtracting methylation in lowest tertile from the highest.

cThe average differential methylation value across the entire DMR. Differential methylation was calculated by subtracting methylation in lowest quartile from the highest.

dGenes associated with PCB or dioxin exposure based on the 2019 Comparative Toxicogenomics Database.

Functional analysis and Comparative Toxicogenomics Database

The GREAT functional analysis tool was used to find possible functional enrichment of molecular signatures of the FDR significant CpGs or DMRs [45]. We observed enrichment of five related immunological molecular signatures (Table S5) in the Molecular Signatures Database [46].

We also searched the 2019 Comparative Toxicogenomics Database (CTD) [47] to determine if CpG (Table 3 and S1) and DMR genes (Table 4 and S2) had been reported in the literature to interact with either PCB or dioxin exposure in experimental studies. For individual CpGs, 13 of the 17 (76%) unique genome-wide significant genes and 199 (70%) of the 286 FDR significant unique genes were listed as interacting with PCB or dioxin exposures. For ACHS-I 37 unique DMRs in Table 4, 19 (51%) of the genes were associated with either PCB or dioxin exposure based on retrieved curated chemical–gene association.

Discussion

The subgroup of participants in ACHS-I and II that were analysed for DNA methylation levels in whole blood DNA had very similar PCB exposure distributions relative to age, race and sex as have been reported [1,48]. Our initial hypothesis for these analyses was that methylation in the AHRR gene, known to be regulated by AHR, might be strongly affected by PCB exposure. Interestingly, while smoking had the expected profound impact on AHRR methylation levels in ACHS-I (cg05575921 p = 7.27E-269), as has been reported many times [17-19,39-42], we detected no association between PCBs (potential AHR ligands) and AHRR methylation in non-smokers (p = 9.40E-01).

Many studies of PCBs and dioxins clearly implicate AHR as mediating the adverse outcomes associated with exposure in model systems and humans [15]. Some PCBs and dioxins interact with AHR in terms of the planar orientation of their two phenyl groups. In dioxin, the two phenyl groups are in a fixed coplanar orientation, which makes it a strong ligand for AHR. The phenyl groups in non-ortho and mono-ortho substituted PCBs can assume a coplanar configuration, allowing these PCBs to have dioxin-like properties including being an AHR ligand. In di-ortho and tri/tetra-ortho substituted PCBs, the chlorines in the ortho positions force the two phenyl groups out of coplanar conformation, leading to non-dioxin-like properties [49]. It should be noted that the most dioxin-like PCBs (non-ortho substituted) were not measured in ACHS-I. Based on our original hypothesis, we expected to observe the most significant changes in methylation among the PCB exposure group that was the most dioxin-like i.e., the mono-ortho substituted PCBs. However, in ACHS-I, exposure to the least dioxin-like PCBs (tri/tetra-ortho substituted) produced the largest number of significant associations (229 CpGs, 82% of all associations). This finding suggests that most of the effects of PCB exposure on whole blood DNA methylation that we observe in this study are not likely to be mediated through the AHR pathway. However, only nine of the 35 PCBs in these analyses were mono-ortho substituted, and these PCBs are weak AHR ligands.

Other studies have reported associations between PCB or dioxin exposures and altered DNA methylation, although most examined global measures of methylation or specific candidate genes. Rusiecki et al. observed that in the Arctic Monitoring and Assessment Program, higher serum levels of 11 PCBs (99, 105, 118, 128, 138, 153, 156, 170, 180, 183, 187) in 71 Greenlandic Inuits were related to global hypomethylation [11]. Two studies found that TCDD [50] and PCB153 [51] exposure were predominantly associated with hypomethylation in sperm DNA. In a study of 368 cancer-free Koreans, researchers found a non-monotonic association between PCB105, PCB138, and PCB153 and MGMT promoter hypomethylation [52]. In a Swedish cross-sectional cohort of 524 elderly men and women, two POPs (PCB126 and OCDD) were connected to global DNA hypermethylation [53]. The research from van den Dungen et al. found in a cohort of Dutch men that serum levels of dioxins and the sum of seven PCBs (28, 52, 101, 118, 138, 180) were significantly associated with primarily hypomethylated DMRs in multiple genes with an average change of 7.4% between the highest and lowest exposed subjects [10].

Some studies have observed possible sex-specific effects. In a subset of a prospective Japanese birth cohort (the Hokkaido Birth Cohort Study on Environment and Children’s Health, n = 169), cord blood methylation levels for H19 and LINE1 were hypermethylated and related to PCB170, PCB 178, PCB180, and PCB183 exposure with the suggestion that this association might be more pronounced in females than males [14]. In the only study utilizing the 450K array, Leung et al. [13] examined a prospective birth cohort of 73 Faroe islanders with high levels of PCB exposure at birth and found sex-specific methylation associations with PCB105 levels. Our look-up analysis of the 195 CpGs reported to be FDR significant in the Leung et al. cord blood study identified 23 CpGs with nominal p-value 5.00E-02, but none reached FDR-corrected significance level. The lack of overlap between a cord blood study with in utero exposures with the present study of older adults is not surprising. Among our most significant CpGs, we tested for sex-specific effects in ACHS-I and observed potential sex-related effect modification for AC068134.2 cg21566642 and GALNT2 cg04991747; however, the size of the ACHS-I cohort (n = 518) limited the statistical power to detect true effect-measure modification.

Among the ACHS-I Bonferroni significant CpGs, the most differentially methylated was PRSS23 cg00475490 (ΔM −5.02 ± 0.60 in di-ortho substituted PCBs), and it showed a strong hypomethylation trend from low to high exposure when examined by quartiles. PRSS23 cg00475490 was the only CpG that displayed persistent and significant alteration in ACHS-II. PRSS23 was first identified as an ovarian protease [54] and has been shown to be upregulated by oestrogen receptor alpha and stimulate growth in the MCF-7 breast cancer cell line [55]. In the Comparative Toxicogenomics Database, it is also associated with dioxin exposures in model systems. Among blood cells, it is expressed only in natural killer cells (Figure 5). The presence of a functional NRF2/MAF binding site in PRSS23 enhancer may indicate that this gene is regulated by NRF2 (NFE2L2), the master transcription factor regulator of the response to oxidative stress, which also modulates inflammatory and metabolic pathways [56]. GALNT2 has been associated with triglyceride levels [57] and may be associated with Type 2 diabetes [58]. The deubiquitinase, OTULIN, was also highly significant in ACHS-I, is an essential negative regulator of the inflammatory transcript factor nuclear factor κB and plays a role in immune homoeostasis [59]. Using genomic location of the FDR significant CpGs and the GREAT functional annotation program we identified a group of immunological molecular signatures. It is established that PCBs can have effects on the immune system [60-62], and our data suggest that PCB exposures may be perturbing the epigenetic profile of immune system genes.

Our study observed 28 CpG/PCB associations (17 unique genes) at genome-wide significance, 369 associations in 286 CpGs at FDR significance, and 115 genes with significant DMRs. Of these genes a large proportion [14/17 (76%) Bonferroni significant genes; 199/286 (70%) FDR significant genes; and 78/115 (68%) FDR significant DMR genes] had a reported interaction with PCB or dioxin in various model system studies reported in the Comparative Toxicogenomics Database. However, surprisingly none of the canonical AHR pathway genes were affected, e.g. genes in the CYP, GST and UGT families.

Two factors that were limitations of our study were the relatively limited genome-wide coverage of the EPIC array (~3% of all genomic CpGs) and the relatively small sample sizes (ACHS-I n = 518 and ACHS-II n = 299) in terms of multiple comparisons correction. The main strength of our study was the opportunity to examine effects in ACHS-I and then later in ACHS-II. We observed half of the ACHS-I genome-wide associations displayed similar directional effect in ACHS-II and only one CpG (PRSS23 cg00475490) among ACHS-I top hits, displayed a persistent, significant association in ACHS-II. The factors that may have contributed to the limited replication in ACHS-II include the smaller sample size (n = 299), the reduced exposure levels (~10% lower), and potential sampling bias (survival bias). In both ACHS-I and ACHS-II individuals recruited long after the occurrence of very high exposures could be more resilient to adverse health effects. From a public health perspective, our observation of reduced PCB effects on methylation after seven years can be considered a reassuring outcome as opposed to a lack of replication. In future research, we will explore methylation associations with dioxins, furans and non-ortho PCBs that were only measured in ACHS-II. We will also examine the associations among methylation, POPs, and adverse health outcomes, e.g. hypertension, liver disease, diabetes, and metabolic syndrome in both cohorts.

We have conducted the first EWAS exploring the associations between PCB exposure and DNA methylation levels in adults and identified numerous highly significant CpGs and DMRs associated with exposure. Our study design has allowed us to examine effects of PCB exposure on DNA methylation at two points in time, seven years apart. It is unclear how biological important these altered methylation levels are, or if CpGs such as PRSS23 cg00475490 are useful biomarkers of exposure. However, it is encouraging that PCB effects on methylation generally appear to be attenuated in ACHS-II after seven years.