Examining the relationship between maternal body size, gestational glucose tolerance status, mode of delivery and ethnicity on mother’s milk microbiota at three months post-partum

Background: Few studies have examined how maternal body mass index (BMI), mode of delivery and ethnicity affect the microbial composition of human milk and none have examined associations with maternal metabolic status. Given the high prevalence of maternal adiposity and impaired glucose metabolism, and the importance of human milk in the colonization of the infant gut, we systematically investigated the associations between these maternal factors and milk microbial composition and functionality. Methods: Women ≥20 years were recruited during pregnancy and milk samples were collected at 3 months post-partum (NCT01405547). Demographic data, weight, height, and a 3-hour oral glucose tolerance test were conducted at 30 (95% CI: 25-33) weeks gestation. Metagenomic DNA extraction and 16S ribosomal RNA gene sequencing of the V4 hypervariable region (Illumina MiSeq) was carried out on 113 milk samples. Results: Multivariable linear regression analyses demonstrated no significant associations between maternal characteristics (maternal BMI [pre-pregnancy, 3 months post-partum], glucose tolerance, mode of delivery and ethnicity) and microbiota alpha-diversity; however, pre-pregnancy BMI was associated with human milk beta-diversity (Bray-Curtis p=0.040). Women with a pre-pregnancy BMI >30 kg/m2 (obese) had a greater incidence of Bacteroidetes (incidence rate ratio [IRR]: 3.70 [95% CI: 1.61-8.48]) and a reduced incidence of Proteobacteria (0.62 [0.43-0.90]), compared to overweight women (BMI 25.0-29.9 kg/m2) as assessed by multivariable Poisson regression. Increased incidence of Gemella was observed among overweight (versus healthy) mothers with gestational diabetes (5.96 [1.85-19.21]) and obese (versus healthy) mothers with impaired glucose tolerance (4.04 [1.63-10.01]). An increased incidence of Brevundimonas (16.70 [5.99-46.57]) was found in the milk of women who underwent an unscheduled C-section versus vaginal delivery. Lastly, functional gene inference demonstrated that obesity was associated Separate statistical models built using pre-pregnancy and 3-month post-partum BMI collinearity. interaction term BMI

Further research is warranted to determine whether this variability in the milk microbiota impacts colonization of the infant gut.

Background
Breastfeeding is the recommended method of feeding for all infants irrespective of whether the country of origin is low, middle, or high-income (1). Mother's milk is a rich source of nutrients and bioactive components, such as the antimicrobial proteins lactoferrin and lysozyme, and contains an array of oligosaccharides which serve as a source of prebiotics (2)(3)(4)(5). It is now well accepted that mother's milk contains a rich supply of bacteria (~ 10 6 bacterial cells/mL), which are believed to play an important role in postnatal colonization of the infant's gastrointestinal tract (gut) (6)(7)(8)(9)(10)(11)(12)(13). Microbial composition of the gut, in turn, is associated with maturation of an infant's gut and immune system, and aberrant gut microbial compositions have been linked to a number of short-and long-term health outcomes including diarrhea, respiratory tract infection, asthma, inflammatory bowel disease, obesity, and metabolic syndrome (14)(15)(16)(17)(18).
There is a limited understanding of the stability of the human milk microbiota in the face of environmental influences. Despite the high prevalence and known impact on other milk constituents, no study we are aware of has examined the association between maternal glucose tolerance status, either gestational diabetes or impaired glucose tolerance, and the milk microbiota. Only a few studies to date have cross-sectionally examined other maternal factors, such as maternal body size and mode of delivery, on the microbiota composition of mature mother's milk (collected from 1-week to 6months post-partum) (19)(20)(21)(22)(23). Moreover, the findings from the few available studies are inconsistent likely due, in part, to the small number of study participants, and differing methods of milk collection and analysis (Supplementary Table 1). Of the limited studies conducted, maternal BMI and mode of delivery have been associated with the mature milk microbiota; however, many of these reports relied on small cohorts and thus their findings require replication in larger scale studies (19)(20)(21)(22)(23).
Further, most studies were unable to employ multivariable statistical modelling to adjust for multiple potential maternal factors of interest simultaneously. Lastly, none of these studies carried out functional inference analyses to determine if maternal body size and mode of delivery also perturb the milk microbiome's potential metabolic activities or functional capabilities. It is vital to determine if maternal metabolic (BMI, glucose tolerance status), obstetrical (mode of delivery), and demographic (ethnicity) factors impact the mother's milk microbiota and metagenome due to the role they play in the colonization of the infant gut and the importance of a healthy gut microbiome in human health.
Therefore, we set out to fill these knowledge gaps in the field. We collected a large cohort of mother's milk samples from women prospectively enrolled in a study to investigate the impact of metabolic abnormalities and maternal nutrition in pregnancy on human milk (NCT01405547). The objective of this current study was to investigate the associations between maternal pre-pregnancy BMI (healthy, overweight, or obese), 3-month post-partum BMI, maternal glucose tolerance status in late pregnancy (gestational diabetes mellitus [GDM], impaired glucose tolerance [IGT], or normoglycemia), mode of delivery (vaginal delivery, unscheduled Caesarean delivery [C-section], or scheduled C-section), and ethnicity (white, Asian, or other) on the microbial community composition and functional capabilities of mother's milk at three months post-partum. This research will help to elucidate the modulatory potential of mother's milk microbiota in the face of physiological perturbations.

Participant Description
Milk samples were collected at 3 ± 1-month post-partum (mean ± standard deviation) (n = 113). Fiftysix (49.6%) mothers fed their infants their own milk exclusively at the time of milk collection, and 61 (53.9%) samples were from a complete breast expression. The mean (± standard deviation) age of the mothers was 34.2 ± 4.2 years (Table 1) with a pre-pregnancy BMI (kg/m 2 ) of 24.3 ± 4.6, which is within a healthy BMI range (18.5-24.9 kg/m 2 ). A modified oral glucose tolerance test (OGTT) administered at 30 weeks' gestation (95% CI: 25-33 weeks) revealed that 24 (21.2%) women had GDM, 20 (17.7%) had IGT, and 69 (61.1%) had healthy glucose metabolism (normoglycemic). Among women with a healthy pre-pregnancy BMI, 12 and 16 had IGT and GDM, respectively. Among women with an overweight pre-pregnancy BMI, 5 had IGT and 5 had GDM (10 total). Lastly, among women with an obese pre-pregnancy BMI 3 had IGT and 3 had GDM (6 total; Supplementary Table 2). Sixtyfour (56.6%) mothers delivered their infants vaginally, compared to the 49 (43.4%) who underwent a C-section. Following 16S rRNA gene sequencing, rarefying, and filtering, the bacterial taxa from 109 milk samples were included in the analyses. The total counts per sample were rarefied to 20,000 reads prior to calculating diversity metrics. Overall microbial composition of mother's milk Twenty-six unique phyla-level and 292 unique genus-level taxa were identified ( Figs. 1 and 2).
We assessed whether taxa abundance was associated with maternal characteristics by using multivariable Poisson regressions and accounting for multiple comparisons (Tables 2, 3, 4, 5;   Supplementary Tables 5, 6). At least one maternal characteristic was associated with the differential abundance of Proteobacteria, Bacteroidetes, and Actinobacteria at the phylum level, and Staphylococcus, Veillonella, Gemella, Corynebacterium, and Brevundimonas at the genus level.
Association between maternal body size and the milk microbiota Pre-pregnancy BMI (i.e., healthy, overweight, obese) was found to be most consistently associated with differentially abundant taxa after controlling for relevant maternal characteristics. Mothers categorized as obese pre-pregnancy displayed a lower incidence of Proteobacteria (incidence rate Association between maternal glucose tolerance status and the milk microbiota When examining an interaction between BMI and maternal glucose tolerance status, Gemella showed an increased incidence among overweight (versus healthy) mothers with gestational diabetes (5.96 [1.85-19.21]). In addition, Gemella was increased in obese mothers with impaired glucose tolerance versus overweight (11 Association between mode of delivery and the milk microbiota Associations between mode of delivery and the differential abundance of select taxa at the phylum and genus level were found for both pre-pregnancy and post-partum BMI models ( Table 2)

Association between maternal ethnicity and the milk microbiota
Lastly, ethnicity (white, Asian, other) was associated with the differential abundance of Corynebacterium and Brevundimonas (Table 3). White mothers had a reduced incidence of both Association between maternal body size, glucose tolerance status, mode of delivery, ethnicity and functional gene expression of the milk microbiota We carried out functional inference analyses using Piphillin to assess if there were any differences in functional capabilities of the milk microbiota based on the maternal clinical data. In contrast with the bacterial taxonomic results, the relative abundance of the 20 top KEGG pathways across all milk samples was fairly consistent (Supplementary Fig. 3). We analyzed the association between maternal clinical data and KEGG ortholog beta-diversity as well as examined the association between maternal clinical parameters and differentially-expressed KEGG pathways (Supplementary Tables 7-9). No significant results were found when examining metadata and KEGG ortholog beta-diversity (Supplementary Table 7); however, one statistically significant differentially-expressed pathway was observed (Supplementary Table 8, Supplementary Fig. 4). BMI, specifically the obese sub-category, was shown to be significantly associated with enrichment of the KEGG category "Biosynthesis of secondary metabolites" (Coefficient 0.00024, p = 0.0079). (Supplementary Fig. 4).

Discussion
Our results suggest that maternal factors, and most consistently maternal pre-pregnancy BMI, are associated with the microbial composition of mother's milk. This is the first study to include maternal glucose tolerance status in the investigations of the association between maternal body size and the milk microbiota (Supplementary Table 1). Gestational diabetes, a well-known risk factor of maternal adiposity, is associated with a number of negative health outcomes, including increased risk of type 2 diabetes and metabolic syndrome in the mother, as well as large for gestational age, congenital malformations and hypoglycemia in the infant (25)(26)(27)(28). These changes in infant health may be partially mediated by microbes transmitted from mother to infant, from both the birthing process as well as during direct breastfeeding. For this reason, it is imperative to study the relationship between these factors in order to best mitigate these maternal and infant outcomes as well as distinguish the role of body size versus glucose tolerance status on the milk microbiota and, hence, the infant.
According to a recent systematic review and dose-response meta-analysis, the risk of GDM increases by 4% for every unit increase in BMI; thus, to investigate BMI without also adjusting for GDM could lead to erroneous results (28). Our results suggest that these associations extend beyond the gut microbiota and that alterations in glucose tolerance status may also perturb the composition (Tables 2, 4) of mother's milk microbiota. Previous studies have reported associations between both type 2 diabetes and insulin resistance and the gut microbiota composition (29). Potential mechanisms of action linking impaired glucose tolerance and the gut microbial community include lipopolysaccharide-triggered inflammation, impairment of GLP-1 and GLP-2 via bacterial production of short chain fatty acids, insulin resistance triggered by bacterial synthesis and absorption of branchchained amino acids, or the metabolism of bile acids by bacteria and their organ-specific effects (29).
It is unclear how impaired glucose tolerance may modulate the bacteria in human milk and the interplay present with maternal body size.
We did not find any differences in alpha diversity based on our maternal characteristics; however, we did find statistically significant differences in beta-diversity, with mother's milk microbiota separating, or non-randomly clustering, based on pre-pregnancy BMI even after adjustment for other covariates (Supplementary Table 4, Supplementary Figs. 1 and 2). The human gut microbiota has been reported to cluster as a function of body size, but this has not yet been reported for the human milk microbiota (30). Our results demonstrating an association between maternal body size and microbial composition is consistent with other smaller scale studies on human milk microbiotas (Supplementary Table 1).
Cabrera-Rubio et al. (2012) examined the association between maternal body size and the differential abundance of mother's milk genera in a study of healthy Finnish women (n = 18) (20). They reported an increase in Staphylococcus in mother's milk collected from obese women, which mirrors the findings in our study ( Table 4).
Mode of delivery was also associated with changes in mother's milk microbiota at both the phylum and genus levels ( Table 2, 4). For example, we observed greater differential abundance of Staphylococcus in mother's milk from women who underwent a (scheduled) C-section versus vaginal delivery (Table 4). Our results are similar to that reported by Cabrera-Rubio et al. (2012, 2016). These two small cross-sectional cohorts of healthy Finnish women (n = 18, 10) showed a non-statistically significant increase in Staphylococcus in milk observed among women who delivered their infant via a scheduled C-section versus a vaginal delivery (Table 4) (19,20). The proposed mechanism whereby mode of delivery alters the milk microbiota is via the infant oral cavity, which is colonized during either vaginal delivery or C-section; from here, retrograde inoculation of bacteria can occur from the infant's oral cavity into the mammary gland via the suckling process with direct breastfeeding (31)(32)(33)(34).
The results of our multi-ethnic cohort revealed associations between ethnicity and specific bacterial taxa. Ethnicity and/or geographic location have been shown to be factors in determining various microbiomes of the body including the gut, oral cavity, respiratory tract, skin, and urogenital tract (35). Ethnicity and geographic location typically come with an overlay of dietary variation, making the impact of each variable challenging to separate. Only a few studies to date have assessed associations between ethnicity and the milk microbiota; however, the ethnic/geographic groups differed from our study as they generally examined Europe, Africa and the United States, making it challenging to compare findings (Tables 3, 5; Supplementary Tables 5, 6) (12,21,36,37).
We used a functional inference approach to characterize the microbial genetic potential in human milk. In agreement with what has been reported for other human-associated microbiotas, the functional capacity of the human milk microbiota is more stable than its taxonomic composition (38).
We then assessed whether there were specific pathways that were associated with maternal characteristics and found that maternal BMI, specifically the obese sub-category, was significantly associated with an increase in the "Biosynthesis of secondary metabolites" KEGG pathway. Microbes produce secondary metabolites, which are small, bioactive molecules, not necessary for growth or development but are instead involved in microbe-host or microbe-microbe interactions (39,40). Indeed, many of the genes in the biosynthesis of secondary metabolites pathway encode for the biosynthesis of antibiotics (41). Increased production of these secondary metabolites could impact the overall microbial composition/function in these feeding infants. These potential alterations could represent a mechanism by which maternal BMI impacts infant health over both the short-and longterm and warrants future investigation.
Human milk is considered a low biomass sample and, for this reason, may be more affected by sample processing than higher biomass samples, such as stool. To address this concern, we used PCoA plots to visualize clustering, or lack thereof, of our milk samples and negative controls ( Supplementary Fig. 5). Our negative controls were seen to cluster away from the samples, suggesting that our results do not arise from technical contaminates, which was confirmed using Adonis analyses to statistically corroborate that our samples clustered away from the negative controls (Weighted UniFrac R 2 = 0.07, p = 0.0001; Bray-Curtis R 2 = 0.10, p = 0.0001, Supplementary   Fig. 5).
Strengths of the current study include its relatively large sample size, the diverse ethnicity of women included, clinical examination via an OGTT, and enrichment of the cohort with women of varying body size who had abnormal glucose tolerance status. These strengths allowed for a more fulsome investigation using multivariable statistics to determine how each maternal factor is independently associated with the milk microbiome. Limitations of the current study include a lack of disinfection of the mother's breast, peri-areolar region, and/or nipple prior to milk sampling, and single time-point sampling. The microbes identified in the milk from the present study likely include bacteria from the skin microbiota. Practically, however, mothers do not disinfect their breast prior to pumping and storing milk for their infant, nor do they disinfect prior to breastfeeding. Therefore, the mother's milk microbiota as collected in the current study is likely a more accurate depiction of what the infant would receive. Secondly, we could not adjust for all variables of interest to due sample size constraints, so we are likely missing additional determinants of both the milk microbiota and functional capabilities. Lastly, our study is limited by its cross-sectional design and thus we cannot assess how the milk microbiome changes over time. It is possible that the associations we identified between maternal factors and microbial composition in milk are not transitory and change across the course of lactation.

Conclusions
Our study found that mother's milk has a highly personalized microbiota with high inter-individual variability. Expressed mother's milk at both the phylum and genus levels appear to be related to maternal metabolic and obstetrical factors. Surprisingly, glucose tolerance status was significantly associated with fewer microbiota parameters than anticipated. Most consistently, maternal prepregnancy BMI, despite glucose tolerance status, was associated with the differential abundance of various taxa in mother's milk and potentially the production of bacterial secondary metabolites as well. To understand the clinical significance of these findings, future research should explore the impact of differences in the microbial composition of mother's milk on colonization of the infant gut microbiome and infant health.

Study participants and design
To address the research objectives of this study, we used maternal metabolic and obstetrical health data and bio-banked human milk samples available from a previously conducted prospective cohort study (ClinicalTrials.gov Identifier: NCT01405547); a detailed description of the study protocol has been previously published (42). Pregnant women (n = 216) were recruited from outpatient clinics at

Mount Sinai Hospital in Toronto, Canada and completed a 3-hour 100 g OGTT between March 2009
and July 2010. In total, 117 women donated a milk sample at 3 months post-partum, with 113 samples available for this study (Supplementary Fig. 6). Women were eligible for inclusion in the original study if they were ≥ 20 years of age and had an intention to breastfeed. Exclusion criteria included pre-existing diabetes diagnosis, current use of insulin, or completion of an OGTT prior to recruitment (43). By design, mothers were recruited from clinics which follow higher risk pregnancies with a greater risk of either GDM or IGT diagnosis. Written informed consent was obtained from all women and the study protocol was approved by the Mount Sinai Hospital Human Research Ethics Board (42).

Collection of Demographic, Anthropometric and Metabolic Data
During the first study visit, which occurred in late pregnancy (30 weeks [95% CI: 25-33 weeks]), demographic and anthropometric data were collected (e.g. age, ethnicity, weight, height); mothers were asked to recall their pre-pregnancy weight. All pregnant women in Canada are screened for GDM by way of a 50 g glucose challenge test (GCT). If the plasma glucose concentration at 1-hour postglucose load is ≥7.8 mmol/L, the patient is then referred for a diagnostic OGTT. Contrary to standard obstetrical practice, all women completed a 3-hour 100 g OGTT in the current study during their first study visit regardless of whether or not they completed a GCT. The OGTT involved having blood samples drawn at fasting, 30, 60, 90, 120-and 180-minutes post-glucose load. Women were then diagnosed with either GDM, IGT, or as normoglycemic based on the following glycemic thresholds: 1) GDM diagnosis = 2 or more of the following: fasting blood glucose ≥ 5.8 mmol/L, 1-hour blood glucose

Data analysis and statistics
The phyloseq package (1.25.2) in R (version 3.4.1) was used to analyze microbiota composition (49).
OTUs that only appeared once or twice (singletons and doubletons) were removed and all OTUs were rarefied to 20,000 reads prior to calculating relative abundances at different taxonomic levels, alpha diversities and beta diversities using phyloseq. Statistically significant differences between the alpha diversities (Chao1/Shannon indices determined in R) and maternal metabolic or obstetrical characteristics were determined using multivariable linear regression models (PROC MIXED) in SAS version 9.4. Independent variables included in the models were: maternal BMI (healthy = 18.5-24.9 kg/m 2 , overweight = 25-29.9 kg/m 2 , obese = > 30 kg/m 2 ), maternal glucose tolerance status (GDM, IGT, normoglycemic), mode of delivery (vaginal, unscheduled C-section, scheduled C-section), DNA extraction batch, and PCR sequencing batch.
Separate statistical models were built using pre-pregnancy and 3-month post-partum BMI as covariates, due to concerns about collinearity. An interaction term between BMI and maternal glucose tolerance status was also tested in each model and removed if it was non-significant. Due to our sample size and the number of covariates we wished to test, separate models for ethnicity (white, Asian, other) were run that adjusted for DNA extraction and PCR sequencing batches, but no other covariates. Multicollinearity was assessed between independent variables in all models, using a variance inflation cut-off of > 5. The significance level was set at p < 0.05. Of note, 6 mothers with a pre-pregnancy BMI between 18.0-18.4 kg/m 2 were placed in the "healthy BMI" group for all analyses.
Beta diversities and principal coordinate analysis (PCoA) were also ascertained in phyloseq and statistical significance based on maternal characteristics was determined using the adonis function in vegan (version 2.5-3) (24). Adonis assesses the amount of variation explained by each metadata variable, such as maternal BMI or glucose tolerance status; all variables were run individually and together in adonis to adjust for one another. The interaction term between BMI and glucose tolerance status was also tested and removed if non-significant. Four patient's samples were missing the postpartum BMI data and thus we used their pre-pregnancy BMI for the post-partum analyses. Again, the significance level was set at p < 0.05.
Multivariable Poisson regression models (PROC GENMOD) were run in SAS version 9.4 to assess differential abundance at the phylum and genus levels based on maternal characteristics. The Benjamini-Yekutieli cut point approach was used to account for multiple testing. A p ≤ 0.022 at the phylum level (5 tests) and p ≤ 0.017 (10 tests) at the genus level were considered statistically significant for the overall group effect. If the overall group-adjusted p-value was significant, pairwise comparisons were conducted and a pairwise p < 0.05 was considered statistically significant.
Piphillin: Functional analysis of human milk microbiota Piphillin, a metagenomics inference tool, was used to infer functional capabilities in milk samples (https://piphillin.secondgenome.com/) (50). In this study, the Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) was used as a reference database to retrieve gene copy numbers and create a gene feature table from the 16S rRNA sequence data. Statistically significant associations between maternal characteristics and KEGG pathways where assessed in two ways: first, examining the association between maternal characteristics and KEGG orthologs beta-diversity using Adonis in R (p < 0.05), and second, investigating metadata associated with differentially-expressed functional pathways using MaAsLin2 in R (p < 0.1).

Availability of data and materials
The dataset supporting the conclusions of this article is available in the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA516669.

Competing Interests
None to declare.     Separate Poisson regression models were run for pre-pregnancy BMI and 3-month post-partum BMI, while adjusting for maternal glucose tolerance status, mode of delivery, DNA extraction batch, and PCR sequencing batch. All models were run testing for interaction terms between BMI and maternal glucose tolerance status. Statistically significant findings shown only (group effect: p≤0.022 for phylum, p≤0.017 for genus; pairwise comparison: p<0.05). Abbreviations: confidence interval, CI; incidence rate ratio, IRR. Ethnicity was investigated for all taxa and models were adjusted for DNA extraction and PCR sequencing batch effects. Statistically significant findings shown only (group effect: p≤0.022 for phylum, p≤0.017 for genus; pairwise comparison: p<0.05). Abbreviations: confidence interval, CI; incidence rate ratio, IRR. Figure 1 Microbial relative abundance in mother's milk at the phylum level. The relative abundances of bacterial phyla in collected mother's milk samples (n=109) are visualized using bar plots.

Figures
For simplicity, only the most abundant 5 phyla are displayed with other phyla merged into the Other category.

Figure 2
Microbial relative abundance in mother's milk at the genus level. The relative abundances of bacterial genera in collected mother's milk samples (n=109) are visualized using bar plots.
For simplicity, only the most abundant 10 genera are displayed with other genera merged into the Other category.