Perfluoroalkyl substances influence DNA methylation in school-age children highly exposed through drinking water contaminated from firefighting foam: a cohort study in Ronneby, Sweden.

Perfluoroalkyl substances (PFASs) are widespread synthetic substances with various adverse health effects. A potential mechanism of toxicity for PFASs is via epigenetic changes, such as DNA methylation. Previous studies have evaluated associations between PFAS exposure and DNA methylation among newborns and adults. However, no study has evaluated how PFASs influence DNA methylation among children of school age. In this exploratory study with school-age children exposed to PFASs through drinking water highly contaminated from firefighting foams, we aimed to investigate whether exposure to PFASs was associated with alteration in DNA methylation and epigenetic age acceleration. Sixty-three children aged 7-11 years from the Ronneby Biomarker Cohort (Sweden) were included. The children were either controls with only background exposure (n = 32; perfluorooctane sulfonic acid: median 2.8 and range 1-5 ng/ml) or those exposed to very high levels of PFASs (n = 31; perfluorooctane sulfonic acid: median 295 and range 190-464 ng/ml). These two groups were matched on sex, age, and body mass index. Genome-wide methylation of whole-blood DNA was analyzed using the Infinium MethylationEPIC BeadChip kit. Epigenetic age acceleration was derived from the DNA methylation data. Twelve differentially methylated positions and seven differentially methylated regions were found when comparing the high-exposure group to the control group. There were no differences in epigenetic age acceleration between these two groups (P = 0.66). We found that PFAS exposure was associated with DNA methylation at specific genomic positions and regions in children at school age, which may indicate a possible mechanism for linking PFAS exposure to health effects.

Introduction

The ubiquitous spread of per- and polyfluoroalkyl substances (PFASs), a group of persistent organic pollutants, in the environment has become a global contamination problem. Drinking water contamination of PFASs occurs not only from PFAS production facilities but also from military and civilian firefighting training facilities where large quantities of aqueous film–forming foams have been used [1-3]. Public concern about PFAS contamination has increased since PFAS exposure has been linked to adverse health effects in all age groups [4]. A growing body of research has indicated that PFASs are associated with several adverse effects in children, such as reduced immune response, dyslipidemia, later age of puberty onset, and decreased bone mineral density [5, 6].

For these different health outcomes, the modes of PFAS toxicity are still not fully known. One possible mechanism is through alteration in DNA methylation that influences gene expression, affecting phenotype and health. DNA methylation occurs predominantly on cytosine (C) followed by guanine (G) residues, referred to as CpG sites. DNA methylation has been proposed as a possible molecular mechanism underlying adverse health effects of various environmental pollutants [7]. DNA methylation has the additional advantage of being a sensitive intermediate effect biomarker since it can be detected prior to developing subclinical changes or pathologies [8].

DNA methylation has been studied at background PFAS exposure levels in newborns from mother–child cohorts [9-12] and adults [13, 14]. Studies at high-exposure levels are still few: one study in adults after high perfluorooctanoic acid (PFOA) exposure [15] and one in adult women from Ronneby with highly contaminated drinking water, dominated by perfluorooctane sulfonic acid (PFOS) and perfluorohexane sulfonic acid (PFHxS) [16]. No study has evaluated how PFASs influence epigenetics in children of school age. Given that age can influence the epigenetic response [17] and children are particularly susceptible to exposure to environmental toxicants [18], studies investigating PFAS-related DNA methylation in children are warranted.

Additionally, data from whole-genome DNA methylation arrays can be used to estimate the biological age based on a number of age-dependent CpG signatures [19, 20]. Accelerated aging, i.e. the difference between biological age and chronological age, has been associated with exposure to organochloride pesticides and polybrominated biphenyls [21, 22], as well as increased mortality [23] and impaired physical and cognitive functions [24] in adults. Although we did not find any association between accelerated aging and PFASs in our previous study on adult women [16], it is still not known if PFAS exposure is associated with epigenetic aging in children. Alterations in epigenetic aging have been reported in children after environmental exposure to other pollutants, such as indoor air pollutants and parental smoking [25].

We conducted an epigenome-wide analysis to examine the associations between PFAS exposure and DNA methylation in children at school age by contrasting a control group with only background exposure to a group with very high PFAS exposure. We also evaluated if PFAS exposure was associated with epigenetic age acceleration.

Materials and Methods

Study Population

The study participants were selected from the Ronneby Biomarker Cohort and its reference group. The Ronneby cohort and reference group details have been previously published [3, 26]. In short, the Ronneby Biomarker Cohort comprises 3297 individuals aged 0–92 years from Ronneby, Sweden, a municipality that had a drinking water supply highly contaminated by PFASs, dominated by PFOS and PFHxS, for decades until December 2013 in one of its two waterworks. The reference group comprised 226 individuals (age range 5–59 years) from Karlshamn, a nearby municipality with an uncontaminated municipal drinking water supply and similar geodemographic characteristics. The average serum PFAS levels in the Ronneby cohort were about 10 to 100 times higher than in the reference group for different PFAS compounds [3].

This study focused on school-age children from 7 to 11 years old. The selection was based on serum PFOS levels and intended to be two distinct groups with a large contrast in exposure: one exposure group with highly exposed children and one control group with only background exposure. Among the 3297 participants from the Ronneby Biomarker Cohort with high PFAS exposure, there were 175 children aged 7–11 years. Additionally, there were 69 children in the same age range in the reference group with background exposure. We first restricted the Ronneby children to those with PFOS levels above the median to ensure a large contrast on exposure levels to the control group. We then performed frequency matching on age, sex, and body mass index (BMI) to selected participants from Ronneby with above-median PFOS and participants from Karlshamn. The high-exposure group consisted of 32 children residing in Ronneby with a mean serum PFOS at 295 ng/ml (range 190–464 ng/ml), and the control group then consisted of 32 children from Karlshamn, with a mean serum PFOS at 3 ng/ml (range 1–5 ng/ml). Unfortunately, data on parental smoking status were not available, nor data on potential confounders such as socioeconomic status, pubertal status, and parental occupation. Prior to the analyses, one sample (from the high-exposure group) was removed due to a too low amount of DNA, resulting in 63 study participants.

The study was carried out in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki). The guardians of the participants gave informed written consent for biomarker studies. The Regional Ethical Review Board approved the study in Lund, Sweden (approved on 22 April 2014; approval number 2014/4).

Serum PFAS Concentration Analysis

Serum PFAS concentrations were analyzed as described by Xu et al. [3]. Briefly, serum samples were collected and stored at −80°C at a biobank in Lund (Sweden) until lab analysis. Total non-isomer-specific PFHxS, PFOS, and PFOA levels were analyzed using liquid chromatography–tandem mass spectrometry at the Division of Occupational and Environmental Medicine in Lund (Sweden). The laboratory is qualified for PFAS analysis as an HBM4EU laboratory and participates in a quality control program for PFAS analysis (University of Erlangen-Nuremberg, Germany).

DNA Methylation Analysis

DNA was extracted from the whole blood using Chemagic Maxiprep (PerkinElmer, Rodgau, Germany). The blood samples were taken simultaneously as the serum samples for PFAS analysis. DNA quality was evaluated using a NanoDrop spectrophotometer (NanoDrop Products, Wilmington, DE), based on the 260/280 nm ratio (all samples had a ratio ∼1.8, which is generally accepted as “pure” for DNA). Next, 500 ng DNA was bisulfite-treated using the EZ DNA Methylation kit (Zymo Research, Irvine, CA). DNA samples were randomized for distribution in a 96-well analysis plate. The entire bisulfite-treated DNA samples were used for the DNA methylation analysis. Genome-wide DNA methylation was determined at SNP&Seq, Uppsala, Sweden, using the Infinium MethylationEPIC BeadChip kit (lllumina, San Diego, CA) to analyze approximately 850 000 specific CpG sites in the genome. All beadchips were from the same batch.

DNA Methylation Preprocessing

The statistical software R was used for the analysis of DNA methylation data. ChAMP [27, 28] and minfi (using ChAMP.load) [29, 30] were used for image processing, quality control, filtering, and normalization. Normalization was performed with beta-mixture quantile (BMIQ) normalization [31]. CpGs with sinngle nucleotide polymorphisms (SNPs), non-CpG probes, and multi-hit probes were filtered out. All 63 samples performed well (the detection P-value was below 0.01 for at least  >98% of CpGs for all samples). CpGs for which the detection P-value was above 0.01 in more than 20% of the samples were removed. Probes in the X or Y chromosome were removed in the analyses including all individuals. 740 630 CpGs were retained in the analyses including all individuals, and 757 346 CpGs were retained for the sex-stratified analyses.

DNA methylation data were used to determine the proportion of different cell types (B cells, CD4+ T cells, CD8+ T cells, granulocytes, monocytes, and natural killer cells) from reference datasets for sorted samples using the function estimateCellCounts in minfi [30] as described by Houseman et al. [32]. These proportions were later evaluated as potential covariates in the statistical models. Singular value decomposition (SVD) [33] was performed in ChAMP to identify the technical and biological variables associated with DNA methylation. The SVD analyses showed that the technical variable “slide” (i.e. physical position in the analysis plate) accounted for a substantial fraction of variation (P < 0.05) in principal component 2; thus, “slide” was then included in ComBat adjustments [34, 35] using ChAMP.

Predicted Epigenetic Age

The biological age was estimated using the epigenetic “skin and blood clock,” a highly robust DNA methylation age estimator based on 391 CpGs [20]. The individual epigenetic age acceleration was then estimated as the residual by regressing the estimated biological age to the chronological age. A positive value of age acceleration indicates a faster aging than expected, i.e. an accelerating epigenetic clock.

Statistical Analysis

The SVD analyses that identified variables associated with DNA methylation showed that estimated fractions of B cells, CD4+ T cells, CD8+ T cells, and neutrophils had a P-value <0.05 in the first two principal components. Due to the strong correlation between the fractions of the different cell types, the model described below was only adjusted for neutrophil fraction, which was the cell type that showed the strongest association with DNA methylation (data not shown).

Differentially methylated positions (DMPs) were evaluated by fitting a linear regression model to each CpG using the R package limma [36], with adjustments for the estimated fraction of neutrophils. M-values (the log2 ratios of the intensities of methylated probe versus unmethylated probe) were used as the dependent variable [37]. Empirical Bayes smoothing was applied to the standard errors with a robust selection against outlier sample variances. Furthermore, P-values were adjusted for multiple comparisons for all CpGs by the Benjamini–Hochberg false discovery rate (FDR) method to obtain q-values. CpGs with a q-value of 0.05 or lower were considered being DMPs. Analyses were performed for 1) all individuals and 2) stratified for sex.

Differentially methylated regions (DMRs) were evaluated using the R package DMRCate [38]. Gaussian kernel bandwidth for smoothed-function estimation was set to 1000 nucleotides, and the scaling factor for bandwidth was set to 2. DMRs that were at least two CpGs long were included. DMRs are ranked by Fisher’s multiple comparison statistics. Regions with a q-value of 0.05 or lower were considered being DMRs.

The correlation between biological age and chronological age was evaluated using Spearman correlation. The difference in epigenetic age acceleration between the high-exposure and control groups was analyzed using the Mann–Whitney U test. A P-value of 0.05 or lower was considered statistically significant.

Gene Ontology and KEGG Analyses

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses [GO/KEGG accessed May 2021] [39-41] were run for genes that had CpGs with a q-value <0.1 in the DMP analyses. GO and KEGG analyses were performed using the gometh function in the R package missMethyl, which considers the differing number of probes per gene present on the array.

Results

Table 1 shows the characteristics and serum levels of PFASs of the 63 children in this study. Even though we selected study participants based on their PFOS level, there were large contrasts between the high-exposure group and the control group in serum levels of PFOA and PFHxS as well.

Table 1:

Descriptive statistics [count (%) or median (range)] in the control group and the high-exposure group

Variable	Control group (n = 32)	High-exposure group (n = 31)
Sex (% girls)	16 (50)	16 (52)
Age (years)	9 (6–11)	9 (6–11)
BMI (kg/m2)	17 (15–21)	17 (14–23)
PFOS (ng/ml)	2.8 (1.0–4.7)	295 (190–464)
PFHxS (ng/ml)	0.5 (0.3–1.6)	277 (147–410)
PFOA (ng/ml)	1.4 (0.7–3.5)	19 (12–32)

Epigenome-wide Analysis of DMPs and DMRs

We first evaluated the DMPs among all study participants. Comparing the high-exposure and control groups, we identified 12 DMPs (q < 0.05, Table 2), and there were 364 CpGs with a q-value <0.1 (Supplemental material, Excel Supplemental Table S1). Most of the 12 DMPs (75%) were hypermethylated in the high-exposure group compared to the control group (Table 2).

Table 2:

The 12 DMPs (q < 0.05) between the high-exposure group and control group

CpG	Chr	Positiona	Gene	Gene name	2^logFC (95% CI)b	q-Value
cg15732078	10	86097061	CCSER2	Coiled-coil serine-rich protein 2	1.36 (1.23, 1.49)	0.025
cg16593554	18	74813401	MBP	Myelin basic protein	1.47 (1.30, 1.67)	0.025
cg02556924	15	101661693	NAc		0.78 (0.72, 0.85)	0.033
cg01899620	13	113689422	MCF2L	MCF.2 cell line derived transforming sequence like	0.77 (0.71, 0.85)	0.034
cg23696710	11	17260094	NA		1.42 (1.27, 1.61)	0.034
cg05371502	8	67522986	MYBL1	MYB proto-oncogene like 1	1.30 (1.18, 1.42)	0.046
cg25840780	20	50156843	NFATC2	Nuclear factor of activated T cells 2	1.42 (1.25, 1.61)	0.046
cg24295885	4	57267050	PPAT	Phosphoribosyl pyrophosphate amidotransferase	1.31 (1.19, 1.45)	0.046
cg02051878	15	74502715	STRA6	Stimulated by retinoic acid 6	1.28 (1.17, 1.40)	0.046
cg07503415	6	44171877	NA		1.26 (1.16, 1.36)	0.046
cg13562278	3	31207280	NA		1.32 (1.19, 1.45)	0.046
cg23821954	18	46497074	NA		0.74 (0.66, 0.82)	0.046

Abbreviations: Chr, chromosome; CI, confidence interval; q-value, FDR-adjusted P-value using the Benjamini–Hochberg method; logFC, binary logarithmic fold change.

Chromosomal coordinates of the CpG (Build 37).

2^logFC, two to the power of logFC. LogFC denotes β1 from the following robust regression model: M-value = β1 × group comparison + β2 × estimated fraction of neutrophils. The control group is the reference group.

NA, not annotated, i.e. the CpG is not present in any known gene.

We then performed a DMP analysis in each sex stratum, but none of the CpGs had a q-value <1.0. Top CpG for girls was cg14111928 in lysine acetyltransferase 6B (MYST1) [2^log fold change (FC) = 1.24, 95% confidence interval (CI): 1.14–1.31, unadjusted P-value = 0.00001, q = 1.0], while top CpG for boys was cg16813892 in metal response element-binding transcription factor 2 (MTF2) (2^logFC = 1.43, 95% CI: 1.26–1.64, unadjusted P-value = 0.000002, q = 1.0). Top 10 CpGs of the sex-stratified analysis are presented in Supplemental material, Excel Supplemental Tables S2 and S3.

We then evaluated DMRs. We found seven DMRs (Fisher’s multiple comparison statistic q < 0.05) (Table 3) in the high-exposure group compared to the control group. To note, the top DMP cg15732078 in coiled-coil serine-rich protein 2 (CCSER2) was also part of one top DMR, although this was a small DMR encompassing only two CpGs. DMR analyses were also performed in each sex stratum, but no regions showed a q-value <0.05.

Table 3:

The seven DMRs (q < 0.05) between the high-exposure group and control group

Gene(s)	Gene names	Number of CpGs	Location (Chr: start, end)	Mean 2^logFC within the regiona	q-Value
LPCAT1	Lysophosphatidylcholine acyltransferase 1	5	chr5: 1503259-1503979	1.022	0.022
CCSER2	Coiled-coil serine-rich protein 2	2	chr10: 86096789-86097061	1.024	0.029
BST2, MVB12A	Bone marrow stromal cell antigen 2, multivesicular body subunit 12A	13	chr19: 17515486-17517762	1.013	0.030
RASAL2	RAS protein activator like 2	2	chr1: 178244359-178245008	1.028	0.033
TLDC2	TBC/LysM-associated domain containing 2	8	chr20: 35503983-35504553	1.020	0.043
MYO7A	Myosin VIIA	3	chr6: 109035178-109036145	1.05	0.045
ITPR1	Inositol 1,4,5-trisphosphate receptor type 1	3	chr3: 4790361-4790790	1.03	0.050

Abbreviations: Chr: chromosome; q-value, FDR-adjusted P-value using the Benjamini–Hochberg method; logFC, binary logarithmic fold change.

2^logFC, two to the power of logFC. LogFC denotes β1 from the following robust regression model: M-value = β1 × group comparison + β2 × estimated fraction of neutrophils. The control group is the reference group.

Fisher’s multiple comparison statistic q-value.

Gene Ontology/KEGG Analyses for DMPs

Next, we tested for potential enrichment of GO and KEGG terms for the CpGs with a q-value <0.1 annotated to genes. We found no significant terms after the correction for multiple testing (FDR < 0.1), and all terms had an FDR of 1.0. The top five KEGG and GO terms, ranked by P-value, are listed in Table 4. We did not perform any pathway analysis for the sex-stratified data since there was no CpG with a q-value <0.1 in the sex-stratified analyses.

Table 4:

Top KEGG and GO terms enriched by the CpGs with a q-value <0.1, ranked by P-value

Term	Term ID	P-valuea	Total number of genes in Term	Number of genes from CpGs with a q-value <0.1 in this study
Gene ontology
Regulation of developmental process	GO:0050793	0.000168	2538	35
Hemopoiesis	GO:0030097	0.000328	839	16
Osteoclast development	GO:0036035	0.000386	17	3
Positive regulation of myoblast fusion	GO:1901741	0.000474	19	3
Hematopoietic or lymphoid organ development	GO:0048534	0.000619	881	16
KEGG
Apoptosis	KEGG:04210	0.01	132	4.5
Basal transcription factors	KEGG:03022	0.02	39	2
Mineral absorption	KEGG:04978	0.04	53	2
Lysosome	KEGG:04142	0.04	123	3
Glycosphingolipid biosynthesis—globo and isoglobo series	KEGG:00603	0.08	13	1

All terms had a q-value of 1.0.

Epigenetic Age Acceleration

Lastly, we investigated whether PFAS exposure was associated with epigenetic age acceleration. There was a high correlation between the epigenetic age and chronological age in the study population (r = 0.58, P < 0.001). We did not observe any difference in epigenetic age acceleration between the high-exposure and control groups (P = 0.82 from the Mann–Whitney test).

Discussion

This is, to our knowledge, the first study to conduct an epigenome-wide analysis to examine the association between PFAS exposure and DNA methylation in a highly exposed group of children at school age. In this exploratory study with 63 children, 12 DMPs and 7 DMRs were significantly associated with high PFAS exposure. However, there was no association between PFAS exposure and epigenetic age acceleration. These observations suggest that PFAS exposure is associated with DNA methylation for some specific genes. However, future and larger studies are needed to confirm these associations and convey potential causality between PFAS exposure and DNA methylation and potentially adverse health outcomes.

Although there have been some studies evaluating PFASs and DNA methylation in newborns at background exposure levels [9-12], no studies involving children beyond the newborn stage have evaluated the influence of PFAS exposure on DNA methylation. Investigating PFAS-related DNA methylation in different age groups is essential, given that DNA methylation can change during the entire life span [42]. The most significant age-related dynamics of DNA methylation manifest in early childhood from 2 to 10 years old [43]. Thus, due to the difference in exposure levels between previous studies and ours and the age-related dynamics of DNA methylation during childhood, it is not known how comparable these findings are to studies in newborns at background exposure levels. Additionally, studies at high-exposure levels are few and restricted to adults. There were no overlaps of DMPs or DMRs with those in our previous study comprising women in the same cohort [16]. One reason for this could be aged-related differences in methylation.

In our study, RASAL2 was found to be hypermethylated in the group with high PFAS exposure compared to controls, indicating possible epigenetic silencing. RASAL2 gene encodes a protein (RASAL2) that belongs to a family of Ras guanosine triphophate (GTP)ase-activating proteins (RasGAPs), stimulating the GTPase activity of wild-type RAS [44]. This gene has been identified as a functional tumor suppressor [45], and downregulation of the RASAL2 protein has been reported in various cancers, e.g. lung, breast, bladder, and colorectal cancer [45-48]. RASAL2 hypermethylation has been associated with higher stages and grades and worse survival rates in patients with renal cell carcinoma [49]. Interestingly, a recently published case-control study observed a higher risk of renal cell carcinoma associated with PFOA exposure [50], suggesting that PFOA is a renal carcinogen [51]. Our research group found modestly increased risks of kidney cancer in the Ronneby population with high PFAS exposure through drinking water [52]. Kingsley et al. found that RAS p21 protein activator 3 (RASA3), another gene coding RasGAPs, was hypomethylated with increased maternal exposure to PFOA in cord blood samples in 44 mother–infant pairs with background exposure levels [9]. Although RASAL2 and RASA3 encode RasGAPs, they may have different biological functions; thus, a direct comparison to our findings is difficult. Future studies with a larger sample size are needed to verify our findings of hypermethylated RASAL2 after being highly exposed to PFASs.

ITPR1, a gene encoding an intracellular receptor for inositol 1,4,5-trisphosphate and mediating calcium release from the endoplasmic reticulum, was also found to be hypermethylated in the high PFAS exposure group compared to the control group. This gene has not been reported in other human studies investigating PFAS-associated epigenetic changes. However, in an animal study with PFOS-exposed rats, alteration in miRNAs encoding ITPR1 was found in neonatal rats’ brain samples, indicating that PFOS may affect calcium actions through interfering with ITPR1 and cause adverse effects on the nervous system [53].

The other genes encompassing DMPs or DMRs in our study have not been previously associated with PFAS exposure, while some of them may be of interest considering their biological function. MCF.2 cell line derived transforming sequence like (MCF2L), a gene encoding a guanine nucleotide exchange factor that interacts specifically with the GTP-bound Rac1 and plays a role in the Rho/Rac signaling pathways [54], was found to be hypomethylated in the high PFAS exposure group compared to the control group, indicating that PFAS exposure may inhibit MCF2L expression. Given that MCF2L was less expressed in osteoporosis patients than controls in one study [55] and PFASs have been associated with lower bone mineral density in mid-childhood [5], an association linking PFAS exposure, MCF2L expression, and bone health warrants further studies. In addition, osteoclast development was one GO term related to PFAS exposure, although it was not significant after FDR adjustments.

In the KEGG pathway analysis, apoptosis was the suggested pathway with the smallest P-value, although the significance disappeared after FDR correction. In animal and cell studies, PFASs were shown to induce apoptosis in cells through the production of reactive oxygen species and oxidative stress [56, 57] and the modulation of gene regulation via peroxisome proliferator-activated receptors (PPARs) and nuclear factor kappa beta (NF-κB) transcription [4, 58]. Our previous study on Ronneby women found PPARα/retinoid X receptor alpha (RXRα) activation as an enriched pathway [16]. However, the interpretation should be cautious due to the lack of statistical significance.

We did not note any association between PFAS exposure and epigenetic age acceleration among these school-age children. In highly PFAS-exposed women aged 20–45 years from the Ronneby Biomarker Cohort, no association with epigenetic age acceleration was seen either [16]. The main limitation of this study is the sample size which reduced the statistical power. Since the PFAS is an endocrine disruptor, sex differences in PFAS-related DNA methylation are plausible, but this study did not have enough power to elucidate sex-specific effects. Another limitation is the lack of information about socioeconomic status, pubertal status, parental occupation of the children, etc., which may be potential confounders. However, the control and highly exposed PFAS groups came from two closely situated municipalities with similar geodemographic statuses. Given the small sample size, we used frequency matching for age, sex, and BMI to control pubertal status. A strength of the study is the large contrast in PFAS exposure, comparing children with a general population background exposure to children with very high exposures, especially for PFOS and PFHxS.

Our study among children indicates that PFAS exposure is associated with DNA methylation changes at specific sites. Further studies with larger sample sizes are needed to validate the findings and couple the PFAS-associated changes in DNA methylation to phenotypical effects of PFASs.

In conclusion, we found that PFAS exposure was associated with DNA methylation at specific genomic positions and regions in children at school age. Most DMPs and DMRs were hypermethylated in the high-exposure group compared to the control group.

Click here for additional data file.