Exposure to polybrominated biphenyl and stochastic epigenetic mutations: application of a novel epigenetic approach to environmental exposure in the Michigan polybrominated biphenyl registry

ABSTRACT

Endocrine-disrupting compounds are associated with altered epigenetic regulation and adverse health outcomes, although inconsistent results suggest that people have varied responses to the same exposure. Interpersonal variation in response to environmental exposures is not identified using standard, population-based methods. However, methods that capture an individual’s response, such as analyzing stochastic epigenetic mutations (SEMs), may capture currently missed effects of environmental exposure. To test whether polybrominated biphenyl (PBB) was associated with SEMs, DNA methylation was measured using Illumina’s MethylationEPIC array in PBB-exposed individuals, and SEMs were identified. Association was tested using a linear regression with robust sandwich variance estimators, controlling for age, sex, lipids, and cell types. The number of SEMs was variable (range: 119–18,309), and positively associated with age (p = 1.23e-17), but not with sex (p = 0.97). PBBs and SEMs were only positively associated in people who were older when they were exposed (p = 0.02 vs. p = 0.91). Many subjects had SEMs enriched in biological pathways, particularly in pathways involved with xenobiotic metabolism and endocrine function. Higher number of SEMs was also associated with higher age acceleration (intrinsic: p = 1.70e-3; extrinsic: p = 3.59e-11), indicating that SEMs may be associated with age-related health problems. Finding an association between environmental contaminants and higher SEMs may provide insight into individual differences in response to environmental contaminants, as well as into the biological mechanism behind SEM formation. Furthermore, these results suggest that people may be particularly vulnerable to epigenetic dysregulation from environmental exposures as they age.

KEYWORDS:

• Endocrine-disrupting compound
• epigenetics
• DNA methylation
• polybrominated biphenyl
• PBB
• environmental health
• epimutations
• SEM
• interpersonal variation
• brominated flame retardant

Introduction

Endocrine-disrupting compounds (EDCs) are a class of chemicals that interfere with the endocrine system by disrupting the synthesis, binding, or transport of hormones [1]. They are common in the industrialized world and are present in pesticides, personal care products, electronics, fabrics, and household dust [2,3]. In human populations, increased exposure to EDCs has been linked to reproductive, hormonal, and developmental health problems [4–8]. However, studies of the health effects of EDC exposure have reported inconsistent results, with some studies finding an association between the same EDC and the same health outcome, and others not finding an association [9–11]. While some of this inconsistency is due to differences in study design, population, specific exposures, and dose, it is also known that there are differences in how individuals respond to and metabolize xenobiotics [12]. Therefore, it is reasonable to expect interindividual differences in the health consequences even among people with similar chemical exposures.

Exposure to EDCs has also been linked to epigenetic variation, both in cell lines and human cohorts, and these variations in epigenetic features have been suggested to be potential mediators between exposure to EDCs and the development of adverse health outcomes [13–19]. However, previous studies on the effect of EDCs on DNA methylation have relied on differences in either mean global DNA methylation differences (as measured by methylation of Alu, LINE, or other repetitive elements) or in the mean methylation at a particular site in the genome. For an association to be detected in these types of studies, the differences in methylation levels have to be shared by a large proportion of the cohort. While this approach is useful and effective and has found associations between exposure to EDCs, particular biological pathways, and health outcomes, it does not fully explain why people are susceptible to different health problems from the same environmental exposure.

It is possible that, in addition to affecting epigenetic marks in a way shared by a population, EDCs can also act randomly throughout the genome, altering different epigenetic marks in susceptible people and promoting distinct health outcomes. If this were so, these epigenetic differences would be missed by traditional analysis methods. Evaluation of rare stochastic epigenetic mutations (SEMs) would provide a measure of random differences in an individual. SEMs are defined as an individual having an extreme methylation value at a site when compared to the rest of the population. SEMs, by nature, are not shared among a population and do not have a substantial influence on the mean methylation level. However, despite their stochasticity, SEMs have been linked to cancer progression, ageing, and age-dependent X-inactivation skewing [20,21]. Therefore, they may be an important mediator between environmental exposures and health outcomes that could help explain the interpersonal variation in response to EDC exposure.

To investigate the association between exposure to environmental EDCs and SEMs, we utilized data collected as part of the Michigan PBB Registry. The Michigan PBB Registry includes approximately 7,000 individuals who were exposed to an EDC when an agricultural accident in the 1970s introduced polybrominated biphenyl (PBB) into the food supply. Individuals in this cohort were exposed to high levels of PBB during the nearly year-long period when contaminated food was distributed. Forty years later, a majority of participants still have higher PBB blood levels than 95% of the general United States [22]. Previous studies in this cohort have linked PBB exposure to earlier age of menarche, genitourinary conditions in males, thyroid dysfunction, and increased risk for certain types of cancer, particularly among those were exposed to PBB at a younger age [23–27]. Additionally, previous epigenetic research has found that higher PBB exposure is associated with methylation proportion at 1890 CpGs across the genome, and the effect that PBB has on the epigenome is suggestive of estrogen signalling and immune dysfunction [28]. However, it is not clear whether environmental stimuli, like PBB, would act randomly on the epigenome and thus associate with the number of SEMs. It is also not known if SEMs would depend on the age when they were first exposed to PBB or if they are sex-specific (like other previous research on the health effects of PBB), and it is not known if SEM count would associate with other age-related health outcomes like epigenetic age acceleration measures.

Results

Study population demographics

The study population was composed of 658 people, many of who were highly exposed to PBB (range: 0.01–236.73 ppb; Table 1). There were more female participants than male participants in this population (381 v. 277; p = 6.87e-05). A majority of the cohort (97%) was of White/Non-Hispanic ancestry, which is representative of the population of rural Michigan in the 1970s. Higher levels of PBB were associated with both being male (ln(PBB) in men vs women: −0.99 vs. −0.36, p = 4.54e-07) and older age (r = 0.25, p = 3.96e-11). PBB level and total lipid levels were not correlated in this cohort (r = −0.07, p = 0.06). Additionally, age of exposure was highly correlated with current age (r = 0.98; p < 2.2e-16), because most of the study population was exposed during the contamination event. Half of this study population was exposed before age 13.

Table 1. Characteristics of the members of the Michigan PBB Registry with epigenetic data.

	Study population demographics
	Total Cohort (N = 658)	Exposed when younger^d (N = 329)	Exposed when older^d (N = 329)	P-value
Number male^a	277 (42.09%)	97 (29.48%)	180 (54.71%)	2.90e-11
Current age (years)^b	54.28 ± 12.74	43.88 ± 6.12	64.69 ± 8.41	< 2.2e-16
Age when exposed (years) ^b	15.18 ± 11.55	5.71 ± 4.02	24.65 ± 8.43	< 2.2e-16
Total PBB level (ppb)^c	0.48 (4.71)	0.33 (4.81)	0.68 (4.24)	3.21 e-9
Total PBB level (ng per g lipids) ^c	73.09 (4.95)	50.79 (5.01)	105.16 (4.48)	3.44e-9
Number of SEMs ^c	471.99 (1.87)	418.15 (1.74)	532.77 (1.94)	5.79e-7
Race/Ethnicity^a				0.64
White/Non-Hispanic	638 (96.96%)	318 (96.65%)	320 (97.26%)
White/Hispanic	20 (3.03%)	11 (3.34%)	9 (2.73%)

^aFrequency and percentage.

^bMean and standard deviation.

^cGeometric mean and geometric standard error.

^dCompared to the median age when exposed (13 years).

SEM calculation and characteristics

SEM values were highly variable in this cohort, with a range of 119 to 18,309, and were not normally distributed (Figure S1A). Therefore, the natural log of the number of SEMs was used for all regression analyses (Figure S1B). SEMs were present in 12.75% of probes (Figure S2). In the participants who had multiple samples, the SEM calls were highly correlated within the longitudinal pair, and not correlated with the other samples (r = 0.55–0.77; p < 2.2e-16; Figure 1), indicating that SEM calls can be replicated within individuals. The number of SEMs in these sample pairs also increased between sampling for a majority of participants (mean difference in SEMs between the later time point and the earlier time point: 22.14 ± 47.64).

Figure 1. SEMs can be replicated in longitudinal repeats. Because DNA methylation levels were measured in longitudinal repeats for 14 individuals, with the first sample taken in 2004, and the second sample taken between 2012 and 2015 (designated below by the Subject ID followed by 1 or 2 to indicate the first or second sample), SEMs were calculated in both of these samples to see if they could be replicated. The SEM calls were highly correlated within the longitudinal pairs (r = 0.55–0.77) and not correlated within samples that were non-longitudinal pairs (r = 0.00–0.10), indicating that SEMs are not just the result of technical error, and can be replicated. In a majority of individuals (64.28%), the number of SEMs increased between samples (mean: 22.1447.64).

Similar to other SEM studies [20], SEMs were positively associated with current age, both when unadjusted (z = 7.11; p = 1.15e-12) and when adjusted for sex and cell type composition (z = 8.54; p = 1.23e-17) (Figure S3). Both the mean and the variance in the number of SEMs were higher in people who were older when they were exposed to PBB (Figure S4). Higher proportions of CD8T cells, B cells, monocytes and NK cells were associated with an increased number of SEMs while a higher proportion of CD4T cells and granulocytes were associated with decreased number of SEMs. While age was also associated with many of these covariates, the association between these covariates and number of SEMs is independent of their association with age (Table S1).

Higher number SEMs associates with higher PBB level

Current PBB levels had a positive, but non-significant association with the number of SEMs, once age, sex, lipid level, and cell type proportion were adjusted for (z = 1.57; p = 0.11; Figure S5). This association was largely consistent when the number of SEMs and each congener of PBB were analyzed individually (p = 0.003–0.16; Figure S6). However, because SEMs were more common in older individuals, and many PBB-associated health problems have been specific to people exposed at younger ages [9], we also tested whether PBB and age of exposure interacted to predict SEM count. PBB exposure occurred primarily during a 10-month period in 1973, making it possible to estimate the age of exposure to PBBs. Further analysis showed that there was a significant interaction between age at exposure to PBB and current PBB level to predict the number of SEMs (z = 2.27; p = 0.02). Stratifying the analysis by median age of exposure showed that there was an association between PBB level and SEM count, but only in the participants who were older when they were exposed to PBB (t = 2.22, p = 0.02 vs. z = 0.10, p = 0.91; Figure 2).

Figure 2. Higher SEM count associates with higher exposure to PBB in those exposed when they were older. The interaction between age of exposure and current total PBB level is associated with the number of SEMs (z = 2.27; p = 0.02). The association between total PBB and SEMs was only significant in participants who were exposed when they were older (black triangles, solid line) (z = 2.22; p = 0.02). For participants exposed to PBB when they were younger, the association between PBB and SEM count was not significant (grey circles, dashed line) (t = 0.20; p = 0.91). Figure was plotted on a log₁₀ scale for ease of interpretation, but the statistical analysis was done as described in methods.

SEMs by sex

SEM count was not different between men and women, once age and cell type proportions were adjusted for (z = −0.03, p = 0.97, Figure S7). There was also no interaction between sex and PBB level to predict SEM count (z = 1.06, p = 0.28). When the population was stratified by sex, the association between PBB and number of SEMs was non-significant in both males (z = 1.23, p = 0.21) and females (z = 0.63, p = 0.52).

SEMs are enriched in biological pathways

To test whether SEMs were more likely to be in certain biological pathways, and thus identify possible biological mechanisms behind SEM formation, an enrichment test was conducted with biological pathways defined by KEGG on the probes with SEMs. Of the 331 pathways tested, the CpGs with the most SEMs were significantly enriched in 94 of them (FDR < 0.05; Table S2). Of these 94, 55.3% are directly related to signalling pathways, and many others are related to cell cycle regulation and autophagy.

Next, to test whether the SEMs were randomly distributed across the genome of each participant or if some participants had clusters of SEMs in particular biological pathways, enrichment tests were conducted with biological pathways defined by KEGG. Of the 658 participants, 129 participants had SEMs enriched in at least one KEGG pathway, while 529 participants had no enrichment of their SEMs in any pathway tested. Participants whose SEMs were enriched in any KEGG pathway had a higher number of SEMs (geometric mean: 783.67 vs. 423.16; p-value < 2.2e-16), and a higher level of PBB (geometric mean: 0.68 vs. 0.44 ppb; p-value = 0.03) compared to the participants who had no enrichment (Figure S8). Men and women were equally likely to have pathway enrichment (p = 0.22). Participants whose SEMs were enriched in at least one KEGG pathway were also older compared to those whose SEMs were not enriched in any KEGG pathway (mean: 57.68 vs. 53.45 years; p-value = 0.002).

Of the 331 KEGG pathways tested, 85 pathways were enriched in at least one participant (Table S3). A majority (51.7%) were only enriched in one person, emphasizing the interpersonal variation of SEMs. However, there were 10 pathways that were enriched in more than 30 of the 129 participants (Figure S9). All of the pathways that were enriched in more than 30 participants were involved in endocrine function or xenobiotic metabolism. However, only one of these pathways was more likely to be enriched in people with high PBB level (compared to the median) (Table 2). Many of these pathways were more likely to occur in male participants compared to female participants (Table 2), indicating that SEMs formation and susceptibility may be under sex-specific influences. Participants with a higher number of SEMs (compared to the median) were not more likely to have enrichment of the common biological pathways (p = 0.42–1; Table S4), suggesting that the association between higher PBB exposure or sex and enrichment in these KEGG pathways was not due to the association between higher PBB exposure or sex and higher SEM exposure. These common pathways were not enriched for the pathways found enriched in the probes with SEMs (OR: 0.30, 95% CI: 0.01–1.69, p = 0.19).

Table 2. People in the Michigan PBB Registry have SEMs enriched in biological pathways related to xenobiotic metabolism and hormone regulation.

KEGG Pathway	Description	Number of people with SEMs enriched	PBB OR (95% CI)	PBB P-value	Male OR (95% CI)	Male P-value
path:hsa00980	Metabolism of xenobiotics by cytochrome P450	69	1.33 (0.81–2.23)	0.30	1.79 (1.08–2.98)	0.02
path:hsa00040	Pentose and glucuronate interconversions	69	1.33 (0.81–2.23)	0.30	1.79 (1.08–2.98)	0.02
path:hsa00830	Retinol metabolism	69	1.33 (0.81–2.23)	0.30	1.67 (1.01–2.78)	0.04
path:hsa00053	Ascorbate and aldarate metabolism	68	1.29 (0.78–2.17)	0.37	1.73 (1.04–2.89)	0.02
path:hsa00982	Drug metabolism – cytochrome P450	68	1.29 (0.78–2.17)	0.37	1.73 (1.04–2.89)	0.02
path:hsa00860	Porphyrin and chlorophyll metabolism	67	1.26 (0.75–2.11)	0.43	1.68 (1.01–2.81)	0.04
path:hsa00140	Steroid hormone biosynthesis	67	1.26 (0.75–2.11)	0.43	1.68 (1.01–2.81)	0.04
path:hsa05204	Chemical carcinogenesis	66	1.30 (0.78–2.20)	0.36	1.63 (0.97–2.73)	0.07
path:hsa00983	Drug metabolism – other enzymes	63	1.27 (0.75–2.17)	0.42	1.58 (0.93–2.67)	0.10
path:hsa04080	Neuroactive ligand-receptor interaction	37	2.45 (1.21–5.29)	0.01	1.18 (0.59–2.30)	0.75

SEM count associates with other well-characterized biomarkers of age

To test for potential health outcomes associated with higher SEM count, we tested the association between SEM count and two measures of age acceleration. Age acceleration was chosen because it associates with many health conditions and lifestyle factors [29–42], and because it is known from previous research that SEM count associates with age and age-related conditions [20]. SEM count was positively associated with both intrinsic age acceleration (z = 3.13, p = 1.70e-3, Figure 3(a)) and extrinsic age acceleration (z = 6.61, p = 3.59e-11, Figure 3(b)).

Figure 3. Higher number of SEMs is associated with age acceleration. The number of SEMs is positively associated with two measures of age acceleration (intrinsic: z = 3.13, p = 1.70e-3, Part A; extrinsic: z = 6.61, p = 3.59e-11, Part B), even when adjusted for age, sex, and estimated cell type proportions. Figure was plotted on a log₁₀ scale for ease of interpretation, but the statistical analysis was done as described in methods.

Discussion

While the association between environmental contaminants, including PBB, and DNA methylation levels shared by a population has been demonstrated before [13–15,28,43–46], this is the first study to show that increased exposure to endocrine-disrupting compounds associates with an increased number of stochastic epigenetic mutations in individuals. This study utilized samples from the Michigan PBB Registry, a cohort of people who were highly exposed to PBB due to an agricultural accident and still have PBB levels well above the national average. Increased level of PBB associated with increased number of SEMs in people who were exposed when they were older, even when controlling for age, sex, lipid levels, and estimated cell type proportions. These results suggest that the number of random SEMs throughout a person’s genome may depend not only on age, as has been shown before, but also on exposure to environmental contaminants [20].

The finding that PBB and SEM count only associated in the people who were exposed when they were older was the opposite of what was expected, given that many of the health consequences associated with PBB exposure, like earlier age of menarche, increased risk of spontaneous abortions, and offspring with genitourinary conditions, were in participants exposed in utero or in early childhood [23,24,47]. This could suggest that the mechanisms of SEM formation and the mechanisms for endocrine-related health problems are distinct and that populations exposed to different developmental time-points may differ in their risk for each. Additionally, this could be a reason that there are fewer health problems statistically associated with PBB exposure in people exposed as adults since each person may have different epigenetic marks affected and thus different health outcomes [9].

Because age at exposure and current age are highly correlated, finding an association with SEM count in older people could also suggest that people are more vulnerable to SEM formation from environmental exposure as they age. If vulnerability increases with age, people who were exposed when they were younger may become more susceptible to health conditions from environmental exposures as they age and gain more SEMs. This would indicate that some of the consequences of environmental contaminants like PBB would not be seen for decades and would be difficult to identify with standard epidemiological approaches. Unfortunately, while this cohort is unique in that the age when they were exposed to PBB can be estimated because the age at exposure and their current age are so correlated, these two possibilities cannot be separated.

We also analyzed whether any of the participants’ SEMs were enriched in any biological pathway or whether they were all randomly distributed throughout the genome. While a majority of participants did not have any biological pathway enriched, 19.6% of participants had SEMs enriched in at least one biological pathway. Participants with higher numbers of SEMs and higher PBB were more likely to have their SEMs enriched in at least one biological pathway. While pathways related to immune function, xenobiotic metabolism, and estrogen signalling have been suggested by previous epigenetic studies in this population [28], and many of the pathways enriched in this analysis are related to xenobiotic metabolism and hormone signalling, most of these 85 pathways enriched in participants have not previously been implicated in epigenetic analyses in these participants. These participants may be more likely to have adverse health effects since the epigenetic dysregulation accumulated within a biological pathway. It is also possible that the extreme difference in epigenetic regulation in certain biological pathways is compensatory. This would indicate that these participants are upregulating certain biological pathways in response to PBB. However, given the diversity of pathways affected in the participants, both of these possibilities are difficult to test with standard statistical methods.

There were several people who had SEMs enriched in biological pathways related to xenobiotic metabolism and endocrine function, even though most of these pathways were not more likely to have high SEM formation. We found that people with higher levels of PBB compared to the median were only more likely to have SEMs in one biological pathway (neuroactive ligand–receptor interaction), and people with higher numbers of SEMs were not more likely to have SEMs in any of these biological pathways. Enrichment of these pathways in multiple participants of the Michigan PBB Registry (especially since they are the only pathways enriched in multiple participants) is of interest because cytochrome 450 was observed to be upregulated in rats after they were exposed to PBB [48–50], and altered steroid hormone levels and endocrine disruption was found to be associated with PBB exposure, both in animal models [51–53] and in epidemiological studies [23,24,54,55]. Even though there was no enrichment of these pathways in people with increased exposure compared to other people in this study, part of this could be because of the extremely high PBB exposure in this population (median exposure is 27 times higher than the national average), thus making it a comparison of high exposure to higher exposure [22]. Because most people in this study have high PBB exposure, it is still possible that these pathways are upregulated following exposure to PBB, and this upregulation makes them more vulnerable to epigenetic mutations. While it has not been tested whether SEMs are more likely to form in active regions of the genome, these pathways have been associated with PBB exposure previously. Enrichment of biological pathways among our results could support the presence of an underlying biological mechanism by which SEMs contribute to health problems.

While SEMs were enriched in some biological pathways in multiple people, it is important to note that a majority of SEMs are randomly dispersed in different genomic locations in different people. Not all participants had the same pathways enriched, with only 129 of the 658 participants having enrichment of SEMs in any biological pathway and 44 of the 85 enriched pathways only being present in one person. The diversity and number of pathways enriched may help explain the individual susceptibility to different diseases following exposure to PBB, where individuals with epimutations in one biological pathway would have different health outcomes than individuals with epimutations in other biological pathways.

Beyond associations with environmental exposure, this study investigated associations with SEM formation and with sex. We found that SEM count, the association of PBB with SEM count, and the odds for enrichment of any biological pathways did not differ by sex. However, the enrichment of SEM in those common pathways involved in xenobiotic metabolism and hormonal regulation was more likely to occur in men compared to women. This could indicate distinct mechanisms behind SEM formation for men and women. It could also suggest that men may be more vulnerable to epigenetic dysregulation following exposure to PBB. While sex-specific epigenetic effects have not been studied with PBB exposure, epigenetic studies of the structurally related polychlorinated biphenyl (PCB) in other populations have reported that men are more vulnerable to epigenetic effects from EDC exposure [56].

We found that the CpGs with SEMs called were more likely to be in biological pathways directly involved with signalling, cell cycle regulation, and autophagy. This could indicate areas of the genome that are active and involved in signal transduction and cell division are more likely to have an error in their epigenetic regulation that leads to an SEM. While this would be consistent with the finding of xenobiotic metabolism and hormone signalling pathways enriched in members of the Michigan PBB Registry (since these pathways are expected to be upregulated following exposure to PBB), again, this has not directly been tested. More research is needed to investigate whether these pathways replicate in other studies of SEMs as well as to test mechanisms behind SEM formation.

This study also found that an increased number of SEMs was associated with age acceleration (as measured by both intrinsic age acceleration and extrinsic age acceleration) [34]. This is the first report of the association between SEM count and more established measures of age acceleration, but it is consistent given that SEM count is associated with older age and age-related health problems like cancer [20,21]. This could also mean that the same mechanisms that lead to the ageing process also lead to SEM formation. Additionally, because age acceleration has been associated with a myriad of adverse health outcomes, such as cancer, cardiovascular disease, and all-cause mortality [36,39,57,58], this could also indicate that an increase in SEM count from environmental exposure to EDCs could increase the risk for health problems later in life. However, more research is needed to test whether this association between SEM count and age acceleration replicates in other cohorts, as well as to test whether increased SEM count leads to increased age acceleration, or if an increase in age acceleration leads to more epigenetic dysregulation.

The findings from this study should be interpreted in light of its limitations. First, DNA methylation was only measured in blood samples, and therefore, it is not known if environmental exposure would associate with increased SEM count in other tissues. Also, it is not known if population-level exposure to other endocrine-disrupting compounds would also associate with number of SEMs, because only PBB levels were tested. Finally, while other studies have found an association between number of SEMs and health conditions, such as X-inactivation skewing and cancer progression [20,21], there were not enough participants with SEMs enriched in the same biological pathways for this to be directly tested in this cohort. Therefore, it is not known whether SEMs would mediate health conditions associated with PBB exposure in this cohort.

However, this study does have several strengths. We were able to implement a very stringent quality control pipeline by both removing probes with potential SNPs and implementing a lower detection p-value so that we could be more confident in our SEM calls. Additionally, because exposure to PBB occurred as part of a unique environmental contamination event, the age at which people were first exposed to PBB can be inferred from the available data, which allowed for the interaction between age at exposure and PBB level to be tested. In addition, in this study population, we were able to show that technical replicates were highly correlated, that the longitudinal pairs were highly correlated with themselves and uncorrelated to other samples, and that SEM calls in multiple tissues from the same person had very low correlation – supporting that SEMs are biological, epigenetic mutations that can be replicated and are not just the result of technical noise or underlying genetic variation. We were also able to replicate other studies finding that SEM count had a strong, positive association with increasing age, further supporting that these SEMs are biological, epigenetic mutations.

In conclusion, PBB exposure level is associated with increased number of SEMs, but only in people exposed when they were older. To our knowledge, this is the first report of environmental exposures associating with the number of rare, random, epigenetic mutations. It is also the first report of these SEMs being enriched in many biological pathways, especially in male participants. In addition, we found that SEMs are more common in signalling pathways, which may give insight into how SEMs are formed and which regions are more susceptible to SEM formation. We also found that SEM count is associated with increased age acceleration, which has not been shown before in previous studies. These SEMs may be important for understanding the variety of adverse health outcomes that are associated with EDC exposure and may mediate downstream effects on an individual’s health risks. More studies are warranted to test whether the adverse health outcomes reported by this population associate with the epigenetic mutations in their genome, and whether other environmental exposures associate with the formation of SEMs in particular biological pathways.

Methods

Participant selection

Participants were selected from the Michigan PBB Registry. The Michigan PBB Registry was started by the Michigan Department of Community Health (MDCH) after the cause of the agricultural accident was identified as PBB. The MDCH recruited individuals that lived on farms quarantined because of PBB contamination, obtained food from quarantined farms, and chemical workers and their families. This registry has recently been transferred to Emory University and has continued to enroll participants, collect biological samples, and survey health outcomes (http://pbbregistry.emory.edu/). For the current analyses, participants were selected from the registry if they met the following criteria: 1) they were exposed to PBB prior to being 50 years old, 2) they had a recent (2004-present) buffy coat or whole blood sample available for DNA extraction, and 3) they had current PBB and lipid level measurements in their serum. Informed consent was obtained from each individual before participation and study protocols were approved by the Institutional Review Board at Emory University. A total of 666 participants met these criteria and were selected for this study, including 14 participants with blood samples from two time points.

Exposure assessment

Two hundred and nine possible congeners of PBB exist and are defined based on the number and position of the bromine molecules around the biphenyl rings [59]. In the technical mixture of PBBs that was added to the food supply in Michigan, the primary congener was PBB-153 [59–61]. 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 [62]. The limits of detection (LOD) were 2 pg/mL for PBB-153, 4.5 pg/mL for PBB-77, 3.9 pg/mL for PBB-101, and 5.6 pg/mL for PBB-180. The extraction recovery ranged from 83.2% to 99.2%. The accuracy ranged from 89% to 119% and the precision ranged from 2.8% to 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 [63]. The congeners were then summed to give a total PBB value per person. This was transformed using a natural log so that the distribution was less skewed. For analyses using single congeners, if the congener was detected in more than 90% of the samples, then the transformed exposure level was used in analyses. If it was detected in fewer than 90% of samples, the exposure was dichotomized for analyses (detected in sample or not).

Lipid measurement

An Abnova Triglyceride Quantification Assay Kit (Abnova Corporation) was used to measure the total triglyceride content in serum, and a Caymen Cholesterol Assay Kit (Caymen Chemical Company) was used to measure total cholesterol content in serum. Both were done according to manufacturer’s instructions. Total lipid amount was calculated based on these components as described elsewhere [64,65].

MethylationEPIC array

Peripheral blood samples were collected from participants as part of the ongoing Michigan PBB Registry between 2004 and 2015. Blood was spun at 3,000 rpm to separate the plasma from the buffy coat. Buffy coats were aliquoted and stored at −80°C. DNA was extracted from buffy coat samples using the QIAamp DNA Blood Mini Kit (QIAGEN, Hilden, Germany). DNA was extracted from the same samples the exposure levels were assessed on in order to limit confounding.

Methylation levels were measured from these DNA samples at >850,000 sites using the Infinium MethylationEPIC BeadChip (Illumina, San Diego, CA) [66]. Briefly, 1 µg of DNA from participants’ buffy coats was bisulfite converted, amplified, fragmented, and hybridized to the BeadChip array according to the manufacturer’s instructions. Samples were randomized into 8 plates and 88 chips by age, sex, and PBB level in order to limit confounding. DNA from a stable lymphoblast line was used as a technical replicate, which was included randomly on the plates and chips (one replicate per half plate) so that technical replicates occurred in multiple positions on the chips. Samples with low signal intensity (detection p-value >1e-10) or missing data in >10% of probes, and probes with <3 beads or missing data in >10% of samples were removed using the wateRmelon package [67,68]. Two samples and 21,103 probes were removed for missing data. Forty-four thousand two hundred and ten previously identified cross-reactive probes, 19,681 probes on the X or Y chromosomes, and 18,100 probes with a common SNP within 5 base pairs of the assessed CpG (MAF > 1% in European populations from the Thousand Genomes Project) were also removed from the dataset [69,70]. To remove potential SNPs present in this population but not marked in the Thousand Genomes Project, probes with gaps (5%) in their methylation proportions in samples were identified using Gaphunter and also dropped [71]. Dye-bias correction and normalization were performed using ssNOOB, and ComBAT was used to adjust for chip-to-chip batch effects [72,73]. Six samples that were mismatched for self-reported sex, genetically identical to other samples, or had a higher number of Mendelian errors were removed. This resulted in a final dataset of 672 samples (from 658 participants, with 14 participants having two samples) and 769,042 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 proportion was estimated for each sample using Houseman’s method for estimating cell type proportion with methylation signals from CpGs that are distinct in blood cell types [74]. This data is publicly available at NCBI Gene Expression Omnibus (GSE116339).

SEM calculation

Samples were considered to have an SEM at a locus if their methylation level at that locus was either greater than three times the interquartile range above the third quartile or less than three times the interquartile range below the first quartile (either < or $Q1-(3\times IQR)or>Q3+\left(3\times IQR\right))$$Q1-(3\times IQR)or>Q3+\left(3\times IQR\right))$. This definition of SEMs is consistent with the approach described and validated by Gentilini and colleagues [20,21,75]. To do this, the IQR for each of the 769,042 probes that passed QC was calculated. Then, a new data matrix was annotated to determine whether each sample’s methylation proportion was an extreme outlier for each probe. The count of all the SEMs per sample was then calculated from the new data matrix. SEMs were also identified and counted in 15 technical replicates from cell lines derived from a single CEPH sample to assess the level of technical noise that could be due to batch effects (Figure S10). Additionally, to investigate the potential effect that genetic variant could have on SEM calls, SEMs were calculated in a publicly available dataset (GSE78743) [76], that contained 148 samples (11 different tissues taken from 16 different people) (Figure S11).

DNA methylation age

Two measures of DNA methylation age were calculated. Intrinsic age was calculated as the linear combination of the beta values of 353 CpGs that have been previously shown to predict chronological age in multiple tissues by Horvath [34]. Extrinsic age was calculated similarly, but from 71 CpGs that have previously shown to predict chronological age in whole blood by Hannum, et. al [77]. Extrinsic age acceleration was designed for use in blood samples, and so may not be well suited for other tissues. Intrinsic age acceleration, because it was developed in multiple tissues, is independent of age-related changes in blood cell composition. For this study, intrinsic and extrinsic age were calculated from the background-corrected beta values using Horvath’s publically available online calculator (https://dnamage.genetics.ucla.edu/home). Both measures first calculate a DNA methylation-based age, and then calculate the age acceleration from the residuals of the regression of DNA methylation age on chronological age.

Statistical analysis

The similarity between SEM calls called in the longitudinal pairs, the technical replicates, and the multiple tissues from the same people was assessed using a Pearson’s correlation coefficient for the presence/absence of all loci considered as part of these analyses. Linear regression models were used to determine the association between PBB serum level (as the independent variable) and the total sum of SEMs for a participant (as the dependent variable) in an analysis of the most recent blood sample from the 658 participants. Age, sex, total lipid levels, and estimated cell type proportions were included as covariates [78]. Because the number of SEMs was more variable in participants who were older or had higher exposure to PBBs, robust sandwich variance estimators were calculated for each of the cross-sectional analyses. Interaction between PBB and age of exposure was evaluated by the inclusion of an interaction term (age at exposure × total PBB level) in the model with main effects for PBB level, current age, sex, total lipid level, and cell type estimates. The study population was then stratified into participants who were exposed when they were younger (age of exposure less than the median age of exposure) and participants who were exposed when they were older (age of exposure greater than the median age of exposure). The association between current PBB level and SEM count was then tested in each subset separately, with the same covariates as before. As a sensitivity analysis, the association between each of the four detected PBB congeners was also assessed using the same covariates. The association between SEM count and sex was assessed with the same model as above, but with just age and cell types as a covariate. Interaction between PBB and sex was evaluated by the inclusion of an interaction term (sex × total PBB level) to the model with PBB level, current age, sex, total lipid level, and cell type estimates. The study population was then stratified into male participants and female participants, and the association between PBB and SEM count was assessed in each sex separately, with age, lipid levels, and cell type proportions as covariates. The association between each age acceleration measure was also assessed via the same model, with age, sex, and estimated cell type proportions as covariates.

Enrichment analysis

MissMethyl was used to assess whether an individual’s SEMs were enriched in particular biological pathways (defined by KEGG) or were randomly distributed throughout the genome [79]. This analysis was conducted in each sample separately. The false discovery rate (FDR) was controlled at 5% in each analysis. For pathways identified in more than 30 participants, a chi-square test was used to determine if the participants with the highest levels of PBB (based on a median split), highest number of SEMs (based on a median split) or who were male were more likely to have each of the common biological pathways enriched. Additionally, using MissMethyl, we assessed whether the probes with SEMs were enriched in any biological pathways.

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; National Institute of General Medical Sciences [T32GM008490].