Study population and data collection
The Norwegian Microbiota Cohort (NoMIC) is a prospective birth cohort [21,22,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.
Sequencing and data processing
We extracted DNA using the Earth Microbiome Project protocol: (http://press.igsb.anl.gov/earthmicrobiome/emp-standard-protocols/dna-extraction-protocol/). We sequenced 100 nt from the V4 region of the 16S rRNA gene with the Illumina HiSeq instrument. We used a recently developed sub-operational-taxonomic-unit approach, Deblur, which uses error profiles to obtain putative error-free sequences from Illumina sequencing platforms [24]. By removing noise, Deblur gives a higher resolution than OTU-based analyses or analyses of raw sequence data, and because it is reference free, it may pick up sequences of novel bacteria that are not represented in existing databases. To control for variation in sequencing coverage, the data were rarified at a depth of 20,000 sequences per sample, which lead to the removal of 8 samples. Data processing was performed in the Quantitative Insights Into Microbial Ecology (QIIME) pipeline version 1.9.1. [61]. More detailed information on DNA extraction, sequencing, and data processing is provided in the Supplemental Material.
Gut microbiota composition: α-diversity, β-diversity, differential abundance of taxa
We used three α-diversity measures: (i) Shannon diversity, the total number of species (species richness) weighted for their relative abundances (species evenness); (ii) Faith’s phylogenetic diversity, the amount or proportion of branch length in a phylogenetic tree that leads to different organisms (species richness); and (iii) the number of observed unique sub-OTUs. To study variation in diversity in the bacterial community in the exclusively breastfed children based on low, medium, or high toxicant concentrations, we calculated β-diversity using unweighted and weighted UniFrac [62]. In order to obtain a phylogenetic tree for diversity computation, we used Qiime2’s fragment-insertion [63] to phylogenetically place the sub-OTU sequences into the reference Greengenes 13.8 99% identity tree [64].
We also tested for differentially abundant taxa (sub-OTUs) as described below.
SCFAs
Two laboratories analyzed fecal samples for eight SCFAs using published analytical methods [65,66,67]. Briefly, distillates of fecal material were analyzed with gas chromatography and quantified using flame ionization detection. We assessed SCFAs with > 50% above LOD: acetic, propionic, n-butyric, i-butyric, and i-valeric acids (Additional file 1: Table S2).
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.
Sensitivity of regression analyses
We tested the sensitivity of regression model estimates: restricting to complete case, breast milk sample collection age < 60 days, exclusive breastfeeding, term births, ln-transformed exposures, and excluding extreme values. We tested the inclusion/exclusion of: household pets, infant antibiotics in the first 2 weeks of life (as those with antibiotics in the 2 weeks prior to sampling were excluded from the analyses), parity, smoking at the start of pregnancy, maternal BMI, and education. We tested interactions between the elastic net-selected compounds sex, maternal BMI and preterm birth for the α-diversity models.