This study was approved by the Biomedical Ethics Committee of the Institute Pasteur of Tunis, under reference number 2018/10/I/LR161IPT09. Written informed consents were obtained from both parents.

All participants originated from Zaghouan, a semirural region in Northwest Tunisia. The mothers had a non-declared pathology and did not receive antibiotics, anti-inflammatory drugs, or immunosuppressive therapy during the last trimester of pregnancy. Neonates were healthy, full-term born (>37 weeks), and did not receive any antibiotic treatment. Only newborns delivered by elective and unlabored cesarean were included in this study. No complications during the delivery were reported. Anonymized mother's health records with their corresponding metadata were compiled. Clinical metadata concerning the newborns were also collected including birth weight, gestational age, gender, weight, and feeding mode (Table 1).

A total of 10 participants were retained, and newborns' stool samples were collected at three time points. The neonatal meconium from VD and ECS newborns was collected during the first 24 h (Day 0) at the regional hospital of Zaghouan and at the Aziza Othmana hospital, located in Zaghouan and Tunis, respectively. The collection was renewed for 2 weeks (Day 15) and 1 month (Day 30) after birth. The samples were transported in ice bags to the laboratory and then stored at −80°C until DNA extraction.

Subsamples from each newborn stool (~200 mg) were used for DNA extraction. The solubilization and liquefaction of the samples were carried out using 1 ml of Tween 20 (10%). The solution was then vortexed horizontally for 20 min and centrifuged at 14,000 rpm for 10 min at 4°C. The obtained pellets were resuspended in a lysis buffer (100 mM NaCl, 25 mM EDTA, 0.5% SDS, 10 mM HEPES, and 1 mg/ml lysozyme) for chemical and enzymatic lysis. Then, bead-beating steps were proceeded using BeadBug and PowerBead Tubes, Garnet 0.70 mm (QIAGEN), with a setting of 3 × 20 s pulses, followed by incubation for 7 min at 95°C. After a centrifugation step, the supernatant was transferred into 2 ml tubes containing 30 μl proteinase K and 400 μl AL buffer (QIAamp Fast DNA Stool Mini Kit) and incubated at 70°C for 10 min. DNA was recovered using phenol-chloroform extraction and purified using the QIAamp DNA mini kit (QIAGEN). In brief, protein precipitation was performed using the phenol/chloroform method, followed by DNA precipitation using absolute ethanol. The supernatant was transferred to the QIAamp spin column, and DNA purification was performed according to the QIAamp Fast DNA Stool Mini Kit manufacturer's instructions. DNA samples were quantified using the Multiskan™ GO Microplate Spectrophotometer (Thermo Fisher™), and samples with >100 ng/μl DNA material proceeded to the paired-end (2 × 150 bp) metagenomic sequencing. Overall, we sequenced 33 samples: 9 participants at 3 time points, 1 participant at 2 time points (Day 0 and Day 30), and 4 technical replicates (2 at Day 0, 1 at Day 15, and 1 at Day 30).

Quality control and preprocessing were undertaken on the raw data as follows: Trimmomatic was run to trim low-quality bases and short reads (SLIDINGWINDOW: 4:20, min length: 70) (Bolger et al., 2014). Following this step, we obtained on average 11 and 14 M reads for Day 15 and Day 30 samples. As expected, Day 0 samples harbored fewer sequences, 28.6 K reads on average (Supplementary Table 1). Then, human DNA and vector contaminants were screened and removed using FastQ Screen software (Wingett and Andrews, 2018). The taxonomic classification of cleaned metagenomic reads was performed using Kraken version 2.0.8 against the Minikraken Database (Wood et al., 2019). Then, Bracken version 2.5 was run on the Kraken taxonomic assignment outputs to increase the accuracy of the abundance estimates at species and genus levels (Lu et al., 2017).

The classification output consisted of a feature table containing 3,805 taxa composed of bacteria 3,489 operational taxonomic units (OTUs) as well as viruses (212 OTUs) and Archaea (103 OTUs). The statistical analysis started with filtering ambiguous and rare bacterial taxa (Supplementary Table 1). First, low-prevalent (<15% samples) and low abundant (<100 reads) phyla were removed. OTUs present in <15% of samples were also excluded. None of the filtered OTUs were specific to a mode of delivery or unique to collection day, thus not impacting the overall rates of abundance (data not shown). Ultimately, 2,609 OTUs belonging to bacteria remained in the taxonomy table.

Alpha diversity metrics (Shannon and Simpson indices and OTU richness) were analyzed using the microbiome R package (McMurdie and Holmes 2013).

A phylogenetic tree was generated and used to estimate Faith phylogenetic diversity (standardized effect sizes = TRUE) using the function “phyloseq_phylo_div” from the R package PhyloMeasures (Tsirogiannis and Sandel, 2016) and metagMisc R package (available at: https://github.com/vmikk/metagMisc). The Friedman test was used to compare differences in alpha diversity between collection days (paired samples), followed by the Nemenyi post-hoc test using the PMCMRplus R package (available at https://cran.r-project.org/web/packages/PMCMRplus/index.html). When comparing alpha diversity between two unpaired groups, the non-parametric Wilcoxon was employed. Kruskal–Wallis test was used to assess statistical differences between more than two groups, followed by post-hoc Dunn's test (with Bonferroni correction). The clustering and heatmap were performed using the microbiomeutilities R package (Shetty and Lahti, 2020). The clustering based on taxa and samples was performed using Bray–Curtis dissimilarity distances.

Beta diversity was evaluated with the Bray–Curtis distance after data normalization using the Cumulative Sum Scaling (CSS) (Paulson et al., 2013). Visualization was performed using the Non-metric Multidimensional Scaling (NMDS) ordination method. Differences between the groups were tested using the permutational multivariate analysis of variance (PERMANOVA) test (number of the permutation set to 999). The linear discriminant analysis effect size (LEfSe) R package was used to identify the differentially abundant taxa between study groups on CSS normalized data (Segata et al., 2011) (parameters: kw_cutoff = 0.01, wilcoxon_cutoff = 0.01, and lda_cutoff = 3). PERMANOVA test was also run to determine covariate impact on newborns' gut microbial community based on Bray–Curtis distance (permutation set to 999).

Functional analysis was processed with SqueezeMeta Pipeline using the merged mode (Tamames and Puente-Sánchez, 2019). First, Megahit was used for the assembly (—min-count 1—k-list 27, 37, 47, 57, 67, 77, 87) (Li et al., 2015). We obtained 139,759 contigs ranging from 200 to 761,139 bp (Supplementary Table 1). Then, open reading frame (ORF)'s prediction and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (KEGG database, release 95, July 2020) prediction were performed by Prodigal and Diamond, respectively (Hyatt et al., 2010; Buchfink et al., 2015). Coverage and abundance estimation for genes and contigs were performed using Bedtools and normalized using the reads per kilobase per million mapped reads (RPKM) method (Quinlan and Hall, 2010). The statistical downstream analyses were performed using the SQMtools and phyloseq R packages (McMurdie and Holmes, 2013; Puente-Sánchez et al., 2020).

Full-term newborns (5 from vaginal delivery and 5 from ECS delivery) were delivered at a mean gestational age of 38.5 weeks (SD ± 1.22). Samples were collected on the birthday (Day 0), Day 15, and Day 30 after birth. The characteristics of the mothers and newborns are listed in Table 1. Overall, 40% of newborns were exclusively breastfed, 10% exclusively received formula milk, and 30% received mixed feeding. Notably, 80% of neonates were breastfed until the second time point (Day 15). By Day 30, this percentage decreased to 40%. There were no significant differences in maternal age [31 years (SD ± 4.70)], mother's weight [70.76 kg (SD ± 10.44)], and birth weight [3.43 kg (SD ± 0.34)] between ECS and VD (Table 1).

We first aimed to determine the significant factors associated with microbiota structure. Therefore, we performed the PERMANOVA test individually for each covariate. The results showed that collection days explained most of the variance in population structure (Supplementary Table 2; p = 0.001, R² = 0.210). Although not impacting the diversity, feeding mode was the second most significant factor explaining variation in neonatal gut microbiota (p = 0.033, R² = 0.115; Supplementary Figure 1, Supplementary Table 2). The other clinical covariates (newborns' weight and gestational age) did not influence the community structure of the early microbiome (Supplementary Table 2).

We then scrutinized the early establishment of Tunisian newborns' gut microbiota by comparing the three sampling time points. This allowed us to identify the temporal patterns in terms of taxonomic diversity. A significant decrease in terms of Shannon evenness and Simpson diversity indices was observed when comparing Day 0 with Day 15 and Day 30 (Figure 1A, top left and top right panel). In addition, the number of observed OTUs and phylogenetic diversity significantly increases over time (Figure 1A, bottom left and bottom right panel). The beta-diversity analysis clearly and significantly distinguished the Day 0 samples from the other time points (Day 15 and Day 30) (p = 0.001, R² = 0.210; Figure 1B). The functional annotation analysis corroborated this ascertainment (Figure 1C; p = 0.001, R² = 0.688). This suggests that the functions potentially brought by bacteria at Day 0 are distinct from those ensured at Day 15 and Day 30. At the phylum level, a pattern of increasing abundance composed of Firmicutes, Proteobacteria, and Actinobacteria is observed (Figure 1D). Clustering at the genus level revealed that the abundance of Escherichia, Klebsiella from the Proteobacteria, Bifidobacterium from the Actinobacteria, and Streptococcus from the Firmicutes increases from Day 15 post-delivery (Figure 1E). The LEfSe analysis highlighted a significant increase for the above-mentioned genera and phyla (Supplementary Table 3).

To address the impact of DM on newborns' gut microbiota, we carried out various taxonomic and functional diversity analyzes. Alpha and beta-diversity showed no significant difference neither in terms of diversity nor in bacterial community structure between ECS and VD groups (p = 0.331, R² = 0.034) (Figures 2A,C). Similarly, we did not observe differences in community function between ECS and VD (Figure 2B; p-value = 0.75, R² = 0.011).

Although the ECS and VD microbiota shares a similar community structure, marginal taxonomic differences do exist. To further characterize these differences, LEfSe analysis was performed to detect specific differentially abundant features at various taxonomic ranks in ECS and VD (Supplementary Figure 2). In the ECS group, species belonging to Actinomyces, Pseudopropionibacterium, and Citrobacter genera were overrepresented compared with the VD group. Interestingly, the microbiota of ECS newborns was enriched in the genus Clostridium. Among these genera, the Clostridium perfringens, an opportunistic pathogenic species, was detected as overrepresented. In addition, the ECS showed an overrepresentation of Enterobacter species, a highly virulent and antibiotic-resistant member of the ESKAPE bacterial group.

The Bacteroidetes phylum, including species belonging to the Bacteroides genus (Bacteroides helcogenes; Bacteroides heparinolyticus), and Parabacteroides were identified as overrepresented in the VD group (Supplementary Figure 2). The ECS group also showed an enrichment in the Bifidobacterium dentium, while the Bifidobacterium pseudocatenulatum was overrepresented in the VD group.

Taken together, these results indicate that although the mode of delivery does not apparently impact the overall diversity, there is nevertheless a significant difference at phylum and genus level between the two groups. These differences are mirrored by the enrichment of opportunistic pathogenic species in ECS and Bacteroidetes in VD.

We then sought to understand how the composition and developmental trajectory of the early gut microbiota occurs in ECS and VD. A transition that occurs within the first 2 weeks of life was observed in both ECS (p = 0.002, R² = 0.246) and VD (p = 0.002, R² = 0.269; Figure 3).

While the gut microbiota of ECS neonates was colonized by Actinobacteria, Proteobacteria, and Firmicutes, the microbiota of the VD group also included the Bacteroidetes phylum. The Bacteroidetes were represented exclusively by the genus Bacteroides, while the genus Raoultella was strictly observed in ECS neonates (Figures 3A,B, Supplementary Table 4). We also observed a gradual increase of Actinobacteria at the expense of Proteobacteria (Figure 3A). This trend appears to be a consequence of the progressive growth of the genus Bifidobacterium, at a lower rate in ECS, and the decline of the genera Delftia, Enterobacter, Sphingomonas, and Stenotrophomonas (Figure 3B, Supplementary Table 3). With a lower relative abundance, the Firmicutes remain constant in the ECS while increasing in the VD, possibly due to expansion of the Veillonella genus (Figures 3A,B, Supplementary Table 4). At the genus level, the most striking difference between the two groups is the presence of Lactobacillus and Collinsella in VD.

Some genera present at birth do not persist in the remaining time points. This was observed for Cutibacterium, Acinetobacter, Delftia, Pseudomonas, Sphingomonas, Stenotrophomonas, and Methylobacterium, regardless of DM. While focusing on ECS newborns' meconium, we noticed that genera Citrobacter, Cupriavidus, Halomonas, and Paracoccus are undetected from Day 15. Similarly, Gardnerella and Methylorubrum are solely present in VD newborns' meconium during the first 2 weeks.

Again, LEfSe analysis was used to detect overrepresented bacterial taxa between ECS and VD at each sampling time point. We found that few distinctive taxa signatures were detected in the neonatal microbiota at Day 0 (Figure 4A). While the Clostridiale order and the Citrobacter genus were overrepresented in ECS, the Methylobacter genus was detected in VD.

On Day 15, Proteobacteria, especially Escherichia and Shigella, were overrepresented in VD newborns compared with the ECS group (Figure 4B). The ECS group had a higher relative abundance of Clostridium and Clostridiaceae. In addition, the ECS group was enriched in species belonging to the Bifidobacterium genus, namely B. dentium and Bifidobacterium actinocoloniform (Figure 4B).

On Day 30, LEfSe analysis revealed an overrepresentation of species belonging to Bacteroides, Parabacteroides, and Gardnerella genera, respectively, represented by Bacteroides ovatus, Parabacteroides distasonis, and Gardnerella vaginalis, respectively, for the VD group (Figure 4C). The ECS samples showed an enrichment of Klebsiella oxytoca and various Enterobacter species (Figure 4C). Given the prevalence of these opportunistic pathogens in the gut of ECS neonates, we extended this analysis and compared the relative abundance of the ESKAPE group. We found that the ECS group had a significantly higher relative abundance at Day 15 and Day 30 (Supplementary Figure 4).

The above-mentioned results indicate that a shift occurs during the same timeframe in ECS and VD. However, ECS showed a disturbed establishment of the Bacteroides genus and an increased relative abundance of ESKAPE bacteria members from the second week.