Abstract

Recently, the intrauterine sterile environment theory has been questioned. Growing evidence shows that microbial in utero pioneer gut colonization could occur prebirth, and this initial colonization may play an important role in the development of the neonate immune system and setting up a niche for the adult-like microbiota. In this study, we compared the microbiota networks from public available meconium datasets from different countries. The findings showed differences at the genera level and were country-dependent. We generated and analyzed bacterial networks, at the genera level of meconium samples from c-section and vaginally delivery modes. Interestingly, bacterial networks from the c-section-delivered meconium samples tended to have a bigger diameter but fewer correlations, whereas the vaginally delivered meconium networks were smaller and with a higher number of correlations. Even more, the networks were similar in the delivery mode, even between countries, at the genera level. The c-section networks suggest incomplete colonization or important lack of bacteria, promoting the susceptibility of the network to receive new members, beneficial or pathogens. These results suggest that the network analysis contributes to the knowledge of microbiota composition, identifying microbial associations, despite the differences between the environment and country habits, and obtaining a better understanding of microbial gut colonization.

1. Introduction

According to the Developmental Origins of Health and Disease (DOHaD) theory, the origin of adult health or disease could be established during fetal development [1]. During the intrauterine developmental phase, the environmental stimuli associated with maternal health like diet and nutritional state, weight, stress, physical activity, pollution-exposed habits like smoking, and alcohol consumption, among others, determine the conditions for developing the fetus. Otherwise, all that stimulus could lead to epigenetic markers that are part of fetal programming for health or disease [2].

Since the report of Aagaard and collegues, in 2014, about the placental microbiota composition [3], the belief of a sterile intrauterine environment was debated. The findings showed a particular and unique characteristic microbial niche with differences to others, such as the vagina, mouth, gut, and skin [3]. These results have been discussed due to possible contamination during the taking and processing of the sample and laboratory contamination [4].

In 2019, Li et al. published a study where using negative controls demonstrated that the first bacterial colonization happens prebirth and before the newborn can even have contact with any surface; therefore, the findings are not a contamination issue. Moreover, they proposed that the way of birth, vaginal or cesarean delivery, could not influence this type of colonization [5].

Kimura and collaborators reported in a mouse model that the maternal gut microbiota confers resistance to obesity in offspring via the SCFA-GPR41 (Short Chain Fatty Acid-G protein-coupled receptor 41) and SCFA-GPR43 (Short Chain Fatty Acid-G protein-coupled receptor 43) axes. SCFAs (Short Chain Fatty Acids) produced by maternal microbiota activate embryonic GPR41 and GPR43 receptors. Activation of these receptors promotes sympathetic, neuronal, enteroendocrine, and pancreatic B cell differentiation. The conclusions were that maternal microbiota influences the offspring’s metabolic phenotype [6].

Meconium is the initial substance present in the intestines of the developing fetus and constitutes the first bowel movement of the newborn. 24 to 48 h following birth, term healthy neonates pass the meconium and represent a noninvasive method for sampling and investigating the gut colonization [7].

Tapiainen and collaborators described interindividual variability among first-pass meconium samples. They found that the most abundant phyla were Firmicutes, with a relative abundance of 44%; in the second place, Proteobacteria 28%; and the third one, Bacteroidetes 15%. They also reported that the microbial characteristics of meconium were not affected by immediate perinatal factors, like the delivery mode and the use of antibiotics, whereas the greater biodiversity of the maternal living environment (such as pet cohabitation) during all three trimesters of pregnancy increases the diversity of the microbiota, which they described as implying an in utero transfer of microbes [8].

The most abundant genera in the gut of a healthy newborn are Bifidobacterium, Veillonella, Streptococcus, Citrobacter, Escherichia, Bacteroides, and Clostridium, according to Milani and colleague's review from 2017 [9], and at the phylum level are mostly Proteobacteria, Firmicutes, Actinobacteria, Verrucomicrobia, and Bacteroides as reported by Del Chierico and colleagues [10].

The neonatal gut microbiota is considered unstable because it undergoes several changes during the first months of life. Facultative anaerobes are the first colonizers to prepare the niches for the anaerobe microbes [11]. Infant gut microbiota shifts due to several factors, such as delivery mode, maternal and infant feeding, the introduction of solid food, the use of antibiotics, and the presence of pets in the household [12, 13].

In 2019, Reyman and colleagues found that Bacteroides phyla were not found in newborns who were delivered by c-section, in contrast with those vaginally delivered (Reyman et al.). Recently, Yassour and colleagues analyzed 73 samples divided into 3 groups: vaginally delivered, postlabor CS (c-section), and prelabor CS. They described the presence of Bacteroides in all 3 groups, but in those born via c-section (pre- or postlabor), the colonization could not be maintained; then, the Bacteroides disappeared at week 2. These results suggest a new hypothesis that both birth delivery modes, vaginal and c-section, promote Bacteroides colonization of the newborn gut; however, the necessary conditions that allow the maintenance of an appropriate niche for Bacteroides are not present in the gut of newborns that were born by CS [14].

The implications of these findings in the changes in the microbiome composition due to delivered mode are discussed and still controversial; however, much evidence supports the main role of the composition, abundance, and function in the microbiome in human health Kulagina et al. [15]. Highlighting the differences of the microbiome over the early steps of life that could determine the immune system development, maturation, and function and how is associated with the susceptibility to diseases during life, even in adults [14, 16]. The gut is associated with neurodegenerative, cardiovascular, and different metabolic diseases. Mainly, the conclusions show a strong association between differences in microbiome between healthy individuals and patients ([17, 18], Méndez-Salazar et al. 2018)

Currently, the controversy around in utero colonization continues mainly because of the lack of reliable methods that allow researchers to confirm whether human bacterial colonization begins during intrauterine development but also due to ethical restrictions and technical issues, like better sample selection and methods to obtain it and contamination issues, among others.

Microbiota composition studies (16S rRNA sequencing) mostly focus on measuring the abundance (raw and relative counts) and the presence or absence of microorganisms in a specific environment. Moreover, there are different platforms to perform 16S rRNA sequencing, but the raw sequencing data should be comparable with the independence of the platform in which it was conducted Allali et al. [19]. Unfortunately, this does not happen, and most of the time, some bioinformatic preprocessing must be done to get comparable data.

Bioinformatic tools bring the possibility of performing various analyses using existing datasets or databases. One of these analyses is the association of microorganisms under different conditions, measured as correlations and generating microbial association networks. Network analysis of microbial association patterns allows getting insights into the characteristics of the bacterial communities and possible interactions, among others [20]. Beyond the abundance of different bacteria, the understanding of the functional role and associations between microorganisms may help us get a deeper understanding of the microbiota composition.

Since the first gut colonization is still a controversy, whether it occurs during intrauterine development or after birth, many articles have reported the sequencing results of the meconium microbiota. In this study, we selected three different publicly available datasets from meconium samples from three different countries, to compare the characteristics between pioneer gut microbiota and between meconium microbiome association networks to better understand the associations between microbial communities. We found that at the phylum level, the microbiota composition is similar between delivering modes, vaginal or c-section. However, at the genera level, all three studies present differences between them. We generated association networks for the meconium samples at the genera level for each study and for each delivery mode. We found common features depending upon the delivery mode, regardless of the study. C-section networks are less compact and do not have a dominant community, compared to vaginal networks. Vaginal networks present a dominant community structure and are highly connected. Our results highlight the utility of network analysis in microbiota studies for assessing a better understanding of the first gut colonization.

2. Materials and Methods

2.1. BioProject Datasets

BioProject datasets PRJNA311499, PRJNA530829, and PRJNA559967 were downloaded from the Sequence Read Archive. We have used the fastq-dump function from the SRAtoolkit (GitHub—ncbi/sra-tools: SRA Tools) with default parameters and with the --split-3 parameter. We performed a quality control analysis using the FastQC v0.11.9 tool (Babraham Bioinformatics—FastQC, a quality control tool for high-throughput sequence data). For the PRJ311499 (Ion Torrent PGM platform) project, the fragments have a mean of Phred score of ~29. The PRJ530829 (PacBio platform) project fragments have a mean Phred score of ~92. And the PRJ559967 (MiSeq Illumina platform) project fragments have a mean Phred score of ~37.

2.2. 16S rRNA High-Throughput Sequencing Data Processing

We used the R package DADA2 v1.16.0 [21] to process the 16S rRNA sequencing data. For the PRJNA311499 (Ion Torrent) and PRJNA530829 (PacBio) datasets, we have trimmed 10 nt from the left and right sides of the fragments. For PRJNA559967 (MiSeq Illumina), we have trimmed 25 nt on the right side. For the PRJNA530829 project, we have removed the forward primer AGRGTTYGATYMTGGCTCAG and the reverse primer RGYTACCTTGTTACGACTT using the removePrimers() function, with default parameters. Error rates were calculated using the learnErrors() function with default parameters for all BioProject datasets; with HOMOPOLYMER_GAP_PENALTY = -1 and BAND_SIZE = 32 for PRJNA311499 and errorEstimationFunction = PacBioErrfun and BAND_SIZE = 32 for PRJNA530829. We have used the dada() function to infer the sample composition using default parameters and the calculated error rates. Then, we removed the chimeras using the removeBimeraDenovo () function with default parameters. The taxonomy for each sequence was assigned using the assignTaxonomy() function and the silva_nr99_v138 database. For all BioProject datasets, we have obtained an amplicon sequence variant (ASV) table.

2.3. Taxon Abundance Analysis

Taxon abundances were processed using custom R v4.0.2 (R Core Team [22]—European Environment Agency) scripts. For each dataset, we have removed samples with no counts in at least a fifth of the samples for PRJNA311499, one sample for PRJNA530829, and at least a fifth for PRJNA559967. We have removed taxa with less than one count in at least a twentieth of the taxa. The normalization for the sample in each dataset (project) was calculated as the logarithm of the counts per million for each sample, .

2.4. sPLS-DA Analysis

Using the normalized abundances, we have performed an sPLS-DA analysis to identify which taxa are more abundant in different sample types. We used the splsda() function from the MixOmics R package [23] with default parameters, and .

2.5. Permutation Test for sPLS-DA Loading

To obtain statistical significance to the loadings from the sPLS-DA, we have performed a permutation test with an in-house R script and by randomizing the variable ids, genera, or pathways according to the analysis. The number of permutations was . For each variable, we calculate the observed score; then, we randomize the variables ids times, and for each randomization, we obtain a randomized score. Finally, we obtain a value by comparing the observed score to the randomized scores.

2.6. Network Inference

We inferred the relationship between taxa for each sample type in each BioProject dataset as the correlation between taxon-normalized abundance. Each taxon is defined as a vertex (vertex/node: a vertex is defined as the fundamental unit of a graph). The correlations are defined as the edges of the network (edge: an edge is the link between any pair of vertices). We used the cor() function from the stats R package with the “Spearman” method. For each correlation , we tested it with the cor.test() function and obtained its corresponding value. We have corrected these values with the p.adjust() function using the method False Discovery Rate (FDR). We kept correlations with in order to reduce the false discoveries. Using these correlations between taxa, we used the igraph R package (Package “igraph” Title Network Analysis and Visualization [24]) to generate an igraph network (network: a multigraph in which more than one edge is allowed between a pair of vertices and edges may be connect a vertex to itself) object with . Network properties (diameter: the maximum distance over all vertices in a network, degree: the number of edges corresponding to each vertex, and hub score: score associated to a node that with a number of links greater than the average) were calculated using the diameter(), degree(), and hub_score() functions from the igraph package. Plots of the networks were generated with the plot() S3 method for igraph objects (Package “igraph” Title Network Analysis and Visualization [24]).

2.7. Pathway’s Inference

We have used Tax4Fun2 [25] to perform the pathway prediction for the 16S sequencing data. First, we used the runRefBlast() and then the makeFunctionalPrediction() functions, with default parameters. The input to these functions was the ASV quantification profiles. We used the reference RF99NR used in both steps. The resulting pathway prediction table was used to perform a sparse Partial Least Squares Discriminant Analysis (sPLS-DA). To perform this, we used the function spls-da() from the R package MixOmics with default parameters and [23]. We have used in-house R scripts to generate the plots to represent the results of the sPLS-DA (R: The R Project for Statistical Computing [26]).

3. Results

A total of three datasets with 16S rRNA-sequenced meconium samples was analyzed using bioinformatic methods (see Materials and Methods). All three datasets included data from meconium samples that were divided into two groups: data of newborn samples delivered by c-section (CS) or vaginally (V). These studies were selected based on the availability of the data, and the type of sample was meconium from newborns from different studies and countries.

3.1. The Phylum-Level Composition Is Consistent across Studies

Relative abundances at the phylum level were calculated to describe the composition by delivery mode in each study. In Figure 1, we observed the previously reported phylum-level composition in all meconium samples with some differences in the order, with Proteobacteria and Firmicutes being the most abundant phyla [9, 10, 13].

In the Oulu, Finland, study (PRJNA311499) [8], at the phylum level, we found 41% Firmicutes, 31% Proteobacteria, 25% Bacteroides, and 1.3% Actinobacteria in the meconium samples from vaginally (V) delivered newborns. However, for c-section delivery (CS), the samples contain 34% Firmicutes, 65% Proteobacteria, 0.015% Bacteroides, and 0.66% Actinobacteria composition. Interestingly, Proteobacteria in CS increases ~2.1-folds, and Bacteroidota reduces ~0.0006-folds.

The Perth, Australia, study (PRJNA530829) [27] uses only c-section- (CS-) delivered newborn samples. The composition was 3.9% Firmicutes, 92% Proteobacteria, 3.3% Actinobacteria, and 0.31% Deinococcota. As in the Finland study, CS microbiome composition was dominated by Proteobacteria and with no detectable fraction of Bacteroidota.

In the Yunnan, China, study (PRJNA559967) [5], the vaginally (V) delivered neonate meconium samples had 11% Firmicutes, 59% Proteobacteria, 11% Actinobacteria, and 17% Deinococcota. The c-section (CS) samples had 8% Firmicutes, 70% Proteobacteria, 4.3% Actinobacteria, and 12% Deinococcota. This Yunnan study has undetectable Bacteroidota evidence for both types of samples, V and CS. Moreover, compared to the other two studies, Firmicutes is not the most abundant phylum. Proteobacteria dominate in both types of samples and are more abundant in CS samples.

Expected phylum distribution was founded in all three studies, with the predominance of Firmicutes and Proteobacteria. These phyla have the characteristics needed to exist in the neonatal gut environment, although they are composed by Gram positive and negative bacteria, and the presence of both aerobic and anaerobic species.

3.2. Differences by Country and Birth Mode

According to the multidimensional scaling (MDS) analysis (Figure 2), there is no tendency to cluster together with the samples according to the birth mode at the genus level. We observed that samples from PRJNA559967 (China) were clustered together (both c-section and vaginal-delivered samples). Samples from PRJNA311499 (Finland) were not clustered with any other study samples, but V and CS samples are differentiated between them. In the PRJNA530829 (Australia) study with only CS samples, these samples are closer to the PRJNA311499 (Finland) V samples. From these results, a clear clustering is not observed given the birth delivery mode. In the PRJNA559967 (China) study, differentiation between birth modes was not observed. Thus, no similarities in the relative abundance composition were observed in the meconium samples analyzed.

3.3. Diverse Genera Are Present by Study and Delivery Mode

Sparse Partial Least Squares Discriminant Analysis (sPLS-DA) and permutation test were used to obtain the significant genera that characterize each delivery mode in the three datasets (Figure 3). In the Oulu, Finland, study (PRJNA311499), the genera from the CS-delivered meconium samples explained ~56% of the total variance. In this study, there was a clear differentiation between V and CS meconium samples (Figure 3(a)). The specific and significant ( value < 0.1) genera for CS (Figure 3(b)) were Bacteroides, Agathobacter, UCG-002, Faecalibacterium, Dialister, Alistipes, and Anaerostipes. The specific and significant genera for V were Bradyrhizobium and Staphylococcus. The genera associated with the differentiation of the CS samples are more mostly Gram positive. Instead, the V sample genera are mostly found in plant, root, and skin microbiota, the dominant genus being Gram negative.

The Yunnan, China, study (PRJNA559967) is shown in Figure 3. In this study, there is no clear differentiation between V and CS samples as the Oulu, Finland, studies (PRJNA311499) (Figure 3(c)), according to the sPLS-DA analysis. Here, the specific and significant ( value < 0.1) genera generating the mild differentiation are Streptococcus, Thermus, Pyrinomonas, Schlegelella, Brevundimonas, and Amaricoccus. Instead, the specific and significant genera for V were only Alcaligenes. However, Deinococcus and Megasphaera were also detectable but not significant (Figure 3(d)). The genera associated with the differentiation of the CS samples are mostly Gram negative, as well as the V samples.

The Australian study (PRJNA530829) only included data from the meconium of CS and amniotic fluid (AF) (Figure 3(e)). There was a mild difference between CS and AF samples. The Pelomonas and Novosphingobium genera were specific and significant findings ( value < 0.1) (Figure 3(f)). The genera associated with the differentiation of the CS samples are Gram negative and possibly related to the Australian water, as was before described [28].

We found common genera between the PRJNA311499 and PRJNA559967. These studies have V and CS samples. For the CS samples, the common genera were Ruminococcus and Lactococcus. For the V samples, the common genus was Bradyrhizobium. For both studies, there are detectable shared genera, but they are uncommon for the delivery mode, CS or V. These genera were Cutibacterium, Streptococcus, Stenotrophomonas, and Staphylococcus. Between the studies PRJNA311499 and PRJNA530829, there are common genera associated with V and AF, respectively. These genera were Ralstonia, Streptococcus, Staphylococcus, and Pelomonas (Figure 4).

Accordingly, these results showed that there is no shared genera for CS or V samples between studies. This suggests that each study has characteristic and specific genera due to the origin of the studies. The geographic, ethnic, and diet characteristics could explain the differences between the studies.

3.4. Network Characteristics

To better understand the relationship between taxa, network analysis was performed for each study and delivery mode. The characteristics of the networks are listed in Table 1.

In the PRJNA311499 study, the network of meconium CS samples has 32 nodes, 246 edges, and a diameter of 0.92. Alternatively, the network of meconium V samples has 393 edges, 32 nodes, and a diameter of 0.62. Besides, the V network has a higher average degree, ~22 for V and ~17 for CS networks, and the average hub score, ~0.74 for V and ~0.60 compared with the CS network. V network has 5 communities, the main community contains 14 genera, and the second one contains 9 genera. In contrast, the CS network has only 2 communities containing 14 and 15 genera. Besides, the V network has a greater number of positive correlations between genera, 88% vs. 84% (Table 1). The results suggest that the V network is more compact, dense, and dominated by only one community, compared to the CS network.

In the PRJNA559967 study, the meconium samples of newborns from the CS delivery network have 28 nodes, 45 edges, and a diameter of 4.1. In contrast, the meconium samples of newborns from the V delivery network have 33 nodes, 126 edges, and a diameter of 2.6, having a higher average degree, ~5 for V and ~3 for CS networks, and an average hub score, ~0.25 for V and ~0.24 for CS networks, compared with the CS network, as well as the PRJNA311499 study. The V network has 5 communities, the main community contains 17 genera, and the second one contains 6 genera. The CS network has 9 communities containing 10 for the main community and 4 for the second to the fourth communities. In this study, the V network has a greater number of positive correlations between genera, 96% vs. 69%. For this study, the results suggest that the V network is also more compact, dense, and dominated by only one community in comparison to the CS network, as the PRJNA311499 study.

The network analysis results of the PRJNA530829 study showed 19 nodes, 166 edges, and a diameter of 1; this project only includes meconium CS samples. In this study, there is no V network to be compared with the CS network. However, the AF network has 3 communities and is dominated by 1 community containing 16 genera and 2 communities with only 1 genus. The CS network has 2 communities with 6 and 13 genera. Here, the AF network is dominated by one community, as for the V networks in other studies.

Networks from the meconium samples from the vaginally delivered mode tended to have a greater number of edges in contrast to those coming from the c-section delivery; however, V networks had a smaller diameter. Interestingly, the findings on meconium samples from newborns vaginally delivered showed the dominance of only one community, given the number of genera that the communities contain, compared to the c-section networks. In contrast, CS networks do not show clear dominance of only one community.

The networks coming from the meconium data of newborns who had a vaginal delivery look centered and tighter, with stronger connections, with a greater number of correlations between genera present in the niche (), and most of the correlations were positive (Table 1). The networks resulting from the c-section-delivered samples seemed more scattered, with greater distance between the nodes, more negative correlations (see Table 1), and more scattered communities, compared to the vaginal networks. The characteristics observed suggest that the CS networks are more dynamic than V networks and prone to integrate new genera, beneficial or prejudicial, into the network.

3.5. Interaction Changes in the Newborn Core Gut Microbiota

The networks of PRJNA311499 are shown in Figure 5. Ruminococcus is one of the generas mentioned as the newborn Core Gut Microbiota and is present in both PRJNA311499 and PRJNA559967 networks (Figures 5 and 6). Ruminococcus represents a bigger node and has a higher number of correlations: 32 and 16, respectively; the majority were positive for vaginal delivery networks in both projects. In contrast, CS networks showed 15 and 3 positive correlations, respectively (Figure 6(a)). Ruminococcus is a genus of mutualistic bacteria and degrades complex polysaccharides, which could be useful in the feeding transition in the neonate.

Alternatively, the Bacteroides genus is only present in PRJNA311499 at V delivery, with 19 positive correlations and only 1 negative correlation with Phyllobacterium. Nonetheless, it is represented as a core node in this network (Figure 5). Moreover, mentioned among the newborn Core Gut Microbiota is the Lactococcus genus. It is present on both PRJNA311499 and PRJNA559967. In the PRJNA311499, it is in both CS and V delivery, while in the PRJNA559967 study, it is only present at the V delivery. In the PRJNA559967 analysis, Lactococcus is at the center as a relevant node in the network (Figure 7). Therefore, in this network, Lactococcus is a core node. Lactobacillus is present only in the PRJNA530829 CS network. All positive correlations were found in the network analysis for the PRJNA530829 study (Figure 6). Lactobacillus is also among the genera mentioned as newborn Core Gut Microbiota. The four genera described as part of the newborn Core Gut Microbiota were found in the networks. However, networks are constructed with the results of the sPLS-DA results; if some genera are not part of the networks, it does not mean that is not present.

3.6. Metabolic Pathway Prediction

Using the 16S rRNA sequencing data, the metabolic pathways related to the main genera for each study were predicted through the R package Tax4Fun2. The predicted metabolic pathways related to abundance at the genus level for each study are shown plotting sPLS-DA (Figure 8).

For the PRJNA311499 CS delivery mode, the significant pathways identified were lipoic acid metabolism, chlorocyclohexane and chlorobenzene degradation, and ether lipid metabolism. Moreover, detectable but not significant pathways were lipoarabinomannan LAM biosynthesis, alpha linolenic acid metabolism, styrene degradation, linoleic acid metabolism, and xylene degradation (Figure 8). For V delivery mode, the significant pathways identified were glycosphingolipid biosynthesis, lacro and neo lacto series, sphingolipid metabolism, glycan degradation, glycosaminoglycan degradation, various types of N glycan biosynthesis, glycosphingolipid biosynthesis globo and isoglobo series, and glycosphingolipid biosynthesis ganglio series (Figure 8). The predicted CS pathways are related to lipid metabolism, metabolism of cofactors and vitamins, and xenobiotics biodegradation mainly. The predicted V pathways are related to glycan biosynthesis and lipid metabolism.

In the PRJNA559967 study of CS delivery, the significant pathways were ascorbate and aldarate metabolism, arginine and proline metabolism, phosphonate and phosphinate metabolism, lipopolysaccharide biosynthesis, and nitrotoluene degradation. Other nonsignificant pathways for CS were taurine and hypotaurine metabolism, styrene degradation, nitrogen metabolism, steroid hormone biosynthesis (lipid metabolism), glutathione metabolism, and drug metabolism cytochrome p450 (Figure 9). For the V delivery mode, the significant pathways were cysteine and methionine metabolism, lysine biosynthesis, citrate cycle TCA cycle, D glutamine and D glutamate metabolism, central carbon metabolism in cancer, X2 oxocarboxylic acid metabolism, peptidoglycan biosynthesis, thiamine metabolism, phenylalanine tyrosine and tryptophan biosynthesis, lipoarabinomannan LAM biosynthesis, sphingolipid metabolism, C5 branched dibasic acid metabolism, carbon fixation pathways in prokaryotes, dioxin degradation, and oxidative phosphorylation. The CS pathways are related to carbohydrate metabolism, amino acid metabolism, glycan metabolism, and xenobiotic biodegradation. The V pathways are related to amino acid metabolism, carbohydrate metabolism, glycan biosynthesis, metabolism of cofactors and vitamins, lipid metabolism, and energy metabolism.

In the PRJNA530829 study of the CS delivery mode, the significant pathways identified were carbon fixation pathways in prokaryotes (energy metabolism), citrate cycle, TCA cycle, glycine, serine and threonine metabolism, thiamine metabolism, phenylalanine, tyrosine and tryptophan biosynthesis, methane metabolism, retinol metabolism, pyrimidine metabolism, glycerolipid metabolism, central carbon metabolism in cancer, pyruvate metabolism, glycolysis glyconeogenesis, peptidoglycan biosynthesis, oxidative phosphorylation from energy metabolism, and oxidative phosphorylation (Figure 10). These pathways are mainly related to energy metabolism, carbohydrate metabolism, amino acid metabolism, metabolism of cofactors and vitamins, nucleotide metabolism, and glycan metabolism.

The original database of the PRJNA530829 study excluded the vaginal delivery data. The significant AF pathways are related to the metabolism of other amino acids, lipid metabolism, metabolism of cofactors and vitamins, and energy metabolism (Figure 10).

There are no clear differences and common pathways associated with the delivery mode, vaginal, or c-section. However, we found that glycan biosynthesis (glycosphingolipids) and lipid metabolism (sphingolipids) were associated with V samples and xenobiotics biodegradation and carbohydrate metabolism in CS samples.

4. Discussion

Omics sciences such as metabolomics, metagenomics, transcriptomics, genomics, and many others generate a large amount of data Yu and Zeng [29]. Commonly, researchers cannot use all information; only one focus and one aim are usually reported. However, the data obtained could be used for many aims. In the microbiota analysis where 16S rRNA high-throughput sequencing is used, relative abundances are usually reported, and most of the time, another kind of analysis is not included. Even more, these studies are from different cities, regions, or countries [30].

Now, the growing progress in computational and bioinformatic tools allows us to dig deeper into the already existing data. The optimization of the information in public databases could be an opportunity for several deeper analyses that bring out new knowledge [31]. For the 16S rRNA sequencing data, in addition to relative abundance and ecological diversity, there are some other analyses to better understand the structure and dynamics of that ecosystem which could be performed. The correlation calculation, cooccurrence, or associations between taxa can be shown as networks, by using the normalized relative abundance between studies, for example. The characteristics such as network structure, network topology, and the potential ecological relationship between taxa could be obtained. Worthwhile information about the metabolic pathways, role, and activity of each genus in the niche could be elucidated with the network and pathway analyses [32].

The composition of the gut microbiota changes rapidly in the first years of development; the practice of breastfeeding leads the child to have an adult-like microbiota, sometimes regardless of the way of birth Kim et al. [33]. The relevance of the meconium microbiota study is the possibility of knowing about possible immunological conditions and general health information developed during the intrauterine stage [34]. The first colonizers made the initial immune imprinting [35]. Much evidence suggests that alteration of the first colonization is associated with the development of diverse diseases during extrauterine life, like autoimmune diseases, allergies, respiratory problems, and obesity, among others that have been previously related [34].

The delivery mode has been associated with specific patterns in the first gut microbiota colonization that may be associated with many diseases in the newborn and during extrauterine life [36, 37]. Previously, diminished bacterial diversity in newborns from CS delivery birth mode has been reported [38].

Here, we selected three studies with 16S rRNA throughput sequencing data that had publicly available databases. This is to compare the meconium microbiota according to birth delivery mode, vaginal (V) or c-section (CS). With the network analysis, several characteristics about the relationship between taxa were obtained and the possible contribution of each one over the niche.

First, we found that all samples had similar relative abundances at the phylum level despite country and delivery modes. All samples present the expected phylum distribution for newborn gut microbiota [10], with mild shifts like having either Firmicutes or Proteobacteria phyla with higher abundance. The fetal gut has specific characteristics like pH, oxygen, and nutrient availability [11]. This is the reason only some microorganisms can survive and remain. At this stage of life, the country and delivery mode of birth seem to not make much difference at the phylum level [5, 39].

At the genus level, groups from the studies were not clustered together (MDS analysis) by delivery mode or country. This could be an indication that, as previously mentioned, the delivery mode might not be the defining factor of the composition of the meconium microbiota, and at the genus level, there are not enough coincidences between delivery modes in the datasets to form clusters based on that characteristic.

From the previously reported in the newborn [40], we identified members of the microbiota from all 3 studies, using the sPLS-DA, that are characteristic genera dependent on delivery mode and country.

From the newborn Core Microbiota, Ruminococcus and Lactococcus are present in the CS of PRJNA311499 and PRJNA559967 studies. Ruminococcus is a strictly anaerobic genus, belonging to the Firmicutes phylum [41]. This bacterium degrades complex polysaccharides and transforms them into nutrients for the host [42]. As part of the newborn Core Microbiota [10], this genus is well adapted to the environment of the neonate’s gut. Lactococcus is a lactic acid genus that previously belonged to the Streptococcus genus until 1985, is Gram positive, and is facultative anaerobic gut bacteria, which is a part of the Firmicutes phylum as well [43].

Coincidental genera, such as Corynebacterium, Cutibacterium, Streptococcus, Staphylococcus, Stenotrophomonas, Lactococcus, and Ruminococcus, were identified in PRJNA311499 and PRJNA559967 (Figure 4). Corynebacterium and Cutibacterium are genera from the Actinobacteria phylum, commonly present in the skin microbiota and have been described as opportunistic pathogens [44]. This genus is characteristic of fermenting lactose into propionic acid in an anaerobic environment [45]. Cutibacterium has been identified as part of the pioneer gut colonizers [46, 47], both delivery modes c-section [48, 49] and vaginal delivery [50]. Streptococcus is a lactic acid producer that belongs to the Firmicutes phylum, like Lactobacillus and Lactococcus [51]. Stenotrophomonas is a Gram-negative bacterium and is an opportunistic human pathogen found mainly in soil and plants [52].

Bradyrhizobium was found in vaginal delivery meconium samples from PRJNA311499 and PRJNA559967. Bradyrhizobium is a Gram-negative soil bacterium with nitrogen fixation capabilities [53]. Stenotrophomonas and Bradyrhizobium are genera that belong to the phylum Proteobacteria, characteristic of the oral cavity [54] and found in soil and plant roots [55], respectively. The leguminous, possibly included in the maternal diet, might be responsible for the presence of the genus [56].

Thermus and Novosphingobium were identified in CS samples from PRJ530829 and PRJ559967 studies. Thermus is Gram negative and is implicated in meningitis, endocarditis, and septicemia. Thermus is considered a potential pathogen [57]. Novosphingobium is a Gram-negative bacterium that degrades aromatic compounds like phenol, aniline, nitrobenzene, and phenanthrene [58].

The differences at the genus level, depending on the delivery mode, could be explained by the differences in feeding habits of the mother [59], previous nutritional status [60], the use of antibiotics during pregnancy, the presence of pets in the household Kates et al. [61], and even the latitude, which is an environment factor [62]. All these factors are different for the country and for the culture. In this analysis, newborns from Finland, China, and Australia were included.

In Finland, potatoes, rye, dairy products, and sausages are commonly consumed Mikkilä et al. [63]. China is a large country with diverse eating habits; depending on the location, still, soy and soy products, rice, tea, seafood, noodles, different kinds of mushrooms, and various oils from different sources are part of the daily diet [57, 58, 64, 65]. Australians frequently consume various cereals, sausages, lamb, poultry, and seafood [66, 67]. The geographical characteristics of each country and the available food play an essential role in the establishment of gut microbiota [6870].

Although the delivery mode seems not to be definitive for the compositing of meconium microbiota expressed as relative abundance, when the correlation networks were analyzed, we found characteristics that differentiate both the CS and V delivery modes. Networks contain similar genera. However, the structure, conformation, edges, and nodes in each network have shown that the correlation between different genera are delivery mode-dependent. The network analysis information might contribute to the understanding of the way microorganisms relate to each other, their influence on niche development, the community organization, and the role of all microorganisms together, conforming the ecosystem for developing the immunological system and contributing to the susceptibility of several diseases.

Del Chierico et al. [10] reported a network analysis of meconium microbiota. Their findings showed fewer correlations in networks from vaginally delivered neonates than the c-section-delivered ones [10]. In contrast, in our analysis of the two datasets, the networks from vaginal delivery have more correlations than those from c-section. However, they obtained a lower proportion of negative correlations for the V network (~12%) compared to CS networks (43% for 1–3 days and ~35% for 7–30 days following), as in our results. Moreover, their V network looks more compact and potentially dominated by only one community.

The most interesting finding from our results is a different pattern of correlation in the networks between vaginal and c-section delivery modes. C-section bacterial networks have a lower number of total correlations; furthermore, a larger diameter (see Table 1) may be due to the incomplete colonization or important lack of bacteria, promoting that the network is “susceptible” to receive new members, beneficial or pathogens. Moreover, c-section networks have a higher percentage of negative correlations than the vaginal delivery mode analysis (see Table 1).

Accordingly, findings from the network analysis could contribute to information about the interaction between different genera and species and the kind of roles, functions, and activities of each one and all together in the niche. Changes in abundance, associations, and niches may be implicated in the susceptibility of individuals to the development of different diseases during adult life Libertucci and Young [71].

From the 16S data, we also predicted the metabolic pathways related to bacteria present in meconium samples. Using the sPLS-DA, the pathways that contribute to the difference between birth modes were identified. For the CS delivery, the xenobiotics biodegradation and metabolism pathway was found; during CS delivery, prescribed drugs are used; therefore, the xenobiotic pathways could be activated [72, 73]. Interestingly, differences at the genus level between delivery modes (CS and V) and in this study were found, but in the pathway analysis, we can observe some characteristic pathways, like xenobiotic metabolism in CS delivery and amino acid metabolism in the V delivery as reported here. The differences between the pathways reported could be due to the ability of different species from diverse genera to perform the same metabolic tasks. The characteristics of the bacteria present at the genus level might be determined by the maternal conditions and habits and the geographical location, but it is possible that in the end, they perform similar and needed tasks inside the neonatal gut.

The treatment in the laboratory of meconium samples is a challenging task. That could be the reason why few studies have reported using this sample. Issues with the DNA yield, extraction methods, and low biomass have been reported, as well as the need for a standardized method to work with these types of samples [7476]. However, meconium is a less invasive sample to analyze prebirth intestinal microbiota. Altogether, these results indicate that the study of meconium microbiota composition and the study of the first colonizers should use a holistic analysis. We must apply abundance compositions, diversity index, predictions of pathways and functions, and the organization of the microbiota via network analysis. All these analyses will help us to determine not only the specificity per study but also the similarities between studies coming from different populations and/or environments.

5. Conclusion

The network analysis of the microbiota obtained from 16S rRNA sequencing could contribute to important information about the niche, composition, and role of each phylum, genus, and species in the ecosystem, as well as the relationship between them. Using normalization steps in the bioinformatic methodology, the results will be more reproducible, and the comparisons are possible even with different sequencing platforms. Meconium samples could represent an accessible opportunity to study the in utero environment effect on fetal development and bacterial prebirth colonization mainly. The findings show that the networks showed different kinds of important conformation and are structure delivery mode dependent. However, more studies incorporating and comparing meconium and feces samples from the same individuals in longitudinal studies may give a better understanding of prebirth colonization, the dynamics of the associations between taxa, and its impact on adult health. Therefore, the network analysis brings evidence of the niche and relationships between the taxa present, increasing the knowledge about the function, activity, and influence of the ecosystem even with similar relative abundances according to the delivery mode.

Data Availability

The three BioProject datasets PRJNA311499, PRJNA559967, and PRJNA530829 are publicly available.

Conflicts of Interest

There are no existent conflicts of interest.

Acknowledgments

The authors acknowledge the Instituto Potosino de Investigación Científica y Tecnológica (IPICYT) and National Supercomputing Center (CNS) computing resources, provided to develop this project under the project grant TKII-R2018- COV1. This work was founded by the FORDECYTPRONACES Ciencia de Frontera 2019/101732 to COV. AKRV acknowledges Consejo Nacional de Ciencia y Tecnología (CONACYT) for the scholarship 2018-000068-02NACF-25236. COV acknowledges the Investigadores por México CONACYT program.