DNA methylation profiling implicates exposure to PCBs in the pathogenesis of B-cell chronic lymphocytic leukemia.

OBJECTIVES: To characterize the impact of PCB exposure on DNA methylation in peripheral blood leucocytes and to evaluate the corresponding changes in relation to possible health effects, with a focus on B-cell lymphoma.
METHODS: We conducted an epigenome-wide association study on 611 adults free of diagnosed disease, living in Italy and Sweden, in whom we also measured plasma concentrations of 6 PCB congeners, DDE and hexachlorobenzene.
RESULTS: We identified 650 CpG sites whose methylation correlates strongly (FDR < 0.01) with plasma concentrations of at least one PCB congener. Stronger effects were observed in males and in Sweden. This epigenetic exposure profile shows extensive and highly statistically significant overlaps with published profiles associated with the risk of future B-cell chronic lymphocytic leukemia (CLL) as well as with clinical CLL (38 and 28 CpG sites, respectively). For all these sites, the methylation changes were in the same direction for increasing exposure and for higher disease risk or clinical disease status, suggesting an etiological link between exposure and CLL. Mediation analysis reinforced the suggestion of a causal link between exposure, changes in DNA methylation and disease. Disease connectivity analysis identified multiple additional diseases associated with differentially methylated genes, including melanoma for which an etiological link with PCB exposure is established, as well as developmental and neurological diseases for which there is corresponding epidemiological evidence. Differentially methylated genes include many homeobox genes, suggesting that PCBs target stem cells. Furthermore, numerous polycomb protein target genes were hypermethylated with increasing exposure, an effect known to constitute an early marker of carcinogenesis.
CONCLUSIONS: This study provides mechanistic evidence in support of a link between exposure to PCBs and the etiology of CLL and underlines the utility of omic profiling in the evaluation of the potential toxicity of environmental chemicals.

Introduction

Chlorinated persistent pollutants (POPs) are a category of environmental pollutants which are causing substantial health concerns (El-Shahawi et al., 2010; Faroon and Ruiz, 2015). They include polychlorinated biphenyls (PCBs), various organochlorine pesticides such as DDT (and its breakdown product DDE) and hexachlorobenzene (HCB), as well as numerous additional chemicals which were previously used for industrial or agricultural purposes. Although the use of these chemicals has ceased since many years, their resistance to degradation results in their wide persistence in the environment, including air, soil and water. Owing to their high lipophilicity, POPs accumulate along the food chain, with the consequence that humans are exposed to them primarily via the diet, especially the consumption of contaminated fish, meat and dairy products.

Significant experimental and epidemiological evidence suggests that exposure to chlorinated POPs may be linked to adverse effects on the immune, endocrine, nervous and reproductive systems, developmental effects and cancer (Crinnion, 2011; Everett et al., 2011; Lind et al., 2012; Perkins et al., 2016). In particular as regards cancer, a recent in-depth evaluation of the epidemiological and mechanistic evidence by the International Agency for Research on Cancer (IARC) concluded that the evidence linking exposure to PCBs with the induction of melanoma is sufficient to allow classification of this group of chemicals as category 1 human carcinogens (IARC, 2016).

The mechanisms by which chlorinated POPs cause their toxic effects are not well understood. Most have low genotoxicity, while many interact with important cellular receptors, including the Ah, estrogen and androgen receptors, and it is possible that such interactions may be important for these chemicals’ toxicity (Mrema et al., 2013). In order to explore the mechanistic basis of possible links between exposure to POPs and disease, a small number of studies have examined changes in genome-wide gene expression in peripheral blood leucocytes of exposed humans. Thus a study on pre-pubertal girls found changes in the expression of genes linked to connective tissue, skeletal muscular and genetic disorders as well as neurological diseases (Mitra et al., 2012), while a more recent follow-up study (Ghosh et al., 2018) on a mixed-sex group of similar age found gene expression changes linked to various types of cancer, including prostate and breast cancer as well as non-Hodgkin’s lymphoma. Recently we examined the association between exposure to a number of PCBs, HCB and DDE, a number of PCBs, HCB and DDE, and miRNA expression profiles in peripheral blood leucocytes of adults, identifying a series of expression changes related to various types of cancer, including lung, bladder, prostate and thyroid cancer, as well as chronic myeloid leukemia (Krauskopf et al., 2017).

Here we report the results of a genome-wide investigation of the associations between the concentrations of 6 PCBs, DDE and HCB in peripheral blood plasma of adult subjects without diagnosed disease and the methylation of CpG sites in peripheral blood leucocytes, which allowed us to characterize exposure-associated epigenetic profiles and to evaluate their significance in relation to the chemicals’ toxicity. In addition, and having in mind the contradictory epidemiological evidence regarding the relationship between PCB exposure and risk of B-cell lymphoma (IARC, 2016; Zani et al., 2017), we compared the exposure-related epigenetic profiles with the epigenetic profile in pre-diagnostic blood leukocytes we recently found to be associated with risk of future B-cell chronic lymphocytic leukemia (CLL) (Georgiadis et al., 2017) as well as with an epigenetic profile reported to characterize clinical CLL (Kulis et al., 2012).

Methods

Study population

The study was conducted in the context of the European EnviroGenomarkers project (http://www.envirogenomarkers.net/). It involved subjects, free of diagnosed disease at recruitment, from two population-based cohorts, the European Prospective Investigation into Cancer and Nutrition study (EPIC-ITALY) (Bingham and Riboli, 2004) and the Västerbotten Intervention Programme within the Northern Sweden Health and Disease Study (Hallmans et al., 2003) (Table 1). Standardized lifestyle and personal history questionnaires, anthropometric data and frozen blood fractions, collected at recruitment (1993–1998 for EPIC-ITALY, 1990–2006 for NSHDS), were available. The Envirogenomarkers project was originally designed as two nested case-control studies, one for B-cell lymphoma and one for breast cancer. No participant was diagnosed with disease within < 2 years of blood sample collection and for this reason in the context of the present study all participants were treated as apparently healthy at recruitment. Incident disease cases, including B-cell lymphoma, were identified through local Cancer Registries (loss to follow-up < 2%) and occurred between 2 and 15.7 years after recruitment. B-cell lymphoma cases were classified into subtypes according to the SEER ICD-0–3 morphology (Fritz et al., 2000). Cases with CLL had a mean age at diagnosis of 59.0 (43.6–75.5) years and were diagnosed 6.9 (2.0–15.5) years after recruitment. Cases with other BCL subtypes had a similar age and time-to-disease distribution, with a mean age at diagnosis of 58.5 (30.1–73.5) years and a time to diagnosis of 6.0 (2.0–15.9) years.

Table 1

Study population.

	Total population	EPIC Italy	NSHDS
All study subjects	659	251	408
Excluded from the current study
Missing data or extreme exposures	20	3	17
CLL cases	28	9	19
Included in the current study
All subjects	611	239	372
Age; mean (SD)	52.2 (7.8)	53.3 (8.1)	51.5 (7.5)
BMI; mean (SD)	25.8 (3.9)	25.8 (3.6)	25.9 (4.1)
Sex
Male (%)	215 (35.2)	59 (24.7)	156 (41.9)
Female (%)	396 (64.8)	180 (76.3)	216 (58.1)
Smoking status
Current smokers (%)	140 (22.9)	61 (25.5)	79 (21.2)
Never smokers (%)	287 (47.0)	111 (46.4)	176 (47.3)
Former smokers (%)	184 (30.1)	67 (28.0)	117 (31.5)
Health status
Controls (%)	316 (51.7)	123 (51.5)	193 (51.9)
Future cases (%)	295 (48.3)	116 (48.5)	179 (48.1)
Disease (future cases)
Breast cancer	91	46	45
BCL	204	70	134
BCL subtypes
DLBL	40	11	29
FL	32	19	13
MM	66	21	45
Other	66	19	47

The project and its associated studies and protocols were approved by the Regional Ethical Review Board of the Umea Division of Medical Research, as regards the Swedish cohort, and the Florence Health Unit Local Ethical Committee, as regards the Italian cohort, and all participants gave written informed consent. The studies were conducted in accordance with approved guidelines.

Analytical procedures and data processing

All analytical and data processing procedures employed, including DNA methylation and gene expression profiling, have been previously described in detail (Georgiadis et al., 2016). Genome-wide analysis of CpG methylation was conducted on the Illumina Infinium HumanMethylation450 platform and, after preprocessing, yielded data on 396,808 CpG sites. Methylation levels were expressed as M-values corresponding to the logarithmic ratio of the methylated versus the unmethylated signal intensities.

Plasma POP concentrations were measured as previously described (Kelly et al., 2017) by a procedure involving protein precipitation with ethanol, extraction of the POPs into dichloromethane–hexane and analysis gas chromatography–mass spectrometry. For quality control purposes in each batch of samples two reagent blanks were additionally prepared and the average result of the blank samples subtracted from the results of the real samples. Furthermore, two control samples of Standard Reference Material 1589a (PCBs, Pesticides, BDEs, Dioxins/Furans in Humans) from the National Institute of Standards and Technology, were also included in each batch (n = 43) of samples. Depending on the POP, mean concentration of SRM 1589a from all sample batches varied from 92% to 106% of certified values and co-efficient of variation from 3.8% to 10.7%.

Statistical analyses

Generalized linear models using the signals corrected for batch effects (date of chip analysis) were ran using the ArrayStudio (Omicsoft, Cary, NC, USA, version 8.0.1.32) software package, with inclusion of the moderated t-test (LIMMA) and filtration (with multiple testing accounted for using FDR Benjamini-Hochberg correction, alpha = 0.05 and maximum iterations = 5).

In the statistical models for the derivation of exposure-associated profiles we used M-values as dependent variables, the plasma concentrations of the different POPs as the independent variable, and sex, age, BMI, cohort, health status (control or future case) as well as the six cell type fractions [CD4, CD8, NK cells, monocytes, B-cells, granulocytes; estimated from the methylation data using a published algorithm (Houseman et al., 2012)] as confounder variables. In some analyses additional parameters were included in the model as confounders, as detailed in the text. Multiple testing was accounted for by using FDR Benjamini-Hochberg correction.

For the derivation of epigenetic profiles associated with future risk of different sub-types of B-cell lymphoma we compared the DNA methylation profiles of subjects who later developed B-cell lymphoma and control subjects who remained free of any disease (Table 1). In the statistical models we used future disease status as the independent variable and the same set of confounder variables as above unless otherwise stated.

Following exploratory evaluations, in the statistical modelling we adopted the plasma POP concentrations winsorised at 1% and 99% in order to control for a small number of subjects with outlier levels of particular POPs (see Supplemental Material, Section 1). We also explored the impact of winsorising the M-value distributions and came to the conclusion that this was not necessary. Venn diagrams were prepared using the software VennPainter (Lin et al., 2016).

Mediation analysis

Model-based causal mediation analysis was implemented using the R package “mediation” (Tingley et al., 2014). A customized R script was developed to iteratively construct the appropriate mediator and outcome models for each selected CG site. Each mediator model consisted of a linear regression fit including exposure (PCB156 plasma concentrations), the confounder variables (sex, age, BMI, white blood cell fraction) and using the methylation M-values of the corresponding CpG site as the dependent variable. Similarly each outcome model comprised a probit regression fit with both PCB156 concentrations and CpG methylation included as independent confounder variables and using the future case/control status as the dependent variable. During each iteration the two constructed models were used as input for the “mediate” R function, declaring PCB156 exposure as the treatment variable (“treat” argument), the CpG methylation as the mediator (“mediator” argument) and running 10,000 simulations (“sims” argument = 10,000). The final results were filtered using a p-value cutoff of 0.05 for the average causal mediation effects (ACME).

Bioinformatics analysis

Gene names obtained from the ArrayStudio output were checked with the on-line HGNC (HUGO gene nomenclature committee) tool (https://www.genenames.org/cgi-bin/symbol_ checker) and the returned names were subsequently used for bioinformatics analysis.

GO term analysis and identification of hub genes (genes linked to multiple GO terms and therefore playing a central regulatory role) were performed using the BioinfoMiner web application (https://bioinfominer.com/) which, thanks to its nonparametric, empirical prioritization approach, can be applied to classes of statistical testing problems that deflect from traditional hypotheses, as is the case for DNA methylation profiles. Pathway and disease connectivity analysis were performed using the “set analyzer” tool of the Comparative Toxicogenomic Database (http://ctd.mdibl.org), which utilizes manually curated information about chemical-gene/protein-disease relationships.

Results

POP exposure assessment

We measured plasma POP concentrations in 659 subjects aged 29.6–74.9 years from two prospective cohorts (Table 1). For the derivation of POP-related epigenetic exposure profiles we excluded 1 subject with outlier levels of all POP exposures and 19 subjects with missing relevant data. We also excluded 28 subjects who later developed CLL because we have previously observed (Georgiadis et al., 2017) that some of these subjects had major perturbations of their epigenetic profiles owing to large increases in their B-cell counts (no analogous effect was seen with other B-cell lymphoma subtypes). Of the remaining 611 subjects, 316 remained disease-free during the observation period (“controls”) while the remaining 295 were diagnosed within 2–15.7 years of recruitment with breast cancer or different subtypes of B-cell lymphoma other than CLL (“future cases”).

PCB exposures, as reflected in the plasma concentrations, were broadly similar for the two cohorts and the two sexes, although small but statistically significant differences were observed for some congeners (Table 2). In contrast, the mean exposures to HCB and DDE were substantially (roughly 3fold) higher in Italy than in Sweden. The exposures to the different PCB congeners were highly inter-correlated, with most Spearman r values > 0.8 (slightly lower for PCB118; Table S1 in Supplementary Text). The exposures to HCB and DDE were moderately correlated to each other and poorly correlated to those of PCBs.

Table 2

POP exposures by cohort and sex, mean ± SD (pg/ml).

	Italy	Sweden	p	Males	Females	p
PCB118	213.7 ± 134.0	145.4 ± 103.9	< 1 × 10−5	152.6 ± 116.1	182.8 ± 122.4	< 1 × 10−5
PCB138	571.2 ± 297.5	632.3 ± 389.8	ns	653.4 ± 438.9	584.0 ± 301.0	ns
PCB153	1112.3 ± 561.7	1116.9 ± 540.6	ns	1162.4 ± 580.6	1089.4 ± 527.7	ns
PCB156	95.6 ± 48.8	101 ± 50.2	ns	107.2 ± 56.2	94.3 ± 45.1	0.0099
PCB170	351.8 ± 187.1	385 ± 198.7	0.0035	414.7 ± 226.5	348.9 ± 170.4	0.00012
PCB180	846.7 ± 477.9	721.3 ± 309.7	0.015	810.6 ± 364.9	748.5 ± 399.0	0.0035
HCB	788.3 ± 634.8	246.1 ± 127.6	< 1 × 10−5	347.8 ± 396.6	518.2 ± 519.5	< 1 × 10−5
DDE	7485.6 ± 5947.1	2447 ± 2331.6	< 1 × 10−5	3551.7 ± 3964.2	4888.2 ± 5149.6	0.0002

Epigenetic exposure profiles

We used generalized linear models to evaluate the relationships between the methylation of different CpG sites and POP exposure levels. Since our aim was to evaluate the impact of POP exposure, as quantitatively reflected in blood-borne POP concentrations, on the epigenetic profile of white blood cells (i.e. a direct interaction between POPs and cells, both present in blood), we employed as measures of exposure the plasma POP concentrations (log-transformed) without correction for lipid concentrations. In addition we have confirmed that correcting for plasma lipid concentrations did not have a major impact on the resulting exposure profiles (see Supplementary Text, Section 2).

Table 3 summarises the numbers of CpG sites whose methylation correlates, at different statistical stringencies, with the exposure biomarkers. It can be seen that a) large numbers of statistically significant signals are observed in males, especially in Sweden, and b) most hits are associated with PCB156. Additional statistical adjustment for education and physical activity, consumption of alcohol and energy, as well as for exposure to DDE and HCB (both much higher in Italy), did not lead to convergence of the cohort- or sex-stratified results (not shown).

Table 3

Number of CpGs associated with exposure to different POPs, at different statistical stringencies.

Exposure	Statistical significance	Mixed cohorts			Italy			Sweden
		All	Males	Females	All	Males	Females	All	Males	Females
PCB118	Bonferroni p < 0.05	1	0	0	0	0	0	5	6	0
	FDR < 0.01	0	0	0	0	0	0	7	5	0
	FDR < 0.05	1	12	0	0	0	0	493	89	0
PCB138	Bonferroni p < 0.05	0	5	0	0	0	0	3	2	0
	FDR < 0.01	0	10	0	0	0	0	0	2	0
	FDR < 0.05	0	226	0	0	0	0	52	238	0
PCB153	Bonferroni p < 0.05	1	10	1	0	1	0	1	6	0
	FDR < 0.01	1	39	0	0	0	0	7	26	0
	FDR < 0.05	1	1303	1	0	56	0	220	1832	0
PCB156	Bonferroni p < 0.05	1	14	1	0	2	1	2	14	0
	FDR < 0.01	1	192	0	0	0	1	0	625	0
	FDR < 0.05	3	4606	1	2	33	2	6	7766	0
PCB170	Bonferroni p < 0.05	1	11	0	0	0	1	2	6	0
	FDR < 0.01	0	21	0	0	0	0	0	115	0
	FDR < 0.05	1	895	0	0	0	2	7	3117	0
PCB180	Bonferroni p < 0.05	2	5	0	0	0	0	0	4	0
	FDR < 0.01	0	6	0	0	0	0	0	29	0
	FDR < 0.05	2	301	0	3	0	0	4	2383	0
DDE	Bonferroni p < 0.05	0	0	0	0	0	0	4	7	0
	FDR < 0.01	0	0	0	0	0	0	3	10	0
	FDR < 0.05	0	0	0	0	0	0	267	213	0
HCB	Bonferroni p < 0.05	0	1	0	0	0	0	7	4	0
	FDR < 0.01	0	0	0	0	0	0	10	4	0
	FDR < 0.05	0	3	0	0	0	0	808	659	76

We carried out a series of additional tests to explore possible reasons for our failure to detect significant signals in the Italian cohort and in females, described in detail in Supplementary Text, Section 3. The results suggest qualitatively similar but substantially weaker responses in the Italian cohort and in females as compared to Swedish males, at least partly accounting for the near absence of statistically significant signals in these sub-groups.

Restriction of the analysis to the group of Swedish male controls, i.e. with exclusion of 72 subjects who eventually developed different subtypes of B-cell lymphoma, yielded 170 signals associated with PCB156 at FDR < 0.01 as compared to 625 signals obtained without this exclusion. As indicated in Fig. S2 in Supplementary Text, the two groups show qualitatively and quantitatively closely similar responses and the top signals in the two groups largely overlap, demonstrating absence of any bias in the profile resulting from the inclusion of case subjects.

Based on the above results, we conclude that the CpG methylation changes observed in the group of all Swedish males reflect qualitatively the effects of POPs on DNA methylation regardless of location, sex or future disease status, and for this reason the discussion which follows is based on the results obtained in this group, unless otherwise stated.

A total of 650 CpG sites are associated at high statistical stringency (FDR < 0.01) with exposure to at least one PCB (656 to at least one POP) (Table 3 and Excel Supplementary Table S1), with most being associated with PCB156 (625 sites) (Fig. 1A). The non-PCB POPs DDE and HCB yielded a much smaller number of significant signals, which largely overlap with PCB-associated signals (Fig. 1B). Based on data from the internal POP standards used in the study, the accuracy and precision in the measurement of the different congeners was similar and cannot explain the preferential association of signals with PCB156. Having also in mind the high inter-correlation of the exposure levels (especially of PCB’s), we conclude that the large number of signals which statistically correlate with specific congeners is unlikely to reflect true chemical-specific effects, rather probably arising from specific characteristics of the exposure distributions or chance. This possibility finds support in the observation in Table 3 of substantial numbers of signals associated with chemicals other than PCB156 when the statistical stringency is relaxed to FDR0.05 (see Discussion). For this reason further discussion is focused on signals associated with any PCB or POP.

Fig. 1.

Venn diagrams illustrating the overlaps between different PCBs (A) and PCBs and the two non-PCB POPs studies (B). Six hundred twenty five signals are associated with PCB156, of which 526 are associated exclusively with this exposure, followed by PCB170 (115, of which 16 are associated exclusively with this exposure).

Approximately equal numbers of CpG sites exhibit hypo- or hypermethylation with increasing exposure, with the mean change in methylation per quartile of PCB156 for the top signals ranging approximately 1–15% of the average methylation value.

Bioinformatics analysis of the POP exposure profile

The 656 differentially methylated CpG sites associated with at least one POP congener are related to 439 unique genes (including 20 hub genes; see Methods), shown in Excel Supplementary Table S2 together with various key characteristics. The list of differentially methylated genes includes a total of 15 homeobox genes (Zhong and Holland, 2011), all of which are hypermethylated with increasing exposure Bioinformatic analysis yields a large number of GO terms (Excel Supplementary Table S3) as well as 11 non-redundant pathway terms (Excel Supplementary Table S4).

Another notable feature of the list of differentially methylated genes is the presence of large numbers of polycomb group protein targets (PcGT’s), a category of genes whose promoter hypermethylation, and consequent expression downregulation, has emerged as a hallmark of the early stages of cancer pathogenesis (Widschwendter et al., 2018). Thus > 25% (121) of the differentially methylated genes belong to the class of PcGT genes (Bracken et al., 2006; Lee et al., 2006), the great majority of which are hypermethylated with increasing exposure at all their differentially methylated CpG sites (Excel Supplementary Table S2). Furthermore, the majority of 45 hypermethylated PcGT genes for which we had expression data showed a decrease in their expression which reached statistical significance for 5. Thus a picture emerges of POPs targeting for hypermethylation and downregulation homeobox and PcGT genes.

Disease connectivity analysis of the set of differentially methylated genes yielded a total of 64 significant non-generic terms (Excel Supplementary Table S5) which embrace, among disease categories, cancer (including melanoma) and diseases of the cardiovascular, nervous, urogenital, respiratory tract and immune systems as well as congenital abnormalities.

Comparison of POP exposure profiles with the profile predictive of CLL risk

We recently reported on an epigenetic profile in prediagnostic blood leucocytes which is strongly associated with future risk of CLL (Georgiadis et al., 2017). This profile includes 4295 significantly (FDR < 0.01) differentially methylated CpG sites and was derived from the comparison of the epigenetic profiles of 28 subjects, who were diagnosed with CLL 2–15.7 years after sample donation, with those of 319 subjects who remained free of disease, 315 of whom were included in the present study, coming from both cohorts and both sexes. Comparison of this profile with the POP exposure profiles described in Section 3.2 reveals overlaps of upto 38 CpG sites (p = 1.86 × 10−16), associated with 30 genes, a “meet-in-the-middle” (MITM) epigenetic profile which potentially represents a mechanistic link between exposure and disease (Tables 4 and 5). Importantly, for all MITM signals, the effects on methylation of a) increasing exposure and b) future CLL case status are in the same direction (Table 5), making the probability of a chance finding even more remote and strongly enhancing the biological significance of this overlap.

Table 4

Numbers of MITM CpG sites significant (FDR < 0.01) for both exposure to POPs and CLL risk using different sub-sets of subjects as well as different sets of statistical adjustments.

Model	CLL risk profile: subjects (number of signals)	Exposure profile	POP hits	Overlap (MUM)	p [a]	Comments
Main analyses (CLL risk profile in all subjects, exposure profilein Swedish males)
1	All subjects (4295) (Georgiadis et al., 2017)	PCB118	5	0
		PCB138	2	0
		PCB153	26	1	0.25
		PCB156	625	37	2.31 × 10−16
		PCB170	115	4	0.037
		PCB180	29	1	0.27
		Any PCB	650	37	7.93 × 10−16
		HCB	4	0
		DDE	10	0
		Any POP	656	38	1.86 × 10−16
Stability analyses
2	All subjects (4295) (Georgiadis et al., 2017)	PCB156 in all males	195	11	1.25 × 10−5	7 of the 11 MITM are also MITM in model 1; remaining 3 have FDR < 0.02 and 1 FDR < 0.05 for PCB156 in Swedish males
3	All male subjects (2893)	PCB156 in all males	195	7	6.32 × 10−4	4 of the 7 MITM is also MITM in model 1; 2 of the remaining have FDR < 0.02 for PCB 156 in Swedish males and for CLL risk in all subjects
4	All subjects, with additional adjustment for PCB156 (4161)	PCB 156 in Swedish males	625	36	5.20 × 10−16	35 of 36 MITM are also MITM in model 1; remaining signal has FDR < 0.05 in CLL risk profile without adjustment for PCB156
5		PCB 156 in Swedish males, adjusting for education and physical activity	496	15	3.07 × 10−4	12 of the 15 MITM are also MITM in model 1
6	All subjects (4295) (Georgiadis et al., 2017)	PCB156 in Swedish male controls	170	6	9.49 × 10−3	2 of 6 MITM are MITM in model 1; remaining 4 have FDR < 0.025 in 1
7	Swedish males (1434)	PCB 156 in Swedish males	625	8	2.24 × 10−3	6 of 8 MITM are in MITM of 1
8	Swedish males, with additional adjustment for PCB156 (1441)	PCB 156 in Swedish males	625	9	5.58 × 10−4	7 of 9 MITM are in MITM of 1

Total population N = 396,808.

Table 5

MITM CpG sites significant at FDR < 0.01 for both exposure to any POP and CLL risk.

CpG	Gene	Name	Associated exposure					Change in methylation with increasing exposure and CLL case status	Significant for CLL long time-to-disease	Gene targeted for extensive epigenetic modification in CLL	CLL risk hub	POP exposure hub	Homeobox gene	PcGT	Differentially methylated in clinical CLL (Kulis et al., 2012)
			PCB138	PCB153	PCB156	PCB170	PCB180
cg00352652	ZFPM1	Zinc finger protein, FOG family member 1			√			Down						√	√
cg00524900	TNFAIP8	TNF alpha induced protein 8			√			Down							√
cg00674365	ZNF471	Zinc finger protein 471			√			Up		√
cg00699993	GRIA2	Glutamate ionotropic receptor AMPA type subunit 2			√			Up		√				√	√
cg01100912	EFNA5	Ephrin A5			√			Up							√
cg01824511	FOXA1	Forkhead box A1			√			Up		√				√	√
cg02312409	RNF217-AS1	RNF217 antisense RNA 1 (head to head)			√			Up
cg03007522	GATA4	GATA binding protein 4			√	√	√	Up						√	√
cg03078269					√			Up							√
cg03646889	PLPPR4	Phospholipid phosphatase related 4			√			Up							√
cg03865667	PCDH17	Protocadherin 17			√			Up	√					√
cg04919489	ARHGEF12	Rho guanine nucleotide exchange factor 12			√			Down							√
cg08215169					√			Down
cg08543028					√			Down							√
cg09321400	SLC6A2	Solute carrier family 6 member 2			√			Up						√	√
cg10196720	PCDH10	Protocadherin 10			√			Up		√				√	√
cg10721834					√			up							√
cg11192895	LATS2	Large tumor suppressor kinase 2			√		√	Down
cg11428724	PAX7	Paired box 7			√			Up			√	√	√	√
cg14247287	NEURL3	Neuralized E3 ubiquitin protein ligase 3			√			Down							√
cg14849237	TLR5	Toll like receptor 5			√	√		Down						√
cg15912800	MIR196B	MicroRNA 196b			√			Up	√	√					√
cg17176573	POU2F3	POU class 2 homeobox 3			√			Up					√	√
cg18235050					√			Up							√
cg18256498					√			Down							√
cg19054524	PAX1	Paired box 1	√	√	√			Up				√		√
cg19384289	HOXD8	Homeobox D8			√			Up		√	√	√	√	√	√
cg19412467	ST6GAL2	ST6 beta-galactoside alpha-2,6-sialyltransferase 2			√			Up		√				√	√
cg19504702					√			Up							√
cg21229268	OLIG1	Oligodendrocyte transcription factor 1			√			Up		√				√	√
cg23111196						√		Down							√
cg23297413	ANKRD33B	Ankyrin repeat domain 33B			√			Down							√
cg23944804	BTBD3	BTB domain containing 3			√			Up						√	√
cg24843380	ZNF454	Zinc finger protein 454			√			Up		√					√
cg25026529	BARHL2	BarH like homeobox 2			√			Up	√	√	√		√	√	√
cg26987597	FOXF2	Forkhead box F2			√			Up	√					√
cg27062243	TCF7L2	Transcription factor 7 like 2			√			Down							√
cg27159979	BCL11A	B cell CLL/lymphoma 11A			√			Down						√	√

We carried out a series of additional tests to check the stability of the above MITM profile (Table4):

Comparison of the PCB156 exposure profile obtained in all males, rather than only Swedish males, with the CLL risk profiles obtained in all subjects or in all males, gave smaller but statistically highly significant MITM profiles which largely overlap with the one described above.

Use of the CLL risk profile obtained with additional adjustment for the level of exposure to PCB156 (to correct for any confounding by this or a correlated parameter) did not substantially change the resulting MITM profile, while adjustment of the exposure profile for education and physical activity yielded a smaller but significant and largely overlapping MITM.

Using the PCB156 exposure profile obtained in Swedish male controls (i.e. with the exclusion of all future cases of B-cell lymphoma) yielded a smaller but statistically significant MITM which largely overlaps with that observed without this exclusion. d) Finally, use of the CLL risk profile derived using only Swedish male subjects, without or with additional adjustment for PCB156, resulted in smaller but still statistically highly significant MITM overlaps.

Biological relevance of the MITM profile

Independent evidence in support of the relevance of the observed MITM profiles to the pathogenesis of CLL comes from the comparison its 38 CpG sites (MITM for exposure to any POP) with 33,653 sites whose methylation status has been reported to distinguish clinical CLL from normal B-cells (Kulis et al., 2012). This reveals an overlap of 28 sites (p = 1.98 × 10−22), for all of which the methylation changes in the same direction with increasing exposure and in clinical CLL (Table 5).

Additional features of the MITM profile shown in Table 5 include the presence of a) 4 CpG sites which we previously found to be significant in the risk profile of CLL cases who were diagnosed with the disease > 7.3 years after sample donation (Georgiadis et al., 2017), b) 10 MITM genes which are among 168 genes we previously reported to be targeted for extensive epigenetic modification in future CLL case subjects, and c) a number of genes which play hub gene roles in the CLL risk or/and the POP exposure profiles. Finally, the MITM profile includes 4 homeobox genes and 18 polycomb group protein target genes, with most of the latter being hypermethylated with increasing exposure at multiple CpG sites within the same CpG islands (coefficient > 0 and hypergeometric p < 0.05 in Excel Supplementary Table S6).

We conducted mediation analysis to evaluate the relationship between exposure to PCB156, future CLL case status and CpG methylation in Swedish males, using the 5 MITM CpG sites with highest statistical association (Bonferroni-corrected p < 0.05) with exposure to PCB156 or CLL risk. As shown in Table 6, significant mediation was found for 3 of these sites, although no statistically significant direct or total effect was observed. The absence of a significant total effect (direct association between POP exposure and CLL risk) is in agreement with our previously reported findings (Kelly et al., 2017) based on the full set of CLL cases of the Envirogenomarkers project (42 subjects), from which the subjects of the present study were drawn, as well as an analogous analysis based only on the cases included in the epigenetics dataset (see Supplementary Text, Section 5).

Table 6

Mediation analysis of the association between exposure to PCB156, DNA methylation and CLL risk.

MITM site	ACME (average causal mediation effects)		ADE (average direct effects)		Total effect
	Estimate	p	Estimate	p	Estimate	p
cg03865667	0.108	0.0052	−0.0943	0.37	0.0140	0.71
cg15912800	0.0015	0.0080	0.0010	0.67	0.00246	0.31
cg25026529	0.0462	0.012	−0.0176	0.90	0.0286	0.45
cg03007522	0.0088	0.085	8.04 × 10−3	0.98	1.68 × 10−2	0.68
cg00352652	0.0086	0.140	0.0120	0.42	0.0206	0.17

Other types of B-cell lymphoma

Comparison of the epigenetic profiles of future cases for the commonest lymphoma subtypes in our study with those of controls (Table 1) yielded risk profiles consisting of 1–3 CpG sites significant at FDR < 0.05 (Table 7), with no overlap between them or with the POP exposure profiles.

Table 7

Epigenetic risk profiles for different BCL subtypes.

Lymphoma subtypes	CpG	FDR	Raw p[a]	Coefficient	Gene	Gene name
Multiple myeloma	cg00036110	0.033	8.19 × 10−8	0.192	HPCAL4	Hippocalcin like 4
DLBL	cg10309377	0.038	9.68 × 10−8	−0.229
Follicular lymphoma	cg13267776	0.012	4.58 × 10−8	−0.434	CNNM2	Cyclin and CBS domain divalent metal cation transport mediator 2
	cg09851981	0.012	6.32 × 10−8	0.506	GOLGB1	Golgin B1
	cg06785701	0.012	8.92 × 10−8	−0.540	LOC407835	Mitogen-activated protein kinase kinase 2 pseudogene

Bonferroni-corrected p < 0.05 corresponds to raw p < 1.26 × 10−7.

Discussion

POP exposure-associated changes in blood leucocyte DNA methylation

In this, the largest epigenome-wide study to-date of the relationship between POP exposure and DNA methylation in peripheral blood leucocytes, we found that in males the methylation of large numbers of CpG sites is strongly associated with the plasma concentrations of at least one of 6 PCB congeners, DDE and HCB, the effect being strongest in Swedish males. While no statistically significant correlations were observed in a smaller group of Italian males or in females at either location, in these groups the response to exposure of the sites significant in Swedish males was qualitative highly similar to, but quantitatively 3–5fold smaller than that seen in the latter group, indicating differential sex- and location-related susceptibilities. A higher male susceptibility to POPs has been previously reported in relation to blood leucocyte LINE-1 DNA methylation (Lee et al., 2017), as well as in relation to a number of developmental effects (Hertz-Picciotto et al., 2005; Kishi et al., 2013; Sonneborn et al., 2008). Such sex-specific responses may result from the well-known interaction of POPs with key nuclear receptors, including the androgen and estrogen receptors (Bonefeld-Jørgensen et al., 2001; Zhang and Ho, 2011). The reason for the lower susceptibility of the Italian cohort is not known. The levels of exposure of the two cohorts to PCBs were generally similar (Table 2), while we have no evidence that the relative contribution of the routes of exposure for the general population (mainly ingestion) (IARC, 2016) differed substantially. We conclude that untested environmental or genetic factors may be responsible for the lower susceptibility of the Italian subjects.

The great majority of significant CpG sites were associated with exposure to PCBs, especially PCB156 (Fig. 1, Table 3). Given the strong inter-correlation of exposure to different PCB congeners (Table S1 in Supplementary Text), such apparently high chemical specificity is likely to be primarily related to the high statistical stringency employed and the exact exposure distribution or measurement error of the particular chemical, although the possibility that this particular PCB congener may possess a higher potency for altering DNA methylation cannot be excluded. PCB156 (2,3,4,5,3′,4′-hexachlorobiphenyl) is a mono-ortho PCB with significant but low dioxin-like activity (IARC, 2016). In a study conducted in Iceland Inuit with high POP exposures, PCB156 showed, among the PCBs examined by us, the highest association with the methylation of Alu repetitive DNA elements in blood cells (Rusiecki et al., 2008), although other studies also using global measures of DNA methylation gave mixed results (Itoh et al., 2014; Kim et al., 2010; Lind et al., 2013). In the only epigenome-wide evaluation of the effects of PCBs reported to-date (van den Dungen et al., 2017), conducted among 34 Danish males, no formally statistically significant associations of site-specific CpG methylation in blood leucocytes were found, while, of 8 differentially methylated regions identified, 4 included CpG sites whose methylation we found to correlate moderately (FDR = 0.025–0.075) with PCB exposure.

PCB-induced epigenetic changes in genes controlling the fate of stem cells

Among the CpG sites exhibiting strongest responses to PCB exposure (large absolute coefficient values; Excel Supplementary Table S1) are sites associated with many genes related to differentiation and development [e.g. ZFPM1 (zinc finger protein, FOG family member 1), erythroid and megakaryocytic cell differentiation; RDH10 (retinol dehydrogenase 10), organ development; TERT (telomerase reverse transcriptase), an antiapoptotic gene and modulator of Wnt signaling]. The importance of the modulation of the epigenetic status of developmental genes is particularly underlined by the large number of homeobox genes affected (15 of 439 differentially methylated genes) (Excel Supplementary Table S2). Homeobox genes act as master regulators in the renewal and fate of stem cells (Seifert et al., 2015), while their altered methylation is associated with cancer pathogenesis (Rodrigues et al., 2016). Therefore modulation of their epigenetic status by PCBs implies potential effects on development and carcinogenesis. Thus, among the differentially methylated hub homeobox genes are HHEX (hematopoietically expressed homeobox) and PAX6 (paired box 6), involved in hematopoietic (Migueles et al., 2017) and neural tissue differentiation (Huettl et al., 2016), respectively, WNT5A (Wnt family member 5A) which regulates pathways related to development, inflammation and cancer (Andersson et al., 2013; Endo et al., 2015; Pashirzad et al., 2017), HOXA9 (homeobox A9) and PBX1 (PBX homeobox 1), associated with myeloid leukemia/myelodysplastic syndrome and pre-B-cell acute lymphoblastic leukemia, respectively (Collins and Hess, 2016; Duque-Afonso et al., 2016), as well as RBP4 (retinol binding protein 4), RDH10 (retinol dehydrogenase 10) and ALDH1A2 (aldehyde dehydrogenase 1 family member A2), all involved in the biosynthesis of retinoic acid, an important signaling molecule in developing and adult tissues (Cañete et al., 2017).

The impact of exposure on stem cells is highlighted by the results of functional analysis (Excel Supplementary Tables S3 and S4) which yields multiple GO terms related to development, especially neurodevelopment, and perturbed pathways related to neurotrophins, a family of proteins which control the development and function of neuronal cells (Huang and Reichardt, 2001). Exposure to chlorinated POPs is well known to be associated with multiple effects on the nervous system, including neurological impairments (cognitive and peripheral nervous system effects, motor and sensory deficits) and neurodegenerative diseases [(Alzheimer’s and Parkinson’s disease, amyotrophic lateral sclerosis) in adults and neurodevelopmental diseases (autism, attention deficit, mental retardation, hearing loss) in children of exposed mothers (Zeliger, 2013)]. Recent evaluations of evidence from experimental and epidemiological studies support the suggestion that epigenetic changes induced by environmental exposures may mediate neurodevelopmental toxicity (Tran and Miyake, 2017).

Among the factors which determine the fate of stem cells are polycomb proteins, which transiently repress the expression of differentiation-promoting genes by binding to their promoters in the form of polycomb-repressive complexes (Mozgova and Hennig, 2015). During recent years strong evidence has accumulated indicating that, during the early stages of the pathogenesis of many types of cancer, including lymphomagenesis (Wang et al., 2015), the promoters of such prolycomb group protein target (PcGT) genes become methylated, and hence silenced, independently of the binding of such complexes, thus locking the cells in an undifferentiated state which predisposes them to malignant transformation (Martin-Perez et al., 2010; Widschwendter et al., 2007). Based on these observations it has been proposed that methylation of PcGT genes is an early hallmark of cancer (Teschendorff et al., 2010; Widschwendter et al., 2018). In the present study we found that a large fraction (121 out of 439) of genes differentially methylated in association with exposure to POPs belongs to the class of PcGTs (Excel Supplementary Table S2). The majority of these genes were hypermethylated with increasing exposure at multiple sites within CpG islands (Excel Supplementary Table S6), while 5 of them were significantly underexpressed, supporting the idea that POP exposure modifies cellular pathways involved in the early stages of carcinogenesis.

POP-induced epigenetic profile and disease

We have previously shown that omic profiles observed in peripheral blood leucocytes of healthy smokers predict with remarkable efficiency diseases caused by tobacco smoking (Georgiadis et al., 2016), suggesting that such profiling has the potential of identifying disease-related perturbations caused by toxic exposures. This potential is further supported by the results of disease connectivity analysis using our list of POP-related differentially methylated genes, which identified melanoma as being linked to this exposure (Excel Supplementary Table S5), in accordance with the conclusions of an IARC evaluation (IARC, 2016). Additional diseases suggested by our disease connectivity analysis include a number of diseases for which there is some supportive epidemiological evidence, including breast cancer (IARC, 2016) as well as diseases of the cardiovascular (Bergkvist et al., 2016; Kippler et al., 2016), digestive (Deierlein et al., 2017) and endocrine (Zong et al., 2018) systems. Furthermore, in agreement with the preceding discussion regarding effects on stem cells, numerous terms related to developmental and nervous system diseases and cancer are obtained.

Overlap of epigenetic profiles associated with PCB exposure and CLL risk

A recent in-depth evaluation concluded that, despite epidemiological and mechanistic data supporting a link between PCB exposure and risk of non-Hodgkin lymphoma (NHL), a definitive conclusion of positive association cannot be drawn (IARC, 2016). Other recent meta-analyses of the epidemiological data found no strong evidence that exposure to PCB increases the risk of NHL (Zani et al., 2017) and a significant positive association of exposure to DDE and HCB with risk of non-Hodgkin lymphoma (Luo et al., 2016).

We recently reported non-significant, positive associations between the plasma concentrations of most of the POP congeners examined in the present study and future risk of CLL and follicular lymphoma [see Kelly et al., 2017 as well as additional analyses in Supplementary Text, Section 5]. In striking similarity with the results of the present study, these associations were substantially stronger in males and in the Swedish cohort. In the current study we explored further the possible links between POP exposure and risk of B-cell lymphoma by comparing the epigenetic profiles associated with exposure to those associated with disease risk. The major finding of this exploration is the discovery of a large, statistically highly significant, overlap between the profile associated with the risk of future CLL and the profiles associated with exposure to any POP (38 MITM CpG sites; 37 sites for any PCB), with the direction of change of methylation in all cases being the same in subjects who ultimately developed CLL and in subjects with higher exposure (Table 5). The plausibility of this MITM profile is further enhanced by the fact that the CLL risk profile had been derived using an independent set of CLL case subjects, its stability to adjustment of the CLL risk profile for exposure to PCB156 and the observation of a smaller but statistically significant and partly overlapping MITM profiles using the exposure profile of control subjects alone (Table 4).

Biological plausibility of the MITM profile

Twenty eight of the 38 MITM CpG sites have been reported to be differentially methylated in clinical CLL relative to normal B-cells (Kulis et al., 2012), with the direction of methylation.

change in CLL being the same as observed in subjects with higher exposure for all 28 sites. This implies that the methylation changes induced at these sites by exposure occur early during disease pathogenesis or are present in clones of pre-clinical CLL-like cells, and are retained all the way to full clinical disease. It is noted that, in the study of Georgiadis et al. (2017) which identified the prediagnostic CLL risk profile employed in the present study, a progressive series of DNA methylation and gene expression changes in white blood cells of future CLL cases was identified, compatible with the presence in prediagnostic blood of CLL-like cells at different stages of progression towards clinical disease. That the DNA methylation changes associated with the MITM probably represent early perturbations on the disease pathogenesis pathway, rather than being present in latent CLL clones, is supported by the fact that 4 of the MITM CpGs (including 2 altered also in clinical CLL) are significant in CLL cases who were diagnosed with the disease > 7.3 years after sample donation (Georgiadis et al., 2017) (Table 5).

The biological plausibility of the MITM profile is further strengthened by the presence of 18 PcGT genes, most being differentially hypermethylated with increasing exposure, in line with the recognized significance of the hypermethylation of PcGT genes in carcinogenesis (Teschendorff et al., 2010; Widschwendter et al., 2018). Finally, a number of MITM genes have been implicated in the mechanism of carcinogenesis in B-cells. For example, BCL11a (B cell CLL/lymphoma 11A) is overexpressed in CLL, where it acts as an oncogene (Satterwhite et al., 2001) and protects CLL cells against apoptosis (Gao et al., 2013). LATS2 (large tumor suppressor kinase 2) is a tumor suppressor and has been found to be underexpressed in CLL (Ouillette et al., 2008). TLR5 (toll like receptor 5) plays a critical role in B-cell homeostasis and has been found to be mutated in CLL (Martínez-Trillos et al., 2014). Finally, MIR196B regulates a number of genes involved in B-cell differentiation and/or CLL, including the oncogene c-MYC (Pozzo et al., 2017), the anti-apoptotic gene BCL2 (Vogler et al., 2017) and the homeobox gene HOXA9 (Gwin et al., 2010).

It is also noted that 3 of 5 homeobox and PcGT genes (HOXA9, PAX6 and NOTCH4), which are hypermethylated and underexpressed at higher exposures, while not in the MITM profile, are known to be involved in lymphomagenesis (Collins and Hess, 2016). Finally, exposure is associated with the perturbation of multiple pathways related neurotrophin signaling which, in addition to its importance in determining the fate of neuronal cells, also plays an important role in carcinogenesis [including B-cell-related cancer (Hillis et al., 2016)], especially in relation to the control of cancer cell stemness.

Mediation analysis and possible causal links between POP exposure, DNA methylation and CLL risk

A statistically significant mediation effect between exposure to PCB156 and disease risk was found for 3 of the 5 MITM CpG sites most significantly associated with exposure or disease risk (Table 6). The involvement of these sites in the pathogenesis of CLL is biologically plausible since they are associated with PCDH17 [protocadherin 17, a tumor suppressor gene (Yin et al., 2016)], miR196B [hypermethylated in leukemia, thus allowing the upregulation of a number of oncogenes (Liu et al., 2013)] and BARHL2 [BarH like homeobox 2, hypermethylated in multiple cancer types (Rauch et al., 2012) and a regulator of proliferation and survival (Juraver-Geslin et al., 2011)]. Given this biological plausibility, the absence of a statistically significant total effect probably reflects study size limitations, in combination with a temporally distal relationship between exposure and disease (in our case 2–15.7 years) (Hayes, 2009), demonstrating the potential of epigenetics-based intermediate biomarkers in the investigation of exposure-disease risk associations.

Risk profiles of other subtypes of B-cell lymphoma

The number of epigenetic signals found to be associated with the risk of future MM, DLBL or FL is very much smaller than that associated with risk of CLL (Table 7). It is likely that this difference reflects, at least to some extent, the accumulation of large clones of pre-CLL cells in future CLL case subjects before they are diagnosed with this indolent disease. The few significant signals observed with the other lymphoma subtypes do not overlap with the epigenetic profile of POP exposure and therefore do not allow any evaluation of the possible association of this exposure with disease risk. Two studies of CpG methylation in clinical samples of follicular lymphoma, using early versions of microarrays, do not allow comparison with our lists (Killian et al., 2009; O’Riain et al., 2009), while a list of 794 CpG sites differentially hypermethylated in multiple myeloma (Agirre et al., 2015) does not include the CpG site we found to be associated with risk of this disease. Finally, there is no reported association of any of the genes of Table 7 with any subtype of B-cell lymphoma, although GOLGB1 has been reported to be involved in chromosomal mutations in hematologic neoplasias (Troadec et al., 2017).

Conclusions

The present study reveals an extensive and biologically plausible overlap between changes in DNA methylation induced by PCB exposure in subjects without diagnosed disease and corresponding changes in prediagnostic blood of subjects who later developed CLL as well as in clinical CLL. The preponderance in the epigenetic profile of PCB exposure of changes in homeobox and polycomb group target genes implies that stem cells may constitute critical targets of these pollutants in relation to their toxicity.

The main limitation of our study lies in our inability to directly replicate in the Italian cohort the effects of PCBs observed in Swedish males, probably owing to the small size of the corresponding population. Another shortcoming relates to the lack of information on the clinical state of the CLL cases at diagnosis, which limits our ability to characterize the CLL risk profile in relation to the possible presence of disease at the prediagnostic stage. However, despite these short-comings, overall our study adds to the weight of the evidence linking exposure to PCBs with the etiology of CLL. In addition our results underline the utility of blood-based profiling for the evaluation of the potential toxicity of environmental chemicals.