Population characteristics

In our study population, 74 children were VD, and 46 children were born by CS. Of those, 36 (78%) were born by planned, and 10 (22%) by emergency CS. There were two cases of pre-/intrapartum antibiotics, one in each delivery mode group, indicated for maternal fever, both of which were included in the analyses. All but three children were born in the hospital. Baseline characteristics and cumulative disease parameters over the first year of life of all children, stratified by mode of delivery, are shown in Table 1. Clinical variables were evenly distributed over both groups with the exception of gestational age (two-sample t test, p = 0.003), duration of ruptured membranes (p = 0.019), hospital stay duration (Wilcoxon test, p < 0.001), and total duration of breastfeeding in the first year of life (p = 0.014), all being intrinsically related to delivery mode. The number of children receiving exclusive formula feeding did not differ between the two groups. Only 36 children (30%) received antibiotics during their first year of life, some receiving multiple courses (56 courses in total for all children), mostly (80%) indicated for RIs.

Table 1 Baseline characteristics Full size table

Microbiota composition and mode of delivery

Of the 1243 fecal samples available from our 120 participants and their mothers, 1139 (92%) passed the quality criteria for further analysis following DNA extraction and 16S rRNA-based sequencing of the V4 hypervariable region (Supplementary Fig. 1), representing 70,886,595 high-quality reads in total. The Good’s coverage of the included children’s samples was high, with a minimum of 99.56% (median 99.96%). The raw Operational Taxonomical Unit (OTU)-table contained 623 OTUs distributed over seven bacterial phyla, with the Firmicutes generally being the most prominent phylum.

We observed that the infants’ overall microbial community composition developed slowly toward an adult-like profile (mothers’), though had not yet reached full maturation to an adult-like composition at 12 months of age (Fig. 1a). This was illustrated by a steady decrease in Bray–Curtis (BC) dissimilarity index between infants’ and mothers’ samples over time, with a median index of 0.999 directly postpartum, which reached 0.739 at the end of the first year (Supplementary Table 1). We observed clear differences in the early development of the overall community composition between VD and CS children, with a maximum effect of delivery mode at 1 week of life (permutational multivariate analysis of variance [PERMANOVA] test, R2 = 0.142, adjusted p-value 0.003, Benjamini–Hochberg method11) and significant differences until the age of 2 months (R2 = 0.021, adjusted p-value 0.055), after which these differences gradually disappeared (Fig. 2). To rule out this finding was due to the indirect exposure of CS children to maternal antibiotics through breastfeeding, we repeated this analysis post hoc on a subset of 11 VD and 11 CS children who received exclusive formula feeding. We observed similar (at some timepoints even bigger) effects (R2) of delivery mode on overall community composition until 2 months of life within this subset. The association between delivery mode and composition was still significant at 1 week (Fig. 2; Supplementary Table 2; R2 = 0.215, adjusted p-value 0.008) and 2 weeks of life (R2 = 0.152, adjusted p-value 0.044) in the exclusively formula fed children. For the whole study population, we also found that the microbial community in VD children was more stable when compared with CS children until 2 months of life (Fig. 1b). The BC dissimilarity between consecutive samples until 4 months was higher in CS, compared with VD children (Mann–Whitney test, p < 0.001, p = 0.001, and p = 0.042, for the intervals 1–2 weeks, 2 weeks–1 month and 1–2 months, respectively). In general, alpha diversity increased directly after birth, and again after 4 months of life, coinciding with the age that solid food was introduced to the children’s diet (median = 128 days, IQR = 119.8–164.2 days). There were no significant differences in alpha diversity between the two delivery mode groups at any timepoint (Supplementary Fig. 2). When testing the effect of delivery mode on alpha diversity longitudinally with a mixed effect model, no significant effect was found (ANOVA, p = 0.511), and neither did feeding type have an effect in this model (p = 0.652).

Fig. 1 Overall gut microbiota community composition development and stability. a Nonmetric multidimensional scaling (nMDS) plot, based on Bray–Curtis (BC) dissimilarity between samples, with data points and ellipses colored by timepoint. Children’s overall gut community composition developed toward a more adult-like pattern in the first year of life, becoming more similar to microbiota of adults (mothers’ samples, n = 87). b As measure of stability, we calculated BC dissimilarities between consecutive sample pairs belonging to an individual per time interval and plotted these at the end of each interval (t + 1). Loess lines were fitted over the data points per delivery mode group, and the gray areas represent the 0.95 confidence intervals. Stability was significantly lower in C-section born infants until 2 months of life. Source data are provided as a Source Data file Full size image

Fig. 2 NMDS plots of children’s samples per timepoint stratified according to mode of delivery. Nonmetric multidimensional scaling (nMDS) plots, based on Bray–Curtis (BC) dissimilarity between samples, visualizing the overall gut bacterial community composition stratified for mode of delivery, per timepoint. Each data point represents the microbial community composition of one sample. The ellipses represent the standard deviation of data points belonging to each birth mode group, with the center points of the ellipses calculated using the mean of the coordinates per group. The stress of the ordination, effect sizes (R2) calculated by multivariate permutational multivariate analysis of variance (PERMANOVA) tests and corresponding adjusted p-values (p.adj) are shown in the plots, and n represents the biologically independent samples per group Full size image

Covariates that were significantly associated with fecal microbiota composition over time, as tested with the adonis2 function12 (PERMANOVA test) were, besides mode of delivery (R2 = 0.013, adjusted p-value = 0.001): age (R2 = 0.034, adjusted p-value = 0.001), breastfeeding (R2 = 0.007, adjusted p-value = 0.001), daycare attendance (R2 = 0.006, adjusted p-value = 0.001), siblings < 5 years of age (R2 = 0.006, adjusted p-value = 0.001), pacifier use (R2 = 0.005, adjusted p-value = 0.003) and antibiotics in the 4 weeks prior to sampling (R2 = 0.003, adjusted p-value = 0.03). Pets in the household (R2 = 0.006, adjusted p-value = 0.061) and duration of hospital stay after birth (R2 = 0.002, adjusted p-value = 0.061) showed a trend toward being associated with fecal microbiota composition over time (Supplementary Fig. 3).

Fecal microbiota seeding from mother to infant

To study the existence of direct maternal fecal microbiota seeding during birth, and to assess the role of delivery mode herein, we studied the concordance of the microbiota composition of children’s fecal samples and their mother’s microbiota over time, in relation to the concordance with other mothers’ samples. Using a linear mixed model, we found that an infant’s fecal microbiota composition was significantly more similar to that of its own mother than to that of other mothers in VD children when studied over the entire first year of life (ANOVA, p = 0.025; Fig. 3), but not in CS children (p = 0.271). This difference between groups seemed independent of the intravenous antibiotics administered to the mothers in the CS group after cord clamping, as the overall fecal microbiota composition of CS and VD mothers themselves did not differ shortly after birth (PERMANOVA test, R2 = 0.013, p = 0.351; Supplementary Fig. 4).

Fig. 3 Comparison of overall composition between children and mothers (own vs. other). Children’s fecal microbiota were compared to the mothers’ fecal microbiota, based on BC dissimilarity and stratified according to mode of delivery. A significantly lower dissimilarity (more comparable microbiota) was observed between a child’s microbiota and its own mother vs. other mothers in children born vaginally throughout the first year of life, but not in children born by C-section (linear mixed models, ANOVA; p = 0.025 and p = 0.271, respectively). Boxplots with medians are shown; the lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles); the upper and lower whiskers extend from the hinge to the largest and smallest value no further than 1.5*IQR from the hinge; outliers are plotted individually. Pp = postpartum, d = day, m = month, n = number of mother-own-infant pair comparison per timepoint. Source data are provided as a Source Data file Full size image

Dynamics of microbiota development

The succession pattern of bacterial taxa in the VD children in our study population was consistent with the description of normal early-life gut microbiota development in previous studies (Fig. 4a)13,14. Facultative anaerobic genera, such as Escherichia, and Staphylococcus were highly abundant in the earliest samples, gradually making way for a predominance of the genus Bifidobacterium. Using smoothing spline analysis of variance (SS-ANOVA), we observed, among others, that Bifidobacterium was more abundant in VD than in CS children from day 1 until day 30, even when correcting for breastfeeding (adjusted p-value 0.003; Supplementary Table 3). Also, Escherichia was more abundant in VD compared with CS born children in the first 85 days of life. In contrast, in CS children we found, among others, higher abundances of Klebsiella from birth to day 139 and Enterococcus between 7 and 35 days (both adjusted p-values 0.003). The dynamics of differences in Bifidobacterium, Escherichia, Klebsiella, and Enterococcus over time are visualized in Fig. 4b, underlining the effect size and duration of differences.

Fig. 4 Mean relative abundance of most abundant OTUs. a Mean relative abundances of the 12 most abundant OTUs are depicted for all samples per timepoint, stratified by birth mode. Pp = postpartum, d = day, m = month. b Mean relative abundances of Bifidobacterium, Escherichia, Klebsiella, and Enterococcus over time. Loess lines were fitted over the data points per delivery mode group and the gray areas represent the 0.95 confidence intervals. Source data are provided as a Source Data file Full size image

To confirm that these results were not the consequence of indirect antibiotic exposure of infants through breast milk, we executed a post hoc SS-ANOVA analysis for the subset of exclusively formula fed children. We analyzed the top five most abundant taxa over the first 2 months of life, and again found an increased abundance of Bifidobacterium (days 5–44, adjusted p-value = 0.003) in VD infants, and an increased abundance of Klebsiella in the CS born infants (days 10–20, adjusted p-value 0.020). Also, Staphylococcus was found to be more abundant in the CS children from 0 to 6 days (adjusted p-value = 0.020).

We used mixed effect models to study the potential associations between delivery mode, age, feeding type, antibiotic use and hospital stay duration, and the five most abundant taxa. We found that the abundance of Bifidobacterium was associated with mode of delivery (ANOVA, p = 0.004), age (p < 0.001), and breastfeeding (p < 0.001). Surprisingly, breastfeeding did not compensate for the lack of Bifidobacterium in children born by CS: children born by CS and receiving breastfeeding had less Bifidobacterium present in their fecal samples than formula fed VD children (at 1 week of life, Wilcoxon test, p < 0.001, median relative abundance 0.016 and 45.1%, respectively). The triad CS birth, age, and formula feeding were also positively associated with the abundance of Enterococcus and Klebsiella. In addition, Klebsiella abundance was positively associated with having received antibiotics in the previous 4 weeks of life (p = 0.001). Escherichia abundance was only associated with vaginal delivery (p = 0.003) and age (p < 0.001). Staphylococcus abundance was associated with age and breastfeeding (both p < 0.001), but not with delivery mode. We did not find an association between duration of hospital stay after birth and any of these taxa.

Among the remaining mode of delivery-associated taxa observed (see Supplementary Table 3), we found to be of particular interest that Bacteroides spp., which are considered to be important regulators of intestinal immunity15, were more abundant in the VD compared with CS children in the first months of life.

Delivery mode-induced microbiota changes and infant health

Since delivery mode is reported to be associated with infant and childhood health, especially regarding respiratory illness16, we defined a secondary research question, namely whether gut microbiota development is associated with health outcome. Although it was not our aim to study differences in health outcomes between the delivery mode groups in our cohort, we did find a trend toward differences in infectious disease and treatment parameters, specifically parent-reported RI events and antibiotic courses over the first year of life (Table 1, chi-square test, p = 0.119 and p = 0.100, respectively). Exploring this further with a temporal post hoc analysis, we additionally found a trend toward a lower hazard ratio for antibiotic prescriptions in VD children in the first year of life (Cox proportional hazard model, HR = 0.606, p = 0.134).

Altogether, these results supported the validity of our secondary aim to investigate the potential role of delivery mode-induced gut microbiota changes on health. To test this, we studied the association between fecal microbiota composition at 1 week of life (where the maximum effect of mode of delivery on microbiota composition was observed; Fig. 2), and all commonly observed health parameters in the first year of life. We categorized the number of RI events into 0–2 vs. 3–7 RIs, based on previous studies of the respiratory microbiome within this same cohort17,18. In these studies, RIs were initially categorized into three groups based on the normal distribution of this variable. The 0–2 RIs group was found to have the most stable development of the nasopharyngeal microbiota when compared with children suffering from >2 RIs in the first year of life, and was defined as the healthy reference group. While there were no correlations between microbiota composition and GI complaints, we observed an association between microbiota composition at 1 week of life and the categorized number of RI events (0–2 vs. 3–7 RI events: PERMANOVA test, R2 = 0.033, adjusted p-value = 0.028) as well as number of antibiotic courses prescribed over the first year (R2 = 0.024, adjusted p-value = 0.055). These two outcomes were related, as the antibiotics prescribed were mostly indicated for RIs. We next aimed to identify the taxa explaining this association between microbiota composition at 1 week and fewer RI events later in life by cross-sectional differential abundance analysis, while adjusting for mode of delivery. We observed, among others, Bifidobacterium to be associated with fewer RI events (0–2 vs. 3–7 RI events, zero-inflated Gaussian mixture model, log2 fold change (log2FC) 2.118, adjusted p-value 0.049, Fig. 5), whereas Klebsiella and Enterococcus were negatively associated with fewer RI events (log2FC −3.242, adjusted p-value = 0.007 and log2FC −2.838, adjusted p-value = 0.009, respectively). Other taxa found to be negatively associated with fewer RI events encompassed genera such as Veillonella and Staphylococcus. Random forest analysis was used to verify these results and identified once again Enterococcus, Bifidobacterium, and Klebsiella as the most important taxa driving the prediction of the categorized RI events in the first year of life (Supplementary Table 4). Furthermore, a stratified analysis for the VD and CS groups separately, showed similar associations between gut microbiota composition at 1 week of life and number of RI events for the two groups (PERMANOVA test, R2 0.003 and 0.005, p-value 0.040 and 0.068, respectively). Taxa associated with the number of RI events were comparable between the overall and stratified analyses, although for VD children we now only found significance for Enterococcus (zero-inflated Gaussian mixture model, log2FC −2.525, adjusted p-value 0.074), whereas for the CS children, Bifidobacterium (log2FC 2.805), Klebsiella (log2FC −6.991) and Enterococcus (log2FC; −4.283) were all significantly associated with number of RIs (adjusted p-values 0.055, 0.010, and 0.055, respectively).

Fig. 5 Differentially abundant taxa between 0-2 vs. 3-7 RI events in first year of life. To identify taxa that were differentially abundant between children experiencing more vs. limited respiratory infection (RI) events over the first year of life, fitZig analysis was performed on the 119 samples obtained at week 1 with the rare-features-filtered OTU-table containing 97 taxa and contrasts set to 0–2 vs. 3–7 RI events in the first year of life. The blue data points indicate taxa that were significantly more abundant in children having 0–2 RI events, while red points represent taxa that were significantly more abundant in children with 3–7 RI events in the first year of life. The results of two data points falling beyond the limits of the plot: *Eggerthella log2FC 3.512, adjusted p-value (log10) 7.357, **Erysipelotrichaceae log2FC 6.912, adjusted p-value (log10) 6.222, calculated using a zero-inflated Gaussian mixture model. Source data are provided as a Source Data file Full size image

Validation of the results with metagenomics and targeted qPCR

To validate our primary findings independently, we executed whole genome shotgun (WGS) sequencing on a subset of 20 randomly selected samples collected at 1 week of life from 10 VD and 10 CS born children. WGS sequencing yielded a total of 119 unique bacterial taxa. The relative abundances of the top 12 OTUs and species of both sequencing methods are represented in Supplementary Fig. 5, and show highly comparable profiles. The most abundant Bifidobacterium species in the WGS data set were B. longum, B. breve, and B. adolescentis. The combined relative abundances of these three species strongly correlated with the most abundant Bifidobacterium of the 16S rRNA data set (Pearson’s r 0.95, adjusted p-value < 0.001). In this way, we could also correlate the E. coli, Staphylococcus, Klebsiella, and E. faecium OTU abundance of the 16S rRNA data set with high certainty to the E. coli, S. epidermidis, K. oxytoca, and E. faecium species abundance in the WGS data set (Supplementary Table 5).

We also used the WGS sequencing data to validate the differences found by 16S rRNA sequencing in overall gut microbiota composition between VD and CS born children at 1 week of life. The results from the ordination using the WGS sequencing data are shown in Supplementary Fig. 6. Again, we found a significant effect of delivery mode on the overall gut microbiota composition at 1 week of life using WGS sequencing data (PERMANOVA test, R2 = 0.125, p = 0.01).

Next, we compared the microbiota profiles obtained by WGS sequencing between VD and CS children and observed that the combination of B. breve, B. longum, and B. adolescentis (median relative abundance 72.2% in the VD vs. 0.074% in the CS born children, Wilcoxon test, p = 0.002), K. oxytoca (<0.001 vs. 0.006%, p = 0.153) and E. faecium (0.014 vs. 0.035%, p = 0.023) were differentially abundant between groups. In this data set, E. coli and S. epidermidis did not differ significantly between the delivery mode groups.

To alternatively validate our results on the overall cohort and in a targeted manner, we performed qPCR analyses for E. coli, Klebsiella spp. and Enterococcus spp. on all week 1 samples (n = 119). The qPCR results confirmed that E. coli is more commonly present in VD children compared with CS children (chi-square test, p < 0.001), whereas CS children are more often colonized with Klebsiella spp. (p = 0.011) and Enterococcus spp. (p = 0.004) than VD children, corroborating the 16S rRNA and WGS sequencing results (Supplementary Table 6). Finally, qPCR also confirmed that colonization with Enterococcus spp. and Klebsiella spp. at 1 week of life was positively associated with more RI events in the first year of life, though this difference was only significant for Enterococcus spp. (p = 0.015, Supplementary Table 7).