Microbiome clusters

We performed community typing analysis based on Dirichlet Multinomial Mixtures (DMM) modeling²⁸ to identify clusters of similar bacterial community structure amongst our samples. Based on their microbiota composition, the infant fecal samples clustered in four groups (Figure 1–2 and eFigure 1). Gut bacterial species richness (Figure 1b), alpha- (Figure 1c) and beta-diversity (Figure 1a) and taxonomic composition (Figure 2) differed between clusters, reflecting broad community differences. Clusters 1 and 4 comprised microbial communities reflecting the well-described effect of temporal maturation during the first year of life; with cluster 1 comprising only three-month (3 M) samples and cluster 4 comprising almost exclusively twelve-month (12 M) samples. Clusters 2 and 3 comprised a mixture of 3 M and 12 M samples. Compared to the other three clusters, cluster 1 showed a higher proportion of exclusive breastfeeding. Cluster 3 included a higher proportion of mothers receiving intrapartum antibiotics, infants born by C-section, and formula feeding (Figure 1). However, there was no difference in maternal ASB consumption between clusters, suggesting that this exposure did not influence the compositional differences that drove cluster classification (Figure 1f). In addition, the clusters did not differ in terms of maternal sugar intake, gestational diabetes, age, ethnicity, education, maternal gestational antibiotics, study site, infant antibiotics, or infant or mother secretor status.

Figure 1.

Discrepancies in covariate distribution, alpha- and beta-diversity between microbiota clusters

(A) Principal component analysis (PCoA) ordinations of variation in beta-diversity of infant gut bacterial communities based on Bray-Curtis dissimilarities among samples. Ellipses represent 95% confidence intervals. (B-C) Box plots showing the alpha-diversity (richness and Shannon’s diversity) per DMM cluster. The central line denotes the median, the boxes cover the 25th and 75th percentiles, and the whiskers extend to the most extreme data point, which is no more than 1.5 times the length of the box away from the box. Points outside the whiskers represent outlier samples. Letters denoted significant differences (non-parametric Kruskal-Wallis test followed by post-hoc test of Dunn with FDR correction following Benjamini-Hochberg method; P < .05). (D-K) Variable distribution between clusters tested with non-parametric Kruskal-Wallis test followed by either a post-hoc generalized linear model (glm) with a binomial/logistic distribution (D-I) or (J-K) a post-hoc Dunn test with FDR correction following Benjamini-Hochberg method. Minuscule letters indicate statistical differences between clusters from post-hoc generalized linear model (glm) with a binomial/logistic distribution. “BF at 3 M” stands for “breastfeeding at three months” and “FF at 3 M” for “formula feeding at three months”. Aside from maternal ASB consumption (F), only the variables that showed a statistical difference in distribution between clusters are presented. No differences were found for maternal age, ethnicity, education, diabetes; study site, household pets, siblings, or introduction of solid foods at 3 or 6 months. Cluster 1 included 48 samples from 48 infants; cluster 2 included 59 samples from 49 infants; cluster 3 included 47 samples from 39 infants; and cluster 4 included 44 samples from 43 infants. See methods for definition of variables. For all clusters with the same letter, the difference between groups is not statistically significant. If two variables have different letters, they are significantly different.

Figure 2.

Differences in relative abundances of the dominant bacterial genera between clusters

(A-J) Relative abundance across DMM clusters of the ten most dominant bacterial genera and (K) of the 15 most dominant bacterial genera. Letters indicate significant differences between clusters (non-parametric Kruskal-Wallis test, post-hoc Dunn test with Benjamini-Hochberg FDR correction). For all clusters with the same letter, the difference between groups is not statistically significant. If two variables have different letters, they are significantly different. Cluster 1 contains only three months of age. Clusters 2 and 3 are composed of a mix three and twelve months of age, and Cluster 4 only 12 M (except two samples).

Relative association of ASB on microbial community structure

Envfit analysis (univariable models) identified thirteen variables as significant drivers of gut bacterial beta-diversity from which we selected eight non-redundant variables to build our models: infant age, maternal intrapartum antibiotics, maternal ethnicity, birth mode, breastfeeding status at three months, presence of older siblings, infant secretor status, and maternal ASB consumption (Figure 3a and eFigure 2). Considering the complete dataset, the significant predictors were infant age, maternal ethnicity, intrapartum antibiotics, and birth mode. The same four variables, plus breastfeeding status at 3 months, were tested in a PERMANOVA (multivariable model), altogether explaining 14.2% of community variance (Table 2). Maternal ASB consumption was a significant predictor of infant gut bacterial composition only in the multivariable model (R² = 0.7%; Table 2). Birth mode (vaginal vs. C-section), intrapartum antibiotics, and breastfeeding status at three months had also a significant influence on community composition (respectively R² = 0.8%, 1.7%, and 1.9%), but to a lesser extent than infant age (R² = 7.3%) and mother’s ethnicity (R² = 2.5%; Table 2).

Figure 3.

Drivers of gut bacterial beta-diversity and indicator taxa associated with maternal consumption of ASB differ between clusters

Univariate models showing significance and explained variance of 10 variables on bacterial community structure across all data and each cluster subset. Horizontal bars show the amount of variance (R²) explained by each covariate in the model as determined by envfit. Asterisk denotes the significant covariates in each data subset (P < .05). All 32 variables considered in this study are shown in eFigure 2. In this figure, ASB represents artificially sweetened beverages and BF at 3 M represents infant’s breastfeeding status at three months (see methodology). (B) 14 bacterial taxa identified as significant features associated with maternal consumption of ASB by DESeq2.

Table 2.

Maternal consumption of ASB during pregnancy is associated with bacterial community assembly during the first year of life. Permutational Analysis of Variance (PERMANOVA) of gut bacterial community composition (Bray-Curtis dissimilarities) testing associations with different explanatory variables. The model on the complete dataset (ALL) accounts for repeated measures. The set of variables to be tested was chosen based on results from univariate envfit models

Variables	All (R2%)	Cluster 1 (R2%)	Cluster 2 (R2%)	Cluster 3 (R2%)	Cluster 4 (R2%)
Infant age (3 M vs. 12 M)	7.3***	8.5*	4.1***	8.0***	3.9**
Ethnicity	2.5***	NS	NS	NS	NS
Breastfeeding at 3 M	1.9***	5.1**	5.0***	6.0*	6.4**
Maternal Intrapartum Abx	1.7***	NS	NS	NS	NS
Birth mode	0.8**	NS	NS	NS	NS
Older siblings	NS	NS	NS	NS	NS
Infant secretor status	NS	NS	NS	NS	NS
Maternal ASB	0.7*	NS	3.2**	NS	NS
Total R2 (%)	15.1	13.6	9.1	14.0	10.3

^NSP > 0.05, * P < 0.05, ** P < 0.01, *** P < 0.001

Next, we repeated the beta-diversity analyses separately within each of the four clusters. Envfit univariable models identified distinct drivers for each cluster (Figure 3a). Interestingly, the drivers of beta-diversity in cluster 1 (only 3 M samples) were mainly maternal factors (i.e. birth mode, mother’s ethnicity, intrapartum antibiotics) whereas the drivers of cluster 4 (mostly 12 M) were infant factors (infant’s secretor status, breastfeeding at three months, and infant age (Figure 3a). Cluster 2 was the only cluster in which maternal ASB consumption was associated with beta-diversity (R² = 3.2%), and this association was confirmed by the univariable (Figure 3a, eFigure 2) and multivariable (Table 2) analyses.

We tested for associations of specific bacterial features in the infant gut with maternal ASB consumption. In the complete dataset, we identified two ASVs associated with maternal consumption of ASB, one species being depleted (Bacteroides sp. ASV45, log2 fold change = −27.2 and another species enriched (Prevotella copri ASV42, 24.2) among infants exposed to high maternal ASB intake (Figure 3b). Repeating this test within each cluster, we identified 15 additional ASVs enriched or depleted. For cluster 2, one ASV was enriched (ASV19, Akkermansia municiphila, 24.9) and four depleted (Bacteroides ovatus ASV27, −25.9; Parabacteroides sp. ASV83, −25.2; Bacteroides sp. ASV45, −24.9; Bacteroides sp. ASV25, −10.7) with maternal ASB consumption (Figure 3b). All adjusted p-values were below 0.001 and corrected with Benjamini-Hochberg for FDR.

Association between ASB exposure and urine metabolites

Using the software MetaboAnalyst,²⁹ we identified twenty metabolites that varied across clusters (eFigure 3). Since functional features in the metabolome are reflective of microbial metabolism redundancy, and less susceptible to interindividual variability typical of taxonomic datasets, we tested for a significant effect of maternal ASB consumption on the complete dataset at 3M and 12M old. Two urine metabolites, spermidine (log2 fold change = 2.27, p = .01) and succinate (log2 fold change = 1.77, p = .001), were significantly higher in 3M old infants exposed to ASB (Figure 4).

Figure 4.

Association between maternal ASB consumption, urine metabolite concentration in 3 M-old infants, and infant BMI at 12 months

Significant changes were detected in the normalized concentrations of (A) spermidine (log2 fold change = 2.27, p = .01) and (B) succinate (log2 fold change = 1.77, p = .001) in infant urine metabolomics by (C) volcano plot (combination of fold change and t-tests). (D) Linear regression of BMI at 12 months against normalized succinate concentration. Metabolite concentrations were normalized by median, log transformed and transformed with pareto scaling (mean-centered and divided by the square root of the standard deviation of each variable) with the software MetaboAnalyst.

Association of ASB, gut microbiome and urine metabolites with infant BMI at one-year-old

Finally, using a multivariable linear model on the complete dataset, we tested the association of maternal ASB consumption, microbial community composition and the two metabolites linked to high ASB consumption with infant BMI z-score at one year of age. In this cohort, infant birth weight was not correlated with BMI at one-year-old (p = .28). However, our multivariable linear model confirmed that daily maternal ASB consumption is associated with higher infant BMI (ß-estimate = 0.42, 95%CI 0.03:0.80, P = .037; Table 3), and showed that BMI was associated with the microbiome composition at 12 months (PCoA1 axis; ß-estimate = −0.71, 95%CI −1.40:-0.01, P = .048; Table 3) but not at three months (not shown). These results suggest that features of PCoA1 (i.e. lower relative abundance of Bacteroidetes and Faecalibacterium, and higher relative abundance of Escherichia, Klebsiella, Bifidobacterium, Haemophilus, Clostridium, and Veillonella; eFigure 4) are inversely associated with infant BMI. Notably, succinate was positively associated with BMI (ß-estimate = 0.45, 95%CI 0.15:0.76, P = .004; Table 3, Figure 4d) but not spermidine (P = .49), revealing that this association may also involve microbial-derived metabolites. We then performed a mediation analysis which showed that urine succinate mediated 29% of the observed association between ASB and BMI at one-year-old (mediation effect of ß-estimate = 0.29, 95%CI 0.03:2.19, P = .04).

Table 3.

Maternal consumption of ASB during pregnancy and urine succinate are associated with higher infant BMI at one-year-old. Linear models showing the explanatory power of maternal ASB consumption and urine succinate at 3 M old on infant BMI z-score at one year old, as well as the two main axes of ordination of bacterial community structure (beta-diversity) on samples acquired at 3 M and 12 M old. The full models are: [1]BMI at 1y ~ Succinate; [2]BMI at 1y ~ ASB + PCoA1 + PCoA2. Microbial variables were transformed (squared root and order quantile normalized respectively) to achieve normality. Here we present only the best model for 12 months fitted by stepwise selection by Akaike information criterion because we detected no association between BMI at one year old and microbiota composition at 3 M old

Variables	Infant BMI z-score at 1 year				Total adj. R2
	ß-est.	95% CI	P-value	R2
Model 1 Urine succinate (3 months)	0.45	[0.15,0.76]	0.004	7.4%	7.4%
Model 2 Maternal ASB (daily vs. no consumption)	0.42	[0.03,0.81]	0.037	4.1%	8.1%
PCoA axis 1	−0.71	[−1.40, −0.01]	0.048	3.9%
PCoA axis 2	NS	NS	NS	NS

Study design and population

We used data and samples collected through the CHILD Cohort Study,^50,51 a Canadian general population birth cohort (3621 families recruited across four provinces) including singleton pregnancies (>35 weeks gestational age with no congenital abnormalities) enrolled from 2008 to 2012. From this cohort, we completed a case-control study by selecting 100 infants divided equally between mothers that reported little or no ASB consumption (less than one per month) or high ASB consumption (one or more per day) during pregnancy. The groups were balanced for six potential confounding factors known to influence the gut microbiome: infant sex, birth mode, breastfeeding at three and 12 months, maternal BMI, and antibiotic use in infants before 12 months (antibiotics before three months old was an exclusion criterion; Table 1). To characterize the gut microbiome, stool samples were acquired at three and 12 months of age for a total of 200 samples. This study was approved by the University of Calgary Conjoint Health Research Ethics Board (CHREB) and ethics committees at the Hospital for Sick Children, and the Universities of Manitoba, Alberta, and British Columbia. Written informed consent was obtained from mothers during enrollment to the CHILD Study.

Maternal diet in pregnancy

Maternal dietary assessment in pregnancy has previously been described.³ Briefly, a food frequency questionnaire (FFQ) was completed during the second or third trimester and ASB consumption was evaluated using reports of “diet soft drinks or pop” (i.e. soda) (serving = 12 oz/one can) and “artificial sweetener added to tea or coffee” (serving = 1 packet). Other dietary variables included: sugar-sweetened beverages, Healthy Eating Index (HEI) total score (see eMethods) added sugar and total energy intake.

Infant BMI

BMI was measured by CHILD staff to the nearest 0.1 kg around one year of age (mean = 12.0 months ± 0.8 [sd]) and height to the nearest 0.1 cm. Age- and sex-specific BMI-for-age z-scores were calculated following the World Health Organization reference.⁵²

Other variables

The following variables were considered in univariable analyses (see eMethods): (1) infant’s sex, age at sample collection, breastfeeding duration (BF duration; months), breastfeeding status at three months (BF at 3 M; yes or no), diet at three and six months (Diet at 3 M and Diet at 6 M; both defined in 8 categories allocated based on the presence in the infant’s diet of breastfeeding, formula, and solids), solids at three and six months (Solids at 3 M and Solids at 6 M), formula feeding at three months (FF at 3 M), number of antibiotic treatments received from six to twelve months (Child 6–12 abx), and secretor status (determined from the single nucleotide polymorphism rs601338 in the FUT2 gene); (2) mother’s gestational diabetes, age, ethnicity, education, oral antibiotics received during gestation (Mother gestational abx), intrapartum antibiotics (Mother intrapartum abx), and secretor status (rs601338 SNP); (3) study site, presence of cats, dogs, and older siblings in the house.

Fecal samples DNA extraction and sequencing

After collection, fecal samples were frozen and stored at −80ºC. We extracted gut microbial DNA from fecal samples using the DNeasy PowerSoil kit (QIAGEN) according to the manufacturer’s instructions and amplified the V4 region of the 16S rRNA gene to generate ready-to-pool dual-indexed amplicon libraries as described previously⁵³ (see eMethods). Using the DADA2⁵⁴ pipeline, the final dataset contained 4,553,000 quality sequences, a mean (range) of 6,509 (22,995–68,265) sequences per sample identified as 954 unique bacterial Amplicon Sequence Variants (ASVs). Samples contained a mean of 40 (10–95) unique ASVs per samples.

Urine untargeted metabolomics

We have applied an untargeted quantitative metabolomics approach to analyze the samples using a combination of direct injection mass spectrometry with a reverse-phase LC-MS/MS custom assay. This custom assay, in combination with an ABSciex 4000 QTrap (Applied Biosystems/MDS Sciex) mass spectrometer, can be used for the identification and quantification of up to 150 different endogenous metabolites including amino acids, acylcarnitines, biogenic amines & derivatives, uremic toxins, glycerophospholipids, sphingolipids and sugars.^55,56 The method combines the derivatization and extraction of analytes, and the selective mass-spectrometric detection using multiple reaction monitoring (MRM) pairs. Isotope-labeled internal standards and other internal standards are used for metabolite quantification. The custom assay contains a 96 deep-well plate with a filter plate attached with sealing tape, and reagents and solvents used to prepare the plate assay. First 14 wells were used for one blank, three zero samples, seven standards and three quality control samples. For all metabolites, samples were thawed on ice and were vortexed and centrifuged at 13,000x g. 10 µL of each sample was loaded onto the center of the filter on the upper 96-well plate and dried in a stream of nitrogen. Subsequently, phenyl-isothiocyanate was added for derivatization. After incubation, the filter spots were dried again using an evaporator. Extraction of the metabolites was then achieved by adding 300 µL of extraction solvent. The extracts were obtained by centrifugation into the lower 96-deep well plate, followed by a dilution step with MS running solvent. Mass spectrometric analysis was performed on an ABSciex 4000 Qtrap® tandem mass spectrometry instrument (Applied Biosystems/MDS Analytical Technologies, Foster City, CA) equipped with an Agilent 1260 series UHPLC system (Agilent Technologies, Palo Alto, CA). The samples were delivered to the mass spectrometer by an LC method followed by a direct injection (DI) method.

Statistical analysis

We used Dirichlet Multinomial Mixtures (DMM) modeling²⁸ on 16S rRNA gene sequencing data to identify clusters of similar bacterial community structure amongst our samples (a technique known as community typing analysis, increasingly used in human microbiome studies^10,57–59). This technique is increasingly employed in microbiome studies for three reasons: (1) identification of unique microbial clusters is unsupervised; (2) cluster size depends on metacommunity variability; and (3) adequate explicit probabilistic model penalizes model complexity to optimize cluster number. The lowest Laplace approximation grouped our samples in four unique clusters (Figure 1–2 and eFigure 1).

The distribution of variables as well as the variation in bacterial richness (Chao 1), alpha-diversity (Shannon index), and community evenness (Shannon index/log_n(species richness)) across the DMM clusters were examined by non-parametric Kruskal-Wallis tests followed by post-hoc Dunn tests or generalized linear models (glm) with a binomial/logistic distribution. To explore the changes in taxonomical community structure at a fine scale, we tested for significant differences in the relative abundance of the 10 most dominant bacterial genera across clusters using non-parametric Kruskal-Wallis tests followed by post-hoc Dunn tests with Benjamin-Holmes False Discovery Rate (FDR) correction. To account for potential heteroskedasticity in bacterial community dispersion between groups and avoid the loss of information through rarefaction,⁶⁰ we performed a variance stabilizing transformation^60,61 prior to any statistical tests on beta-diversity. To select variables that could be drivers of infant gut bacterial community structure, we tested for correlations between our variables and community scores on the Principal Component Analysis (PCoA) ordination axes in univariable models (envfit function of vegan⁶²). The relative influence of the significant drivers of gut bacterial community structure was then assessed statistically in multivariate models using a Permutational Multivariate Analysis Of Variance (PERMANOVA based on Bray-Curtis dissimilarities; adonis function of vegan⁶²) with 999 permutations and visualized using PCoAs. We used DESeq2 (with Benjamini-Hochberg False Discovery Rate (FDR) correction) to test for differentially abundant bacterial taxa according to maternal ASB consumption on the 100 most relatively abundant bacterial taxa to limit spurious significance driven by very rare ASVs. Finally, we used linear models on the 3 M and 12 M-old samples to test for the influence of maternal ASB consumption, urine metabolites, and microbial ordination axes (PCoA1 and PCoA2) on infant BMI z-score. We performed the mediation analysis using “Succinate” as mediator and “ASB” as mediated with the mediation package. All analyses and graphs were computed in R version 3.6.1 (R Development Core Team; http://www.R-project.org) and MetaboAnalyst²⁹ for the metabolomics.