Exposure to polybrominated biphenyl (PBB) associates with genome-wide DNA methylation differences in peripheral blood

ABSTRACT

In 1973, Michigan residents were exposed to polybrominated biphenyl (PBB) when it was accidentally added to farm animal feed. Highly exposed individuals and their children have experienced endocrine-related health problems, though the underlying mechanism behind these remains unknown. We investigated whether PBB exposure is associated with variation in DNA methylation in peripheral blood samples from 658 participants of the Michigan PBB registry using the MethylationEPIC BeadChip, as well as investigated what the potential function of the affected regions are and whether these epigenetic marks are known to associate with endocrine system pathways. After multiple test correction (FDR <0.05), 1890 CpG sites associated with total PBB levels. These CpGs were not enriched in any particular biological pathway, but were enriched in enhancer and insulator regions, and depleted in regions near the transcription start site or in CpG islands (p < 0.05). They were also more likely to be in ARNT and ESR2 transcription factor binding sites (p = 3.27e-23 and p = 1.62e-6, respectively), and there was significant overlap between CpGs associated with PBB and CpGs associated with estrogen (p < 2.2e-16). PBB-associated CpGs were also enriched for CpGs known to be associated with gene expression in blood (eQTMs) (p < 0.05). These eQTMs were enriched for pathways related to immune function and endocrine-related autoimmune disease (FDR <0.05). These results indicate that exposure to PBB is associated with differences in epigenetic marks that suggest that it is acting similarly to estrogen and is associated with dysregulated immune system pathways.

KEYWORDS:

• Endocrine-disrupting compound
• epigenome-wide association study (EWAS)
• DNA methylation
• polybrominated biphenyl
• environmental health
• epigenetics

Introduction

Endocrine-disrupting compounds (EDCs) are a diverse class of chemicals that interfere with the secretion, transportation, binding, function, or synthesis of hormones [1]. EDCs have been found in pesticides, plastics, electronics, and household dust, and as such, they are ubiquitous in industrialized countries [1–4]. As the number of chemicals classified as EDCs is increasing, EDCs are being linked with long-term health consequences. Animal models and human cohorts have associated increasing exposure to EDCs with decreased fertility, altered thyroid hormone levels, and increased risk of cancer in people directly exposed, as well as increased risk of adverse health outcomes in the children of those that were directly exposed [5–16].

To further study the effect that EDCs have on humans, we utilized the Michigan Polybrominated Biphenyl (PBB) Registry, which was recruited after an agricultural accident introduced an EDC, PBB, into the food supply in the 1970’s [17,18]. An estimated 6 to 8 million people living in Michigan were exposed to PBB primarily through the consumption of contaminated animal products, and because PBB is both lipophilic and biologically stable, people exposed over 40 years ago still have detectable PBB levels in their blood [17–22]. Additionally, because PBB can pass through the placental barrier and is present in breast milk, children born to exposed women are also exposed [23]. This cohort has reported many endocrine-related and reproductive health concerns, such as thyroid disease, spontaneous abortions, earlier age of menarche, and genitourinary conditions in male offspring, as well as breast cancer, digestive cancer, and lymphomas [24–33]. In general, the people who were exposed in utero or in childhood were found to be particularly susceptible to endocrine-related conditions [34].

The mechanisms behind the increased health risks with exposure are not yet fully understood, nor is it known what the health risks are for future generations. One potential mechanism is through EDCs binding steroid hormone receptors, which could alter gene expression levels and epigenetic marks, and thus impact health outcomes [35–38]. In animal models and cell lines, exposure to polychlorinated biphenyl (PCB) and polybrominated diphenyl ethers (PBDEs), EDCs with a chemical structure similar to PBB, is associated with global differences in DNA methylation and with alterations of DNA methylation and histones at particular locations in the genome [39–41]. In cell lines, alterations of epigenetic marks from PCB exposure has led to increases in inflammation and adipocyte differentiation [39,41]. In human cohorts, PCB exposure has been associated with differences in global DNA methylation levels, although there may be population-specific differences that have not been resolved. In adult populations from Greenland and Korea, increased exposure to PCB was associated with decreased global DNA methylation [42,43]. However, in an older, Swedish cohort, exposure to PCB was associated with increased global DNA methylation [44].

No studies have evaluated whether PBB exposure associates with epigenetic marks, and epigenome-wide studies of other EDCs similar in chemical structure to PBB are lacking. This study used Illumina’s Infinium MethylationEPIC BeadChip to interrogate the association between DNA methylation and exposure to PBB in the members of the Michigan PBB Registry. We then used existing databases to determine if the CpGs associated with PBB exposure were enriched in any particular biological pathways or genomic regions. The results will improve our understanding of how EDCs like PBB influence the epigenome and the health risks of people exposed and their children.

Results

The study population is composed of 658 people of primarily White/Non-Hispanic ancestry who were exposed to PBB (range: 0.01–236.73 ppb; Table 1). There were more female participants than male participants (N = 381 vs. 277). The cohort’s racial composition is indicative of the population structure of rural Michigan in the 1970s. Higher age was associated with increased current levels of PBB (R = 0.25; p = 3.96e-11), and males had significantly higher exposure than females (p = 4.55e-07). Because a majority of the cohort was exposed during the year-long exposure incident, age at exposure is highly correlated with current age (R = 0.98; p < 2.2e-16). Current PBB levels were not associated with their current lipid level (R = −0.07; p = 0.06).

Table 1. Characteristics of subset of Michigan PBB registry that was selected for epigenetic analysis.

Cohort Demographics
N	658
Percent male^a	277 (42.1%)
Current Age (years)^b	54.3 ± 12.7
Age when exposed (years)^b	15.2 ± 11.5
Total PBB level (ppb)^c	0.48 (4.71)
Total PBB level (ng per g lipid)^c	73.08 (4.94)
Race/Ethnicity ^a
White/Non-Hispanic	638 (97%)
White/Hispanic	20 (3%)

^aFrequency and percentage

^bMean and standard deviation

^cGeometric mean and geometric standard error

Association of DNA methylation and PBB exposure

Global DNA methylation level was negatively associated with total PBB, controlling for age, sex, lipid level, and cell type proportions (p = 0.04; Figure S3). When the association between the total PBB exposure level and the methylation proportion at each of the 816,999 CpG sites was calculated individually, 1890 CpGs were associated with total PBB (Figure 1, Table S1). Of these CpGs, 83.91% had lower methylation levels in those with higher levels of PBB (Figure S4). Results of the analysis adjusting for total lipid level were highly correlated with results obtained without adjustment for lipid level (R = 0.99; p < 2.2e-16; Figure S5). To ensure that summing the congeners to total PBB did not alter the test statistics in this analysis, the test statistics from the total PBB analysis were compared to analyzes of each congener individually. The results were all highly correlated (R = 0.78–0.99), and indicated that analyzing each congener individually would lead to inflated test statistics (Figure S6 and S7). Of the top 20 significant results, fourteen are annotated to genes, including nuclear hormone receptors (cg04664328, p = 6.52E-09) or linked to endocrine-related diseases (cg13991617, p = 2.09E-10; cg08800183, p = 1.16E-07), and some are involved in immune response (cg00157963, p = 6.90E-08; cg20143841, p = 8.01E-08). Many were also members of the cytochrome P450 family (cg02879655, p = 9.11e-6; cg07148145, p = 3.94e-5; cg05372113, p = 8.37E-05; cg19201144, p = 1.05E-04). Scatterplots for the most significant CpGs were plotted (Figure S8). These plots showed that several samples had lower DNA methylation levels at some of the CpGs compared to the rest of the samples. As a sensitivity analysis to evaluate whether these samples influenced the results from the overall analysis, these samples (N = 5) were dropped and the analysis was repeated. The test statistics from these analyzes were highly correlated (R = 0.95, p < 2.2e-16), especially in the 1890 significant CpGs (R = 0.98, p < 2.2e-16) suggesting that the inclusion of these samples did not bias the results (Figure S9). There was no enrichment of KEGG pathways among PBB-associated CpGs [45].

Figure 1. Total PBB level is associated with DNA methylation differences genome-wide. A Manhattan plot of the association of total PBB level with DNA methylation proportion at 816,999 sites. The x-axis is the location of each site across the genome. The y-axis is the – log₁₀ of the p-value for the association with PBB. The blue line indicates statistical significance (FDR < 0.05). 1890 CpGs had a statistically significant association with total PBB level.

Enrichment of PBB-associated CpGs in positional regions

Because there were no specific biological pathways enriched, we tested whether PBB-associated CpGs were overrepresented in the regions near CpG islands or genes. In regions annotated to gene position, PBB-associated CpGs were depleted in the first exon (OR = 0.49; 95% CI = 0.37–0.65; p = 2.85e-7), and near the transcription start site (within 200 base pairs: OR = 0.56; 95% CI = 0.46–0.68; p = 7.27e-9; within 1500 base pairs: OR = 0.83; 95% CI = 0.72–0.95; p = 0.01), and were enriched in regions in the gene body (OR = 1.12; 95% CI = 1.02–1.23; p = 0.01; Figure 2A). In regions annotated to CpG islands, PBB-associated CpGs were depleted in the CpG island (OR = 0.34; 95% CI = 0.28–0.40; p < 2.2e-16) and were enriched in the south shelf (OR = 1.44; 95% CI = 1.17–1.78; p = 0.0006; Figure 2B).

Figure 2. PBB-associated CpGs are depleted in regions near transcription start sites. Enrichment tests were conducted to determine if PBB-associated CpGs were enriched in certain regions of genes (Part A), or in certain regions of CpG island (Part B). PBB-associated CpGs are depleted in regions near transcription start sites and CpG islands, but are enriched in regions within the gene body and in the south shelf of a CpG island. * indicates statistical significance (p < 0.05).

Enrichment of PBB-associated CpGs in functional regions

Because CpGs in the gene body and south shelf can have a variety of functions [46,47], CpGs were annotated to functional regions (as defined by ChromHMM) in order to better understand the effects of PBB on the epigenome. PBB-associated CpGs were depleted in active promoters (OR = 0.76; 95% CI = 0.67–0.87; p = 4.78e-5), inactive promoters (OR = 0.56; 95% CI = 0.41–0.76; p = 0.0002), and heterochromatin regions (OR = 0.60; 95% CI = 0.52–0.69; p = 2.90e-13), and were enriched in strong enhancers (OR = 3.37; 95% CI = 2.88–3.94; p < 2.2e-16), weak enhancers (OR = 2.00; 95% CI = 1.73–2.32; p < 2.2e-16), and insulator regions (OR = 1.39; 95% CI = 1.13–1.72; p = 0.001; Figure 3). To better assess the biological meaning of these enhancers, enrichment for transcription factor binding sites near the PBB-associated CpGs was also tested. Many transcription factor binding sites (N = 56) were enriched or depleted in the results (Table S2). Some of these enriched transcription factors were involved in the AHR pathway (HIF1A::ARNT – p = 2.09e-37; ARNT – p = 3.27e-23), sex-determining pathways (SRY – p = 7.57e-18; SOX2 – p = 2.54e-26), estrogen-related pathways (ESR2 – p = 1.62e-06; STAT3 – 3.00e-10; PPARG::RXRA – p = 8.26e-12), or nuclear hormone receptor signaling (Myc – p = 6.09e-43; RORA_1 – p = 2.33e-02).

Figure 3. PBB-associated CpGs are enriched in enhancer and insulator regions. Enrichment tests showed that PBB-associated CpGs were enriched in certain functional regions in the genome. Specifically, PBB-associated CpGs were more likely to be in enhancers regions (both strong and weak) and insulator regions, and were less likely to be in heterochromatic regions or areas near the promotors of both active and inactive genes. * indicates statistical significance (p < 0.05).

PBB-associated CpGs correlated with estrogen signal

Because PBB is an endocrine-disrupter and ESR2 transcription factor binding sites were enriched in the results, it was of interest whether the epigenetic signal that associated with PBB was similar to the epigenetic signal from natural hormones. Complementary to estrogen related transcription factor binding sites being enriched near the PBB-associated CpGs, the effect size for CpGs associated with estradiol levels is correlated with the effect size of those CpGs with PBB, both in the total cohort (R = 0.20, p < 2.2e-16) and in the female participants (N = 381, R = 0.16, p < 2.2e-16. Figure S10) [48]. Additionally, of the 700 PBB-associated CpGs that were also in the estradiol analysis, 211 of them of them were also significantly associated with estradiol. This is considerable enrichment over what was expected by chance (OR = 6.50; 95% CI = 5.53–7.64; p < 2.2e-16), suggesting that many of the epigenetic marks that associate with estradiol levels also associated with PBB levels.

PBB-associated CpGs are enriched for immune-related gene expression

Because many of the PBB-associated CpGs were in enhancer regions, it was of interest to investigate which genes may be regulated by those CpGs, even in the absence of RNA data. Therefore, the PBB-associated CpGs that were also on the Infinium HumanMethylation450 BeadChip (N = 735) were compared to previously published sets of CpGs that associate with gene expression (eQTMs) [49]. PBB-associated CpGs were more likely to be eQTMs in both datasets (MESA: OR = 2.54; 95% CI = 2.19–2.93; p < 2.2e-16; GTP: OR = 1.45; 95% CI = 1.15–1.82; p = 0.001; Figure 4, Table S3). The genes whose expression associated with these PBB-associated eQTMs were also analyzed (the genes in 5,042 CpG-gene expression pairs). Functional annotation with DAVID [50,51] showed that these genes were enriched for several immune functions, including autoimmune thyroid disease (Table 2). Because of the enrichment for immune functions and autoimmune pathways in the eQTMs, the epigenetic signal from PBB was compared to the other epigenetic immune phenotypes. The effect size for CpGs associated with C-Reactive Protein (CRP) levels, a marker of low-grade inflammation, was correlated with the effect size of those CpGs with PBB (R = 0.30, p-value = 9.58e-06, Figure S11) [52]. However, while PBB was associated in univariate analyzes with CD8T, CD4T, NK, and monocyte proportions, after controlling for age, sex, and lipid levels, PBB was only associated with the proportion of B cells (p = 0.02; Table S4). PBB was also not associated with other markers of immune function and inflammation like the neutrophil to lymphocyte (NLR) ratio (p = 0.36) and CD4T to CD8T ratio (p = 0.42) calculated from the cell types called from the methylation data (Table S4) [53].

Table 2. Pathways enriched in PBB-associated eQTMs. PBB-associated CpGs associate with gene expression (eQTMs), and expression of these genes were more likely to be located in KEGG pathways implicated in immune function and autoimmune disorders (FDR < 0.05).

KEGG Pathway	FDR
Hematopoietic cell lineage	0.0046
Graft-versus-host disease	0.0091
Allograft rejection	0.012
Type I diabetes mellitus	0.015
Antigen processing and presentation	0.017
Asthma	0.018
Primary immunodeficiency	0.028
Autoimmune Thyroid Disease	0.029

Figure 4. PBB-associated CpGs are more likely to be CpGs linked to expression. PBB-associated CpGs were more likely to be eQTMs than expected by chance (identified in blood in two datasets: GTP and MESA). * indicates statistical significance (p < 0.05).

Discussion

In this study we investigated the associations between PBB and the epigenome by conducting a EWAS in a cohort that was uniquely exposed to PBB through contamination of the food supply in Michigan in the 1970’s. We then investigated the role of the associated epigenetic marks by testing for enrichment in biological pathways, positional regions, functional regions, locations linked to gene expression, and transcription factor binding sites, as well as investigating the correlates between the signal from PBB and the signal from estrogen and inflammatory markers. We found that 1890 CpGs across the genome associated with higher current PBB levels. A majority (83.91%) had lower methylation proportions with higher PBB levels. This is consistent both with our finding that global methylation was lower in people with higher PBB levels, and with a majority of the previous findings in the structurally-related EDC, PCB, where two studies found that higher PCB was associated with global hypomethylation (as measured using ALU and LINE1 methylation) [42,43]. Our finding is only inconsistent with a study that found that higher PCB was associated with higher global methylation (as measured using Luminex) [44].

While no biological pathways were enriched among the PBB-associated sites, it is unclear whether this is because annotation databases do not take into account the function of CpGs and that enhancers may be acting on distant genes. This means that currently CpGs in enhancers are either not annotated to any gene (for CpGs not within a gene) or are incorrectly annotated to a gene (for CpGs that are within one gene but increase the transcription of a different gene). This is especially a concern because the MethylationEPIC array added 350,000 probes for known enhancers, and a majority of our results (61.1%) were for probes that were not on the previous array. Therefore, it is possible, as more is known about the function of these CpGs and which genes these CpGs enhance the transcription of, that specific pathways would be enriched. That there were immune-related pathways enriched in the genes from the eQTMs supports that these PBB-associated CpGs do have a function, but that appropriate annotation may be lacking. However, it is also possible that PBB does not act on a specific biological pathway and that instead it is acts more broadly on gene transcription and cellular function.

We found that PBB-associated CpGs were less likely to be in positional regions near the transcription start site or in CpG islands, but were more likely to be in the body of a gene. The function of gene body methylation is often unclear [46,47,54,55], however, functional enrichment showed PBB-associated CpGs were more likely to be in enhancer and insulator regions, and similar to positional enrichment, were less likely to be in promoter or heterochromatin regions. This could be because methylation at promoters is more commonly determined by genetic factors and thus is protected [56,57], and that enhancers are more vulnerable to environmental exposures [57]. However, this is difficult to assess, given that enhancers were only interrogated by the most recent MethylationEPIC array, and previous genome-wide studies with EDCs have used the Human Methylation450 array that has fewer probes in enhancer regions [7,58,59]. It is possible that as more studies with environmental exposures are done on arrays that more completely assess enhancer methylation, enrichment at enhancer regions will be more commonly found. This analysis itself is also limited because we could only test for enrichment of areas probed by the array, and, as more of the epigenome is assessed, other functional regions could be enriched.

PBB-associated CpGs were more likely to be in transcription factor binding sites for a member of the AHR pathway, and many of the significantly associated CpGs were in members of the cytochrome P450 family. This finding is consistent with early studies in rat models that found that PBB increased activity of other members of the AHR pathway like the aryl hydrocarbon hydroxylase (AHH) and members of the cytochrome P-450 enzyme family [60–62]. Additionally, PBB-associated CpGs were more likely to be near nuclear hormone receptor transcription factor binding sites. This is consistent with the altered hormone metabolism seen in rat models exposed to PBBs [63–65]. PBB-associated CpGs also massively overlap with estradiol-associated CpGs, and the effect sizes of CpGs associated with estradiol is correlated with the effect sizes seen with PBB. PBB exposure being associated with estrogen-related transcription factor binding sites and with the epigenetic signal from estrogen may help explain the epidemiological findings in this cohort, which found an increase in genitourinary conditions in males exposed in utero through PBBs crossing the placental barrier [29], and an early age of menarche and an increased risk of spontaneous abortion in women exposed in utero [30,31], all of which are influenced by estrogen levels [66–68]. This could potentially indicate that PBB is acting as a weak estrogen on the epigenome and increasing exposure affects people’s risk for endocrine related disorders.

Consistent with the association of PBB-associated CpGs with enhancers, we also found that our results significantly overlapped with CpGs that were known eQTMs in two separate datasets, and those eQTMs were associated with genes involved in immune function and autoimmune disease. This is consistent with early studies that that showed that 35–40% of Michigan farmers who had eaten PBB-contaminated food had a significant reduction in in vitro immune function [69]. While no studies have looked at PBB’s association with the development of autoimmune disorders, studies in this cohort have found PBB exposure is associated with the development of hypothyroidism (which is commonly caused by an autoimmune disorder) in both men and women, although neither of these analyzed hypothyroidism by type [24,25]. Additionally, when the signal for the CpGs associated with low grade inflammation was compared to the signal from PBB, they were correlated, and PBB was associated with the proportion of B cells, even after controlling for age, sex, and lipid levels. This indicates both that the PBB-associated CpGs do have a connection to gene expression and immune function. However, a limitation was that the GTP eQTMs were identified in whole blood from people of African American ancestry, while this cohort is predominantly people of European ancestry, meaning that eQTMs specific to people of European ancestry could have been missed. Additionally, the MESA eQTMs were identified in isolated monocytes not whole blood. For this reason, eQTMs were tested in this dataset separately from the eQTMs in the GTP dataset because, first, any eQTMs specific to other blood cell types would have been missed by this analysis, and second, the CpGs associated with PBB may not be associated in the monocyte proportion of whole blood which could mean that some of the eQTM overlap may be a false positive. Furthermore, while it is difficult to discern whether PBB and immune cell type composition have an association outside the association between age and cell type composition, finding an association between PBB levels and estimated B cell proportion is consistent with previous findings in the literature which found an increased risk for non-Hodgkin’s lymphomas in people exposed to PBBs [28] and altered proportions of B- and T-lymphocytes [69].

Finding association with both estrogen and immune function highlights the role that estrogen has on the immune system. Estrogen regulates differentiation of immune cells, is required for appropriate immune response to pathogens, and is postulated as one of the reasons that women have a stronger immune response and higher incidence of autoimmune disease (compared to men) [70–72]. Other estrogenic EDCs, like bisphenol A (BPA), are also associated with inflammation and increased incidence of type 1 diabetes in mice models [73–75]. Therefore, PBB may be affecting DNA methylation at certain sites in the epigenome like an estrogenic compound, and thus be altering gene expression and immune function more broadly.

This study presents the first evidence for PBB affecting the epigenome. We also found evidence that that PBB could be acting similar to an estrogen, could be impacting gene expression, and could be potentially dysregulating the immune system. However, the findings of this study should be interpreted in light of this study’s limitations. Because the people in Michigan were uniquely exposed due to a factory accident, there is no other cohort in which to replicate our results. Additionally, only blood samples were available for DNA methylation analysis, so it is not known if or how PBB would affect other tissues. This study also includes several related participants. However, because PBB exposure is highly collinear with family structure, it could not be accounted for in statistical models without eliminating all the variation in PBB levels. This could inflate our results, as could other unmeasured confounders, like smoking status and BMI, which were also not taken into account in our models since not all participants reported current health information. Additionally, we were limited in our ability to find enriched transcription factor binding sites near PBB-associated CpGs due to the limited availability of databases of transcription factor binding sites in human tissues. As it was, many pathways only had a few of the transcription factors involved in that pathway in the database, meaning that we could not test for every member of a pathway. Also, many of the transcription factor motifs for the members that were tested were based on studies of non-human vertebrates, meaning that we could be missing human-specific binding motifs. It is possible that as more transcription factor binding sites are identified and added to validated databases, more transcription factor binding sites would be enriched in our results. Furthermore, CRP, gene expression, and estradiol levels were not directly measured in this cohort. This meant that associations with inflammation, eQTMs, and estrogen required the PBB results to be subset to just the CpGs on the previous array types. This makes us unable to determine whether the enhancers associated with PBB exposure would also be associated with enhancers that respond to estrogen or inflammation. Therefore, there may be key biological pathways that we are missing at this time. However, using the results from these previous studies, we were able to find associations between PBB and the epigenome, estrogen, and the immune system. Further studies into associations between PBB exposure and estrogen levels, other markers of immune function like cytokines, and the epigenome of other relevant tissues are warranted.

In conclusion, PBB is a persistent organic pollutant that contaminated the food supply in Michigan in the 1970s, and has been associated with many endocrine-related health problems among the exposed individuals. Because Michigan residents still have detectable PBB levels, the children of exposed mothers are still exposed during development, and PBB is similar to other common environmental pollutants like PCB and PBDE, the effects of PBB on the epigenome were investigated. We found that 1890 CpGs across the genome associated with higher exposure to PBB, and that these CpGs were depleted in promoter regions and heterochromatic regions, but were enriched for enhancer and insulator regions. PBB-associated CpGs were near transcription factor binding sites for members of the AHR pathway and estrogen-related transcription factors. These CpGs were also enriched for eQTMs associated with immune function and autoimmune disease, and correlated with both estrogen and inflammation signals. Taken together, our results suggest that exposure to PBB affects specific sites in the genome, and may be affecting immune- and estrogen-related function through gene expression more generally rather than through a specific pathway.

Methods

Participant selection

Participants were selected from the Michigan PBB Registry. This cohort was started by the Michigan Department of Community Health (MDCH) after the agricultural accident and is now managed by Emory University. The MDCH enrolled individuals believed to have had the highest exposure to PBB: families that lived on quarantined farms, people who ate food from quarantined farms, and chemical workers and their family members. Biological samples and health information from the original registry participants, their children, and other members of the community who were exposed to PBBs are still being collected and added to the registry (http://pbbregistry.emory.edu/). Informed consent was obtained from each individual before participation. Study protocols were approved by the Institutional Review Board at Emory University. Participants were selected from the registry for this study if they were exposed to PBB prior to being 50 years old, had given information on their health outcomes, had a recent (2004-present) buffy coat or whole blood sample available for DNA extraction, had lipid levels measured, and had detectable PBB levels in their serum. A total of 666 participants matched these criteria and were selected. All samples were taken between 2004–2015 (31–42 years after the exposure incident) and PBB levels, lipid levels, and DNA methylation levels were all measured on the same sample for each participant. Current and past health and demographic information on all of these participants is still being collected.

Exposure assessment

209 possible congeners of PBB exist based on the number and position of the bromine molecules around the biphenyl rings [19]. The primary congener in the technical mixture that contaminated the food supply in the 1970s was PBB-153 [17–19]. Exposure to four congeners of PBB (PBB-153, PBB-101, PBB-77, and PBB-180) was previously assessed in members of this registry using gas chromatography-tandem mass spectrometry [20]. The limit of detection (LOD) for PBB-153 was 2 pg/mL (0.002 ppb); for PBB-77, it was 4.5 pg/mL (0.0045 ppb); for PBB-101, it was 3.9 pg/mL (0.0039 ppb); for PBB-180, it was 5.6 pg/mL (0.0056 ppb). The extraction recovery ranged from 83.2–99.2%. The accuracy ranged from 89–119% and the precision ranged from 2.8–8.5%. For the purposes of this study, the value for congeners below the LOD in a sample was imputed as the LOD divided by the square root of 2 [76]. The four congeners were variably correlated with each other (R = 0.23–0.99; Figure S1). The congeners were then summed to give a total PBB value per person. This was then transformed using a natural log so that the distribution was normal (Figure S2).

Lipid measurement

Total lipid level was calculated in this cohort as part of previous studies [25]. A Triglyceride Quantification Assay Kit (Abnova Corporation) was used to measure the total triglyceride content in serum, and a Cholesterol Assay Kit (Caymen Chemical Company) was used to measure total cholesterol content in serum. Both were done in accordance with the manufacturer’s recommendations. Total lipid amount was calculated based on these components as described elsewhere [77].

DNA extraction

Peripheral blood samples were collected from participants as part of the ongoing Michigan PBB Registry activities between 2004–2015. DNA was extracted from buffy coat samples using the QIAamp DNA Blood Mini Kit (QIAGEN, Hilden, Germany). DNA was extracted from the blood sample from which the exposure levels and lipid levels were assessed on in order to limit confounding.

MethylationEpic array

Methylation levels were measured across the genome at 866,895 sites using the Infinium MethylationEPIC BeadChip (Illumina, San Diego, CA) [78]. Briefly, 1µg of DNA from participants’ buffy coats was bisulfite converted, amplified, fragmented, and hybridized to the BeadChip according to the manufacturer’s instructions. DNA from a female, lymphoblast line was used as a technical replicate. Samples with low signal (detection p-value greater than 0.01) or missing data (missing in greater than 10% of probes), and probes with missing data (missing in greater than 10% of samples) were removed with CpGassoc [79]. Two samples were removed for low signal, and 5,982 probes were removed for missing data. Previously identified cross-reactive probes (N = 44,210) were also removed from the dataset, but 9,993 probes with a common SNP (MAF > 5%) in the CpG were not removed [80]. BMIQ was used to adjust for probe type and ComBAT was used to adjust for batch effects [81,82]. 6 samples that were mismatches for sex, genetically identical to other samples, or flagged for sample contamination were removed. The 57 SNP probes were removed from the dataset. This resulted in a final dataset of 658 participants and 816,999 probes. For each individual sample at each probe, the methylation proportion (β) at that site was calculated from the methylated (M) and unmethylated (U) signal as β = M/(U +M). Cell type estimation was calculated for each sample using Houseman’s method for estimating cell type proportion with methylation signals from CpGs that are distinct in blood cell types [83,84]. The DNA methylation data can be accessed on NCBI’s Gene Expression Omnibus (GSE116339).

Epigenome-wide association testing

A linear regression model was used to interrogate the association between the combined average methylation level across all 816,999 CpGs that passed QC (global methylation) and total PBB level. Age at sample collection, sex, lipid level, and cell type proportions estimated earlier were used as covariates. An alpha level of 0.05 was used to determine statistical significance. Linear regression models in CpGassoc then were used to interrogate the association between total PBB level and methylation proportion at every probe individually that passed QC (816,999 sites) [79]. Age at sample collection, sex, lipid level, and cell type proportions estimated earlier were used as covariates. A Benjamini-Hochberg False Discovery Rate (FDR) of less than 0.05 was used to adjust for multiple testing [85]. MissMethyl was used to determine whether any of the 329 KEGG pathways were enriched in the results in order to correct for the number of probes in each gene [86].

Enrichment tests for specific regions

Chi-square tests were used to determine whether any PBB-associated CpG sites were enriched based on either their position within a gene or their proximity to a CpG island. CpG positions were determined from the annotation provided for each probe by the manufacturer. Functional regions were determined by the ChromHMM tracks publically available on the UCSC Genome Browser [87–89]. Briefly, these tracks used data from multiple chromatin immunoprecipitation sequencing (ChIP-seq) experiments to annotate the genome into 25 categories based on the function of that region in that cell type. These 25 categories were then further condensed into 10 region types based on that category’s position relative to a gene, its function, and its effect on gene transcription. For the enrichment tests, the tracks from the most similar cell type to whole blood (the GM12878 lymphoblastoid cell line) were used. An alpha level of 0.05 was used to determine statistical significance for the positional and functional enrichment. Enrichment of transcription factor binding sites was determined by oPOSSUM-3, which is based on data from ChIP-seq data for various transcription factors. This reports a Z score for enrichment of transcription factor binding sites within 10 base pairs on either side of each CpG that was associated with PBB compared to 10 GC-matched CpGs that were not associated with PBB [90]. A Benjamini-Hochberg False Discovery Rate (FDR) of less than 0.05 was used to adjust for multiple testing.

Enrichment tests for eQTMS

Chi-square tests were used to determine whether any significant results were enriched in expression quantitative trait methylation (eQTM) regions. For this, eQTMs previously identified in blood were used. These eQTMs were identified through analysis of methylation and expression data from the Grady Trauma Project (GTP) and the Multi-Ethnic Study of Atherosclerosis (MESA) [49]. Because these were identified using data from the Infinium HumanMethylation450 BeadChip, the results from this study were first subset to the probes interrogated by both arrays (N = 422,298). An alpha level of 0.05 was used to determine statistical significance. Enrichment of the genes these CpGs associated with the expression of was conducted with DAVID [50,51].

Correlations with estradiol and immune signaling

Because PBB is a known endocrine-disrupting compound and the participants of the Michigan PBB Registry have reported endocrine-related health problems, whether the epigenetic signal from PBB overlapped with the epigenetic signal from estradiol was tested. Enrichment of estradiol associated CpGs was tested using a chi-square test and results from a previously published epigenome-wide study [48]. Because these were identified using data from the Infinium HumanMethylation450 BeadChip, the results from this study were first subset to the probes interrogated by both studies (N = 401,387). Additionally, the effect sizes between the CpGs that overlap in both datasets and were significant (FDR < 0.05) in the analysis with estradiol (N = 25,143) were compared using a Pearson’s correlation coefficient.

Because DNA methylation from immune cells was tested, the correlation between the epigenetic signal from PBB was correlated to epigenetic-based immune phenotypes. To test for correlation with immune phenotypes, results from a previously published study of low-grade inflammation (C-Reactive Protein, or CRP) were used. Similar to the analysis with estradiol, the effect sizes between the CpGs that overlap in both datasets and were significant (FDR < 0.05) in the analysis with CRP (N = 207) were compared using a Pearson’s correlation coefficient [52]. No CpGs were found significant in both analyzes, and so enrichment tests could not be done. The estimated cell type proportions, as well as the neutrophil to lymphocyte ratio and CD4T to CD8T ratio (calculated from the cell types estimated from DNA methylation) and regressed on total PBB, both univariate and controlling for age and gender [53]. An alpha level of 0.05 was used to determine statistical significance.

Acknowledgments

We are grateful to the members of the Michigan PBB Registry for their participation and engagement with research studies over the past 40 years, to the Michigan Department of Health and Human Services, which had the foresight to create the Registry, and to our community partners (PBB Citizens Advisory Board, Pine River Superfund Citizen Advisory Group, and the Mid-Michigan District Health Department) who continue to provide guidance and insight to the Michigan PBB Research.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here

Additional information

Funding

This work was supported by the National Institute of Environmental Health Sciences (NIEHS; 5R01ES024790, 5R01ES025775; 5P30ES019776) and the National Institute of General Medical Sciences (T32GM008490). This study was supported in part by the Emory Integrated Genomics Core (EIGC), which is subsidized by the Emory University School of Medicine and is one of the Emory Integrated Core Facilities. Additional support was provided by the Georgia Clinical & Translational Science Alliance of the National Institutes of Health under Award Number UL1TR002378. The content is solely the responsibility of the authors and does not necessarily reflect the official views of the National Institutes of Health.