Maternal Human Milk Oligosaccharide Profile Modulates the Impact of an Intervention with Iron and Galacto-Oligosaccharides in Kenyan Infants

There is little data on human milk oligosaccharide (HMO) composition in Sub-Saharan Africa. Iron fortificants adversely affect the infant gut microbiota, while co-provision of prebiotic galacto-oligosaccharides (GOS) mitigates most of the adverse effects. Whether variations in maternal HMO profile can influence the infant response to iron and/or GOS fortificants is unknown. The aim of this study was to determine HMO profiles and the secretor/non-secretor phenotype of lactating Kenyan mothers and investigate their effects on the maternal and infant gut microbiota, and on the infant response to a fortification intervention with 5 mg iron (2.5 mg as sodium iron ethylenediaminetetraacetate and 2.5 mg as ferrous fumarate) and 7.5 g GOS. We studied mother–infant pairs (n = 80) participating in a 4-month intervention trial in which the infants (aged 6.5–9.5 months) received daily a micronutrient powder without iron, with iron or with iron and GOS. We assessed: (1) maternal secretor status and HMO composition; (2) effects of secretor status on the maternal and infant gut microbiota in a cross-sectional analysis at baseline of the intervention trial; and (3) interactions between secretor status and intervention groups during the intervention trial on the infant gut microbiota, gut inflammation, iron status, growth and infectious morbidity. Secretor prevalence was 72% and HMOs differed between secretors and non-secretors and over time of lactation. Secretor status did not predict the baseline composition of the maternal and infant gut microbiota. There was a secretor-status-by-intervention-group interaction on Bifidobacterium (p = 0.021), Z-scores for length-for-age (p = 0.022) and weight-for-age (p = 0.018), and soluble transferrin receptor (p = 0.041). In the no iron group, longitudinal prevalence of diarrhea was higher among infants of non-secretors (23.8%) than of secretors (10.4%) (p = 0.001). In conclusion, HMO profile may modulate the infant gut microbiota response to fortificant iron; compared to infants of secretor mothers, infants of non-secretor mothers may be more vulnerable to the adverse effect of iron but also benefit more from the co-provision of GOS.


. Fecal Samples
The following protocol was used for DNA extraction of infant and maternal fecal samples. Samples were first thawed at 4 °C. Then in a 2.0 mL screw-cap tube containing 0.5 g of 0.1 mm sterilized zirconia beads, 250 (± 10%) mg of feces and 700 µL S.T.A.R. buffer (Roche, Indianapolis, IN, USA) were added. The FastPrep instrument (MP Biomedicals, Santa Ana, CA, USA) was used for lysis at 5.5 ms for 3 times 1 min at room temperature. Thereafter samples were incubated, while shaking at 100 rpm and 95 °C for 15 min. The samples were then centrifuged at 16000 g for 5 min at 4 °C. The collected supernatant was kept on ice, while another lysis round as described above, except that only 350 µL S.T.A.R. buffer was added, was done with the remaining stool pellet. The supernatant kept on ice was then pooled with the supernatant from the second lysis round. Purification of DNA was performed on the automated Maxwell instrument (Promega, Madison, WI, USA) by applying the Maxwell 16 Tissue LEV Total RNA Purification Kit (Promega) according to the manufacturer's instructions. To the first well of the Maxwell cartridge 250 µL of the supernatant was added and finally DNA was eluted with 50 µL of RNAse/DNAse free water.
Using a 2-step PCR, barcoded amplicons from the V3-V4 region of 16S rRNA genes were generated. For initial amplification of the V3-V4 part of the 16S rRNA we used universal primers with the following sequences: forward primer, '5-CCTACGGGAGGCAGCAG-3' (broadly conserved bacterial primer 357F); reverse primer, '5-TACNVGGGTATCTAAKCC' (broadly conserved bacterial primer (with adaptations) 802R) appended with Illumina adaptor sequences. The PCR amplification mixture contained: 1 µL fecal sample DNA, 1 µL barcoded forward primer, 15 µL master mix (1 µL KOD Hot Start DNA Polymerase (1 U/µL; Novagen, Madison, WI, USA), 5 µL KOD-buffer (10×), 3 µL MgSO4 (25 mM), 5 µL dNTP mix (2 mM each)), 1 µL (10 µM) of reverse primer and 33 µL sterile water (total volume 50 µL). PCR conditions were: 95 °C for 2 min followed by 35 cycles of 95 °C for 20 sec, 55 °C for 10 sec, and 70 °C for 15 sec. We then purified the approximately 500 bp PCR amplicons using the MSB Spin PCRapace kit (Invitek, Berlin, Germany). Subsequently we checked concentration and quality with a Qubit fluorometer (Thermo Fisher Scientific). For the second PCR in combination with sample-specific barcoded primers, purified PCR products were shipped to BaseClear BV (Leiden, The Netherlands). PCR products were purified, checked on a Bioanalyzer (Agilent) and quantified. This was followed by multiplexing, clustering and sequencing on an Illumina MiSeq with the paired-end (2x) 300 bp protocol and indexing. We analyzed the sequencing run with the Illumina CASAVA pipeline (v1.8.3) with de-multiplexing based on sample-specific barcodes. From the raw sequencing data we removed the sequence reads of too low quality (only "passing filter" reads were selected) and discarded reads containing adaptor sequences or PhiX control with an in-house filtering protocol. On the remaining reads we performed a quality assessment using the FASTQC quality control tool version 0.10.0.
Selected enteropathogenic bacteria (Clostridium difficile, Clostridium perfringens, enterohemorrhagic Escherichia coli with shiga toxin 1 (EHEC stx1), enterohemorrhagic E. coli with shiga toxin 2 (EHEC stx2), enteropathogenic E. coli with the attaching and effacing gene (EPEC eaeA), enterotoxigenic E. coli with heat-stable enterotoxin (ETEC ST), enterotoxigenic E. coli with heat-labile enterotoxin (ETEC LT), Salmonella spp. and Staphylococcus aureus) were targeted in infant and maternal fecal DNA using quantitative PCR (qPCR). Primers used for qPCR are given in Table  S2. Pre-amplification for preparation of the DNA for qPCR was performed using the Fluidigm Preamp Master Mix (Fluidigm, PN 100-5580, California, United States), primers shown in Table S2, and the recommended 12 cycles as described by the manufacturer (Fluidigm Quick Reference for Pre-amplification of cDNA for Gene Expression with Delta Gene Assay, PN 100-5875 C1, Fluidigm). We diluted the final reaction products 5-fold. Sample pre-mix and assay mix were prepared and the qPCR was run with BioMark 96.96 Gene Expression Dynamic Arrays (PN BMK-M-96.69, Fluidigm, South San Francisco, USA) as described by the manufacturer (Advanced Development Protocol 41: Using EvaGreen DNA Binding Dye for Gene Expression with the 48.48 and 96.96 Dynamic Array IFCs, Fluidigm). Target genes were amplified from representative strains to generate standards. Amplicons were then purified, cloned into PGEMT Easy Vector and heterologously expressed in E. coli according to supplier instructions (Promega). Concentrations of the plasmids were measured using a Qubit dsDNA BR assay kit (Q32850, Thermo Fisher Scientific) in triplicate on a Spark M10 plate reader (Tecan Group Ltd., Maennedorf, Switzerland). 10-fold dilutions of linearised plasmids harbouring the gene of interest were used to generate standard curves.
We analyzed 16S rRNA gene sequences using a workflow based on Qiime 1.8 [7]. We performed operational taxonomic unit (OTU) clustering (open reference), taxonomic assignment and reference alignment with the pick_open_reference_otus.py workflow script of Qiime, using uclust as clustering method (97% identity) and GreenGenes v13.8 as reference database for taxonomic assignment. Reference-based chimera removal was done with Uchime [8]. The RDP classifier version 2.2 was performed for taxonomic classification [9]. For the Bifidobacterium genus specific analyses we re-clustered OTUs assigned to Bifidobacterium at the 99% identity level and we analyzed these further on composition and diversity. We performed statistical tests as implemented in SciPy (https://www.scipy.org/), downstream of the Qiime-based workflow.
Using the Fluidigm Real-Time PCR Analysis Software we processed qPCR data, including melting curve analysis. For further analyses of the qPCR data we used Excel (Microsoft Office 2010).
Linear regression analysis of the CT (threshold cycle) values versus the amounts of the template DNA of the standard (log gene copies/µL) was performed to generate the standard curves. For each linear regression we calculated the goodness of fit (r 2 ) and for each target we calculated primer efficiency. We set the limit of quantification (LOQ) at the last point of the standard line falling in the linear range. We assigned samples below the LOQ the log gene copies/g feces defined as LOQ/2 and defined the sample as below detection limit. We created a summary variable for all assessed virulence and toxin genes (VTGs) of pathogenic E. coli (eaeA, LT, ST, stx1, stx2) and for the sum of all pathogens (C. difficile, C. perfringens, EHEC stx1, EHEC stx2, EPEC, ETEC LT, ETEC ST, Samonella spp. and S. aureus). Gene copy number of C. difficile was divided by 11 to correct for multiple 16S rDNA genes.
We performed multivariate redundancy analyses (RDAs) on the infant and maternal gut microbiota composition as assessed by 16S rRNA gene sequencing in Canoco version 5.11 using default settings of the analysis type "Constrained" [10]. For the RDA at baseline, we used relative abundance values of OTUs from infant and maternal fecal samples, respectively, as response data, and maternal secretor status as explanatory variable. To investigate the interaction between secretor status and intervention groups in the RDAs, we used relative abundance values of OTUs from infant fecal samples at baseline, week 3 and 4 months, respectively, as response data, and the six intervention-secretor-status groups as explanatory variable. For visualization purposes families (and not OTUs) are plotted as supplementary variables. RDA calculates p values by permutating (Monte Carlo) the sample status.
We tested for between group differences at baseline (secretor vs. non-secretor mothers and infants of secretor mothers vs. infants of non-secretor mothers, respectively) in alpha diversity and abundance of the taxa of primary interest (Bacteroidetes, Bifidobacterium, Clostridiales, Enterobacteriaceae and Lactobacillus) and of all taxa, using Wilcoxon rank-sum tests with FDR correction for multiple testing. We tested for differences in alpha diversity (PD whole tree, and for Bifidobacterium diversity analysis also Shannon index), phylogenetic distance (weighted UniFrac within individuals between baseline and 3 weeks sample, and between baseline and 4 months) and abundance of the taxa of primary interest and of all taxa, at baseline, 3 weeks and 4 months between the intervention-secretor-status groups (control group infants of secretor mothers (control-S), control group infants of non-secretor mothers (control-NS), Fe group infants of secretor mothers (Fe-S), Fe group infants of non-secretor mothers (Fe-NS), FeGOS group infants of secretor mothers (FeGOS-S), and FeGOS group infants of non-secretor mothers (FeGOS-NS)), using Kruskal-Wallis tests with Dunn's post hoc test to adjust for multiple comparisons. For longitudinal analysis, change of taxon relative abundance over time, 2log ratios were calculated, in which the relative abundance of a taxon at 4 months or 3 weeks was divided by the relative abundance of the same taxon at baseline. Ratios were compared among groups by Wilcoxon rank-sum tests with FDR correction for multiple testing.
We analyzed demographic and anthropometric data, abundance of HMOs, enteropathogen abundances as assessed by qPCR, fecal calprotectin, intestinal fatty acid binding protein (I-FABP), hemoglobin (Hb), plasma ferritin (PF), soluble transferrin receptor (sTfR), body iron stores (BIS), C reactive protein (CRP) and alpha-glycoprotein (AGP) using the R statistical programming environment (R. 3.4.1 software; R Core Team). We run descriptive statistics for all variables and assessed normality by testing the distribution of continuous variables against a normal distribution using the Shapiro-Wilk W test. We transformed the variables with the use of log(x), sqrt(x), or −1 divided by x to correct positive skewness or x 2 , x 3 , or antilog(x) to correct negative skewness until achieving W > 0.97 before proceeding with further data analyses if departing significantly from normality. We show normally distributed data as mean ± SD and non-normally distributed data as median (IQR).
We tested for changes over time of lactation in HMO concentration using Wilcoxon signed-rank tests. We tested for differences in HMO concentrations between secretor and non-secretor mothers using t-tests or Wilcoxon rank-sum tests. Between group differences at baseline (secretor vs. non-secretor mothers and infants of secretor mothers vs. infants of non-secretor mothers, respectively) of enteropathogen abundances as assessed by qPCR, fecal calprotectin, I-FABP, Hb, PF, sTfR, BIS, CRP, AGP and infant anthropometrics were tested using t-tests or Wilcoxon rank-sum tests. Using linear regression models we investigated if maternal secretor status predicts abundance of enteropathogens, infant anthropometrics, hematological and inflammation status, fecal calprotectin and plasma I-FABP. We investigated the interaction between maternal secretor status (secretor and non-secretor) and intervention groups (control, Fe and FeGOS) on Bifidobacterium abundance, enteropathogen abundance, fecal calprotectin, I-FABP, anthropometrics, Hb, PF, sTfR and BIS, by fitting linear mixed-effect models using CRAN package lme. We defined the fixed effects on the variance as time, secretor status, intervention group, time-by-secretor-status, time-by-intervention-group, secretor-status-by-intervention-group and time-by secretor-status-by-intervention-group, and the random structure as the subject to control between-subject differences. We performed model diagnostics by visually inspecting the Tukey-Anscombe plot and the residual plot to verify the assumptions of homoscedasticity and normality of residuals. Post hoc analyses were performed to investigate the effect of time within group (intervention-secretor-status groups: control-S, control-NS, Fe-S, Fe-NS, FeGOS-S, FeGOS-NS) by applying the mvt method. We assessed infant morbidity as longitudinal prevalence, calculated as weeks with illness (RTI or diarrhea and/or mucus in the stool) from baseline to 4 month, divided by the total weeks of observation. We calculated longitudinal prevalence ratio (LPR) and 95 % CI using Poisson regression with robust standard errors [11]. We calculated p values as described by Altman DG and Bland JM [12].
In all analyses, p values < 0.05 were considered statistically significant. Table S3. Absolute and relative concentration of human milk oligosaccharides (HMOs) in breast milk samples of Kenyan mothers (n = 75), all and by secretor status. Relative concentration given as % of total HMOs (%total). Between group differences (secretor vs. non-secretor) were tested using Wilcoxon rank-sum tests; *p < 0.05, °p < 0.10. Table S4. Absolute and relative concentration of human milk oligosaccharides (HMOs) in breast milk samples of Kenyan mothers (n = 16) at two time points of lactation.   Figure S1. Gut microbiota composition as assessed by 16S rDNA sequencing among Kenyan mothers at baseline. A) All mothers, B) secretor mothers and C) non-secretor mothers. The fraction of 16S rDNA reads (%) attributed to specific taxonomic levels is given below the taxon name.  Figure S2. Gut microbiota composition as assessed by 16S rDNA sequencing among Kenyan infants at baseline. A) All infants, B) infants of secretor mothers and C) infants of non-secretor mothers. The fraction of 16S rDNA reads (%) attributed to specific taxonomic levels is given below the taxon name.