Environmental toxicants in breast milk of Norwegian mothers and gut bacteria composition and metabolites in their infants at 1 month

Characteristics of the study cohort

PFOS and PBDE-28 associate with infant gut microbiome α-diversity

To test which of the 28 chemicals were associated with α-diversity (Shannon diversity, Faith’s phylogenetic diversity, and number of observed sub-OTUs), we used elastic net regression, a penalized regression method to select among correlated chemicals [25, 26], adjusting for C-section, preterm delivery, maternal α-diversity 4 days after birth, and proportion of meals given through breastfeeding at 1 month. Elastic net selected 2 of 28 toxicants as the best predictors associated with less α-diversity, although more PBDEs were associated in single pollutant models (Additional file 1: Figure S3). In the unpenalized model, a one standard deviation (SD) increase of 0.5 ng/g milk lipid of PBDE-28 was associated with less Shannon diversity (− 4%, 95% confidence interval [CI]: − 7% to − 2%, relative to mean Shannon), explaining 2% of variance. In the same model, preterm delivery was associated with less Shannon diversity (− 15%, 95% CI − 23% to − 7%), any formula feeding with more (11%, 95% CI − 1% to 24%), while none of the other potential confounders were associated. A 1SD increase of 63 ng/L of perfluorooctane sulfonate (PFOS) was associated with less phylogenetic diversity (− 5%, 95% CI − 9 to − 1%), explaining 4% of variance. Of the potential confounders in that model, only C-section was significantly associated with phylogenetic diversity (− 9%, 95% CI − 18% to − 1%). PFOS was also associated with number of observed sub-OTUs (− 7%, 95% CI − 12 to − 1%) (Additional file 1: Figure S3).

Sensitivity analyses (restricting to complete case, breast milk sample collection age < 60 days, exclusive breastfeeding, and excluding extreme values) attenuated effect estimates or altered precision in some cases, but in general did not affect the overall interpretation. When restricting to term births, elastic net additionally selected PCB-167 associated with significantly less phylogenetic diversity and number of sub-OTUs (Additional file 1: Table S5). We found no material influence on results from the inclusion/exclusion of additional variables or interactions (including preterm delivery) as described in the methods.

PFOS and PCB-167 associate with infant gut microbiome β-diversity

For each toxicant, we tested whether infants exposed to “low” (< 20th), “medium” (≥ 20–≤ 80th) or relative “high” (> 80th percentile) breastmilk concentrations were more or less similar based on β-diversity distances (weighted or unweighted UniFrac), using permutational multivariate analysis of variance (PERMANOVA) to test significance, with Bonferroni correction for multiple testing. There were significant differences in community composition between PFOS exposure groups, with greater dissimilarity between high and medium exposed (p = 0.003), and between high and low exposed (p = 0.010) using unweighted UniFrac (Fig. 1a). Community composition in infants relatively highly exposed to PCB-167 was more diverse than in low exposed infants (p = 0.001) with unweighted UniFrac (Fig. 1c). There were no other significant differences in unweighted UniFrac, nor any differences in weighted UniFrac among other toxicant exposures.

Relative abundance of sub-OTUs belonging to Firmicutes differ in infants in the “high” toxicant exposure group

Next, we tested differential abundance of microbes in infants with relatively high vs. low exposure to individual chemicals, using analysis of composition of microbiomes (ANCOM) with a Benjamini-Hochberg correction, and adjusted for gestational age [27]. We detected differential abundance of some sub-OTUs in those exposed > 80th vs. < 20th percentile, and assigned lineages following the position of the sub-OTU sequence in the Greengenes reference tree and collecting taxonomic labels along the path to the root. Infants in the high exposure group had some differentially abundant sub-OTUs within Firmicutes, particularly the genus Lactobacillus (Fig. 2). Infants in the high perfluorooctanoic acid (PFOA) exposure group lacked a sub-OTU from the lineage of Lactobacillus zeae (p < 0.001) and 1.1-fold more of a sub-OTU of the genus Enterococcus (p = 0.03). Infants in the high dioxin-like PCBs exposure group had 1.8-fold lower Lactobacillus gasseri (p = 0.005) and 1.8-fold greater abundance of Clostridium perfingens (p = 0.013). High organochlorine pesticide exposure showed 0.9-fold more of a sub-OTU within the genus Streptococcus (p = 0.001), while the high PBDE-28 exposure group had 2.8-fold less Veillonella dispar (p = 0.017). The deblurred FASTA sequences are in Additional file 1: Table S3. Using a closed-reference OTU table for the same raw reads, more taxa were differentially abundant in the high breast milk-toxicant exposure group, notably PCB-167 (Additional file 1: Table S4). In addition to differential OTUs within Firmicutes, we detected lower Actinobacteria (Bifidobacterium bifidum, Corynebacterium, Eggerthella lenta) and Bacteroidetes (Bacteroides fragilis and other unidentified species of Bacteroides).

Dioxin-like PCB-167 associates with relative abundance and more α-diversity of the predicted functional profile

Since a taxonomic group might be replaced by another while preserving its functional contribution to the microbial community, we used PICRUSt to infer the microbiome functional profile based on clusters of orthologous groups of proteins (COG) and alternatively by players in metabolic pathways of the Kyoto Encyclopedia of Genes and Genomes (KEGG) [28]. We tested relative abundance of the various features in “high” vs. “low” exposure groups, and compared α- and β-diversity of the functional profile across exposure groups. Infants exposed to milk with higher PCB-167 had a larger functional spectrum: significantly higher Shannon diversity in COG (Fig. 3a) and KEGG (Fig. 4a). The microbiome of higher exposed infants also had 0.2-fold enrichment of “cell motility” associated predicted proteins (p = 0.003) and 0.05-fold less abundant “carbohydrate transport and metabolism” (p = 0.006) (Fig. 3b). They also had upregulated predicted proteins of “cellular processes,” “human diseases,” and “unclassified” pathways, and less abundance of predicted proteins involved in “genetic information processing” and “metabolism” pathways (Fig. 4c). Seven enzymes involved with amino acid synthesis, one carbon pool by folate, and lipid metabolism pathways were less abundant in the high PCB-167 exposure community compared with low (Fig. 4d). Other toxicants were less consistently associated with differences in the predicted metagenome (Additional file 1: Figures S4 and S5).

Toxicants associate with lower concentrations of short-chain fatty acids, except PCB-167 and PFOA, which associate with higher concentrations

We studied SCFAs measured in fecal samples (n = 70). At 1 month of age, infant fecal samples were dominated by acetic (91.7 ± 7.4%) and propionic acids (5.4 ± 5.1%) (Additional file 1: Table S2). Elastic net selected a number of toxicants as predictors of acetic and propionic acid in models adjusting for confounders (C-section, preterm delivery, maternal diversity 4 days after birth, and proportion of meals given through breastfeeding at 1 month). Multipollutant models explained 20% and 34% of variance in acetic and propionic acid, respectively, increasing to 25% and 48% when the confounders were included (Fig. 5). PCB-209 and PBDE-47 were associated with less acetic acid (− 15% [95% CI − 29% to − 0.4%] and − 11% [95% CI − 31% to 9%], respectively, statistically non-significant for the latter). Brominated flame retardants were also associated with less propionic acid (− 24% [95% 95% CI − 35% to − 13%] for PBDE-28, and − 16% [95% CI − 35% to 3%] for PBDE-47), as was PCB-170 (− 40% [95% CI − 102% to 21%], statistically non-significant for the latter two). Conversely, dioxin-like PCB-167 and PFOA were associated with more acetic acid (22%, 95% CI 8% to 35%) and propionic acid (61%, 95% CI 35% to 87%). Associations were generally imprecise for other SCFAs (low concentrations and large proportion below LOD) (Additional file 1: Figure S6).

PBDEs

PBDE-28 was associated with 4% less Shannon diversity. Additionally, infants highly exposed to PBDE-28 had samples with relatively less abundance of Veillonella. PBDE-28 and PBDE-47 were also associated with less propionic and acetic acids, the signaling metabolites that play an important role in immune system development [29]. We speculate that these lower levels of SCFAs may in part be explained by relatively lower abundance of Veillonella in the high exposed, as these bacteria are known to utilize lactate in the intestine, fermenting it to both propionic and acetic acids [30]. The potential of brominated flame retardants for disrupting microbiome composition and function has been demonstrated experimentally. In a small study, mice fed PBDE-47 or PBDE-99 had decreased microbial richness, differential abundance of some taxa, and disrupted bile acid metabolism compared with control mice [31]. The acute doses of PBDEs by oral gavage (48.5 mg/kg body weight for BDE-47), although higher than environmental exposures, were estimated to result in circulating levels at concentrations similar to that found in human populations. PBDE-71 exposure at environmentally realistic concentrations led to decreased bacterial diversity of the zebra fish gut microbiome and disrupted metabolic functions such as energy metabolism, virulence, respiration, cell division, cell signaling, and stress response [32]. However, we also note that for Shannon diversity, the differences associated with PBDE-28 were much less than that from either preterm delivery (15% less than term delivery) or any formula feeding (11% more compared to exclusive breastfeeding).

PFASs

PFOS was associated with 5% less microbiome α-diversity. By comparison, C-section, known to disrupt the microbiome, was associated with 9% less phylogenetic diversity in this population (while preterm delivery and full formula feeding were not associated). The PFOS finding was robust to sensitivity analyses including the addition/exclusion of other potential confounders and restriction to term births. Furthermore, there was greater dissimilarity between the communities in the low and high PFOS exposure groups than within them. In an experimental setting, mice fed PFOS through oral gavage had a significant decrease in abundance of bacteria [33]. However, another study of mice with dietary exposure to PFOS at doses corresponding to general population and occupational exposure found no significant differences in gut microbial diversity relative to the control group [34]. They did observe differential abundance of bacteria within Firmicutes and Bacteroidetes, and high-dose PFOS exposure significantly induced butanoate metabolism. Here, we did not find a statistically significant association between PFOS and the SCFA metabolites. In contrast, higher PFOA exposure was associated with both more propionic acid, absence of a sub-OTU within the genus Lactobacillus, and greater relative abundance of a sub-OTU of Enteroccocus. In rodents, propionic acid enhances adipocyte differentiation of 3 T3-L1 pre-adipocytes via increased expression of GPR43 and peroxisome proliferator-activated receptor γ (PPARγ) [35], as does PFOA [36]. PFASs continue to be a concern for human health, as evidenced by the recent lowering of the European Food Safety Authority’s tolerable weekly intake [37], and the interaction between these compounds and the gut microbiome requires further investigation.

Dioxin-like PCB-167 was associated with greater β-diversity of the gut microbiome, enriched metabolic activity, and more acetic acid. This may indicate a more heterogeneous response to exposure, leading to a larger functional spectrum in more exposed communities. In experimental studies, the AhR mediates dioxin-like compounds’ toxicity and exposure increased butyrate and propionate in AhR^+/+ but not AhR^−/− mice [38], while larval exposure to PCB-126 increased phylogenetic diversity in the frog gut [39]. Exposure may perturb community structure, selecting for specialized/more tolerant microbes able to degrade the chemical [40]. This could allow colonization by potential opportunistic microbes; here, infants with higher PCB-105 exposure had greater abundance of Clostridium perfringens. In mice, dioxin-elicited changes in the host decreased B. fragilis [41], which, in our study was relatively lower in abundance in the higher breastmilk toxicants exposure group, as were some Lactobacillus. Due to financial constraints, we only had measured levels of the less toxic of the dioxin-like PCBs; however, we expect these to be moderately correlated with dioxins and other more toxic dioxin-like compounds [3]. Given our findings and experimental studies, a more detailed investigation of dioxins and the human gut microbiota is warranted.

Organochlorine pesticides

There were fewer associations with organochlorine pesticides. High oxychlordane and dichlorodiphenyltrichloroethane (p,p′-DDT) exposed communities were associated with greater relative abundance of a sub-OTU of the genus Streptococcus, while the metabolite dichlorodiphenyldichloroethylene (p,p′-DDE) was selected as a predictor of less propionic acid. In rats, coliform bacteria metabolize p,p′-DDT to p,p′-DDD [42], which is probably relevant for humans.

Non-dioxin like PCBs

These compounds were associated with gut microbial function (decreased acetic and propionic acid), with less evidence for a disturbed community composition. Limited experimental evidence reports that mice orally exposed to non-dioxin-like PCBs had decreased abundance of gut bacteria, primarily Proteobacteria [43].

Strengths and limitations

This study is based on a prospective birth cohort with rich questionnaire data to assess potential confounding. We had extensive exposure assessment of persistent chemicals at the time when fecal samples were collected, providing substantial information on the multiple breast milk toxicants to which gut bacteria were exposed. Our methods allowed for adjustment for confounding from co-occurring toxicants in linear regressions, although we could not consider the exposure profile of co-occurring toxicants when assessing low, medium, and high exposure groups of individual toxicants. PCB-167 influenced a number of microbiome metrics; however, given the high correlation with other PCBs (r = 0.77 – r = 0.92), these may represent common effects. Additionally, other unmeasured compounds may be more influential confounders (i.e., arsenic [44]). We assessed breast milk concentrations as direct exposure for the bacteria and not child blood concentrations (although they are correlated in early life [45]), which could influence gut microbiota through host physiology.

Using Deblur increases resolution, and reduces false positive annotations [24]. We detected more differentially abundant taxa using a Greengenes closed-reference table, thus some taxonomic differences may be expressed in species- to genus-level rather than sub-OTU level. This could be followed up with in vitro studies of toxicant effects on strains and species from the same genus.

We had a reasonable sample size for the α-diversity analyses; however, SCFA analyses were in 70 infants, and should be interpreted cautiously. Arsenic and diazinon perturb the gut microbiata in a sex-specific manner in mice [46, 47]; there was no interaction between toxicants and sex on α-diversity, but for SCFAs we could not test this.

All toxicant classes were associated with some alterations in composition and function, but not consistently across all metrics, or toxicants. This could be due to specific chemicals only affecting particular aspects of the microbiome, misclassification of exposure, statistical methods, or chance findings. The sensitivity and resilience of gut microbiota to environmental toxicants has been demonstrated in fish [48], and it may be difficult to detect small, transient effects in an observational design.

Preterm babies, whose gut microbiota could be more susceptible to the effect of toxicants due to immaturity of their immune system, were over-represented. In the linear regressions, adjusting for preterm delivery did not influence the toxicant effect estimates, although restricting to term births (22.5% reduction in study population) affected the interpretation of PCB-167, which became associated with decreased α-diversity.

Forty percent of mothers in the NoMIC cohort did not deliver milk, and these women were not breastfeeding (15%), or preterm (39% vs. 26% of those who delivered milk), or if they were breastfeeding exclusively breastfed for a shorter period (2.3 vs. 4.2 months), possibly indicating difficulties breastfeeding. However, there were no significant differences in infant gut diversity in the full cohort compared with our study population, so we do not expect this to bias our results.

Breast milk is an evolutionary development containing numerous specialized bioactive substances. Oligosaccharides, milk lipids, secretory IgA, and hormones are released into the milk and are uniquely adapted to the individual baby in response to the mother’s living conditions. Protection against infections and a small beneficial effect on IQ are well documented health benefits of breastfeeding [2]. Furthermore, there is a growing understanding of the role of breast milk bacteria in seeding the infant gut [49]. Here, maternal-reported formula feeding at 1 month (i.e., non-exclusive breastfeeding) was associated with more Shannon diversity in their infants (with stronger effect size than from the toxicants), and the implications of such on child health should be investigated in a targeted study. Although the toxicant concentrations were lipid-normalized, we recognize that there is a complex relation between lipophilic chemicals and fatty acids, with a possible effect on the gut microbiome. Potential interactions between lipids, bioactive substances, bacteria, and environmental toxicants in breast milk are an avenue for future research.

The milk was sampled between 2002 and 2006 from women in Norway. In World Health Organization surveys 2005–2010, breast milk from Norway had higher levels of dioxin-like PCBs (expressed as toxic equivalency factor, 3 pg TEQs/g lipid) and sum of 6 indicator PCBs (62 ng/g lipid) than the less industrialized countries of the southern hemisphere (i.e., Australia 1.8 pg TEQs/g lipid and < 20 ng/g lipid), but not among the highest (i.e., Czech Republic 7 pg TEQs/g and 380 ng/g lipid). However, regardless of level, surveyed countries had levels of dioxin-like PCBs and sum of PCBs in human milk at one to two orders of magnitude above those considered toxicologically safe in early childhood [3, 50]. By contrast, the sum of DDT in Scandinavian breast milk was the lowest (< 100 μg/kg lipid), with other European countries comparatively higher (i.e., Czech Republic 130 μg/kg lipid, and the highest in the tropical countries using DDT for vector control i.e. India > 1000 μg/kg lipid). PBDE levels are also relatively low in Norway [51]. Due to restrictions, these chemicals are in decline [52], although exposure continues through dust and food, especially in countries with lower environmental controls. These findings are relevant for the general population due to continued contamination of fish and meat.

Study population and data collection

The Norwegian Microbiota Cohort (NoMIC) is a prospective birth cohort [21–23]. Mothers were recruited at the maternity ward of Østfold county hospital (2002–2005), two consecutive term births per preterm delivery. Fluency in Norwegian and residency in the county were inclusion criteria. Mothers were asked to collect and freeze one fecal sample from themselves at 4 days postpartum, as well as samples from their infants when they were 4, 10, 30, 120, 365, and 730 days old. Participants were asked to collect by hand a 25-ml breast milk sample each morning for eight consecutive days, between 2 weeks and 2 months postpartum [20], but minor changes in sampling protocol were also accepted. Sampling was undertaken on multiple days to reduce within-subject variability in estimated level of exposure. The milk was stored in a 250-ml container in the freezer. When the mothers had filled the container, the milk samples and fecal samples were collected by study personnel, kept frozen during transport to the Norwegian Institute of Public Health (NIPH), and stored at − 20 °C upon arrival. DNA was extracted after all samples were collected [21, 53]. Six hundred one women agreed to participate, 89% returned fecal samples, leaving a cohort of 552 children. Three hundred twenty-one mothers also delivered breastmilk samples with measured toxicants, corresponding to 333 children (including multiple births); 5 did not have microbiome information due to lost samples leaving 328 children with measurement of toxicants and gut microbiome diversity at any time point. We focused on 1 month, a sensitive period when the microbiome undergoes rapid development [54], and 87% of women were exclusively breastfeeding; whereas feces sampled at later time points could be influenced by other factors such as antibiotic use, introduction of solid food, and diet. Three hundred seven infants had both chemicals measured and fecal samples at 1 month. We excluded twins and triplets who may have different feeding patterns and thus the toxicants sampled in milk were not representative of their exposure (n = 26), infants whose mothers reported no breastfeeding at 1 month (n = 3), or infants with antibiotics use 14 days prior to fecal sampling (n = 3). Two hundred sixty-seven infants were in our α-diversity analysis (Additional file 1: Figure S7). To assess differences in microbiome composition between infants grouped by low (< 20th), medium (≥ 20th–< 80th), and high (≥ 80th percentile) breast milk toxicant exposure, we restricted to exclusively breastfed babies, since we could not adjust for confounders in those analyses (n = 239). SCFAs were not available for all children due to lack of sample volume, thus we studied microbiota function in a subset of participants (n = 70).

We obtained information on gestational age, maternal smoking, and birth weight and length through the Norwegian Medical Birth Registry, and additional important covariates, including maternal education, antibiotics use, breastfeeding, and C-section, from the 1-month questionnaire.

Exposure variables

Mothers sampled their breast milk at a mean (SD) age of 31.4 (19.9) days using the WHO protocol [20]. Due to financial constraints, 28 chemicals were analyzed in breast milk samples in 3 laboratories: non-dioxin-like PCBs, mono-ortho dioxin-like PCBs, organochlorine pesticides, PBDEs, and PFAS (Additional file 1: Table S6). University of Life Sciences-NMBU measured PCBs and organochlorine pesticides in 15 samples using liquid-liquid extraction, gravimetrical lipid determination, and clean-up with sulfuric acid [53, 55, 56]. Following this, the laboratory at The Department of Environmental Exposure and Epidemiology, Norwegian Institute of Public Health (NIPH) established methods and analyzed the lipophilic chemicals using liquid-liquid extraction and gas chromatography–mass spectrometry (GC/MS) with negative chemical ionization [57, 58]. The majority of PFAS samples were measured in breast milk using high-performance liquid chromatography/tandem mass spectrometry (LC-MS/MS) at the NIPH [57, 59], with additional samples measured at Vrije University, Institute for Environmental Studies [60]. Measured concentrations of the lipophilic chemicals were normalized by dividing by total lipid content in the specific milk sample. We replaced values below the limit of detection (LODs) by a randomly imputed number between zero and LOD.

Outcomes

Mothers were in close contact with health personnel and reminded to collect fecal samples at 1 month using a standard protocol.

Covariates

We selected potential confounding factors a priori using directed acyclic graphs (Additional file 1: Figure S8). The minimum adjustment set to assess the effect of breast milk toxicants on gut diversity/SCFAs at 1 month was proportion of meals given through breast milk (vs. formula feeding, continuous 0–1), preterm delivery (Yes/No), maternal gut α-diversity, and C-section (Yes/No).

Statistics

For microbiome α-diversity analyses, we imputed missing values for exposures and covariates using multiple imputation by chained equations to generate 20 imputed data sets [68, 69]. Correlations between exposures were assessed using Spearman’s rank correlation coefficients. To assess associations between breastmilk toxicants and gut microbiota α-diversity, we adopted two regression approaches, in which we standardized exposures to one SD, and adjusted for identified covariates. First, to select among individual toxicants, we used elastic net regression modeling, a hybrid penalized method robust to extreme correlations among the predictors [25, 26]. We selected α = 0.9 and optimized λ using tenfold cross-validation repeated 10 times based on minimum standard error, and unpenalized covariates [70], and repeated in each of the 20 multiply imputed datasets, considering the exposures which were selected (β ≠ 0) in more than half of the models as noteworthy [71]. We then used generalized linear models to obtain unbiased (unpenalized) estimates and assessed all pollutants individually for comparison, and then fitted multipollutant models with the elastic net-selected exposures.

We investigated group differences for low (< 20th), medium (≥ 20–≤ 80th), and high (> 80th percentile) breast milk toxicants. First, we assessed β-diversity using weighted and unweighted UniFrac, testing significance of pairwise groups with PERMANOVA [72]. Second, we investigated differences in sub-OTU abundance between low vs. high groups using the analysis of composition of microbiomes (ANCOM) framework [27]. ANCOM accounts for the compositional nature of the taxa relative abundances and is based on the analysis of difference in pairwise log-ratios of microbial OTU abundances/relative abundances, between comparison groups of interest. For each taxon, we computed a statistic indicating the number of significantly different pairwise log-ratios while controlling for false discoveries. We applied ANCOM with a Benjamini-Hochberg correction at 5% level of significance, and adjusted for gestational age. For comparison with other studies, we assigned lineages to the identified differentially abundant sub-OTUs. Instead of using machine learning approaches like classifying against the RDP, we used the phylogenetic tree produced for diversity computation and its assigned Greengenes taxonomy labels to obtain lineages for the sub-OTUs: For every sub-OTU sequence, we started from the inserted sub-OTU tip and followed the path up to the root while collecting taxonomic labels along this path. Third, we analyzed the predicted metagenome. We tested for differentially abundant functions using a discrete false-discovery rate correction [73]. For the same relative abundances, we computed α-diversity of the normalized relative PICRUSt abundances by applying the Shannon metric and used two-sided Mann-Whitney tests to check for significant differences between the three exposure groups. We then computed β-diversity distances for the same relative abundances via the Bray-Curtis metric and performed PERMANOVA tests with 9999 permutations to check for statistically significant differences within vs. between the groups of “high,” “medium,” and “low” labeled samples. We used Bonferroni correction for multiple hypothesis testing with p < 0.05 to consider two groups as different.

Finally, we tested the relation between toxicants and SCFAs, using elastic net regression and generalized linear models, with a natural logarithm to transform the SCFAs and adjusting for confounders. We did not include maternal gut diversity in the SCFA analyses, as there were 56% missing and the multiple imputation models for the SCFAs would not converge. All regression models were tested for and met the assumptions of normality, homoscedasticity, and linearity.

We used STATA 14.0 for multiple imputation and generalized linear regression and R programme version 3.2 [74] for ANCOM and elastic net (using the glmnet package [25]). We used scikit-bio 0.5.1 for PERMANOVA tests.