Characterizing the gut microbiome of early school-age Dutch children

To investigate the gut microbial characteristics of healthy early school-age Dutch children as well as the relationship between their gut microbiota and multiple phenotypic parameters, we collected stool samples from 281 children at 6–9 years of age (mean age 7.3 years) enrolled in the KOALA Birth Cohort Study [14], subjected the samples to shotgun metagenomic sequencing, and analyzed the results against additional measures of 45 phenotypic parameters classified into four categories: (1) early events; (2) pre-school lifestyle, including diet, both collected through repeated questionnaires up to 5 years of age; (3) blood parameters collected in parallel with fecal samples at 6–9 years of age; and (4) anthropometric measurements collected at both 4–5 years and 6–9 years of age (Fig. 1a, Additional file 1: Table S1). After performing shotgun metagenomic sequencing and quality control, we acquired a total of 1.28 Tb of high-quality non-human clean reads, corresponding to 49.24 million reads per child. Genes were identified by aligning the clean reads to the 9.9M human gut microbiome integrated gene catalog (IGC) [15]with an average of 80.1% reads in each sample being mapped (Additional file 1: Table S2). Screening the Maastricht Irritable Bowel Syndrome Cohort (MIBS-CO) [16] for healthy adult controls who fitted our inclusion criteria provided 62 metagenomic datasets from healthy Dutch adults for our age-based comparison (Additional file 1: Table S3). We applied the same pipeline as used for processing the children’s samples to analyze the published adult samples, resulting in 26.9 million clean reads per adult and a 77.4% IGC mapping ratio per adult.

Fig. 1 Comparison between early school-age Dutch children and adults. a Categories of phenotypic data collected within the KOALA cohort. b Box plot showing the gene-based α-diversity (Shannon index) in early school-age Dutch children (n = 281) and healthy Dutch adults (n = 62). c Genus-based principal component analysis (PCA) of children and adults. d Box plots of intra- and inter-group beta diversity based on genus profiles in children and adults. The “Intra-children” and “Intra-adults” indicate the genus-based beta diversity in children and adults, and the “Inter-groups” indicates the genus-based beta diversity between children and adults (***P < 0.001; Wilcoxon rank-sum test). e Box plots displaying the ten most abundant genera among children and adults. Genera indicated with red font are enriched in children, and genera in blue are enriched in adults. Boxes are ordered according to median relative abundance of genus in children Full size image

Comparison of gut microbiota in Dutch children and adults, and overweight and lean children

Comparison of gut microbial gene diversity in Dutch adults and children revealed an adult-like alpha diversity in the school-age children (Wilcoxon rank-sum test, P > 0.05, Fig. 1b). Principal component analysis (PCA) based on genus profiles showed no separation between Dutch children and adults (Fig. 1c). Within- and between-group β diversity of genus profiles further revealed a slightly more similar microbial composition within children compared to within adults (Fig. 1d). In line with previous studies comparing healthy infants, children (1–12 years), and adults [3, 17, 18], Dutch children showed a relative enrichment in abundance of the genus Bifidobacterium compared to adults (adjusted P < 0.05, Fig. 1e, Additional file 1: Table S4). Of note, the abundance of Bifidobacterium adolescentis, a reported adult-type Bifidobacterium species with no or poor ability for human milk oligosaccharides (HMOs) degradation [9, 19, 20], showed no differences between children and adults (Wilcoxon rank-sum test, P > 0.05, Additional file 1: Table S5). Comparisons further revealed that children were enriched in bacteria from the phylum Bacteroidetes including Bacteroides and Prevotella, while Firmicutes assigned to the genera Eubacterium, Clostridium, Dorea, and Coprococcus were more abundant in adults (Additional file 2: Figure S1a) (adjusted P < 0.05, Fig. 1e). A recent large gut microbiome cohort study on 1135 Dutch adult participants also reported higher abundances of Firmicutes (63.7%) than Bacteroidetes (8.1%) [21]. By contrast, findings from a USA cohort revealed greater abundances of Firmicutes and significantly lower abundances of Bacteroidetes in healthy children aged 7–12 years than in healthy adults [17], indicating that the composition and development of bacterial communities varies in populations with different geographic and genetic origins.

At the functional level, a total of 6771 KEGG (Kyoto Encyclopedia of Genes and Genomes) orthologues (KOs) were identified in the childhood samples (median of 4695 per individual) and assigned to eggNOG (evolutionary genealogy of genes: Non-supervised Orthologous Groups) functional categories. An adult-like stable composition based on 25 eggNOG functional categories was observed in children (Additional file 2: Figure S1b).

In order to determine if the previously reported relation between bacterial gene richness and BMI in adults [22] was also observed in children, we examined whether bacterial gene richness differed between overweight children (BMI z-score ≥ 1.04, n = 23) and lean children (BMI z-score < 1.04, n = 258). A BMI z-score ≥ 1.04 was used as the threshold for identification of overweight in children as recommended by the Dutch National Growth Study [23]. Interestingly, we observed a bimodal distribution of bacterial gene counts in the overweight group (Additional file 3: Figure S2) where the children with a gene number < 600,000 (n = 8) exhibited significantly higher BMI z-score (Wilcoxon rank-sum test, P = 0.016). By contrast, no such gene distribution pattern was observed in lean children. We did not identify significant differences in relative species abundances between overweight and lean children (Wilcoxon rank-sum test, adjusted P > 0.05).

Stratification of Dutch children based on their gut microbiome

Although, the concept of gut microbiota enterotyping has been highly debated [24], a consensus concerning enterotypes, including guidelines for rational enterotyping, was recently achieved [25]. Previous studies have suggested that human adult fecal metagenomes can be stratified into enterotypes that are independent of nationality, health status, age, BMI, or sex [26, 27], but associate strongly with long-term dietary patterns [28]. To examine this phenomenon in children, we here conducted a Dirichlet multinomial mixtures (DMM)-based enterotype analysis [29] to investigate the presence and characteristics of enterotypes in children based on the underlying compositional structure of their gut microbiome (see details in “Methods” section).

Three enterotypes were identified and found to be driven by a relatively high abundance of the genera Bacteroides (enterotype 1 (E1), n = 143), Prevotella (enterotype 2 (E2), n = 74), and Bifidobacterium (enterotype 3 (E3), n = 64), respectively (Fig. 2a, b, Additional file 1: Table S6). Surprisingly, the abundances of 50 out of 81 detected genera and 132 out of 214 detected species (species with ≥ 100 detected genes in any enterotype) differed significantly between enterotypes (Kruskal-Wallis test, adjusted P < 0.05, Additional file 1: Tables S7, S8). All Bifidobacterium spp. exhibited higher abundances in enterotype 3 (E3), with the adult-type species B. adolescentis and B. catenulatum dominating this genus (58.8% of the genus) (Additional file 1: Table S8). SparCC analysis revealed strong positive correlations between species enriched within enterotypes while negative correlations were found between species enriched in different enterotypes (Fig. 2c). For instance, the relative abundance of Prevotella copri (E2) negatively correlated with both the abundance of Bacteroides uniformis (E1) and Bifidobacterium longum (E3) (SparCC, pseudo P < 0.01, Fig. 2c).

Fig. 2 Stratification of early school-age children into three enterotypes based on their gut microbiome. a Scatter plot representing the three enterotypes identified using Dirichlet multinomial mixtures (DMM)-based clustering among early school-age Dutch children. MDS multidimensional scaling. b Genus abundance box plots showing the main contributors of each enterotype. c Correlations between enterotype enriched species, with the log10-transformed relative abundance of each species indicated by the circle area. Only the top 10 most abundantly enriched species in each enterotype are displayed. Red line indicates positive correlation and gray line indicates negative correlation (SparCC, pseudo P < 0.01) Full size image

Comparative analyses further revealed distinct compositional and functional differences of the gut microbiota between the three enterotypes. Interestingly, children with a gut microbiota belonging to E1 and E2 showed similar gene count and diversity (Dunn’s post-hoc test, P > 0.05) while children with E3 harbored about 110,000 (15%) fewer genes (median value: 604,692 genes/child, Dunn’s post-hoc test, P < 0.05) compared to the children with the other two enterotypes (Additional file 1: Table S6, Additional file 3: Figure S2a). Kruskal-Wallis tests on differences between the numbers of genes in each genus further demonstrated that E3 harbored significantly fewer genes from unannotated bacteria and several abundant genera including Bacteroides, butyrate-producing Faecalibacterium, Eubacterium, and Roseburia (Additional file 1: Table S9). E3 also exhibited lower Shannon diversity (Dunn’s post-hoc test, P < 0.05, Additional file 4: Figure S3b), and higher reads mapping rates (median of 64.2%) to taxonomically annotated genes than the other two enterotypes (Dunn’s post-hoc test, P < 0.05, Additional file 4: Figure S3c). At the functional level, E3 showed a slightly lower number of KOs than E1 (Dunn’s post-hoc test, P < 0.05, Additional file 4: Figure S3d). However, compared to the other enterotypes, E3 displayed similar functional diversity (Dunn’s post-hoc test, P > 0.05, Additional file 4: Figure S3e) and higher reads mapping ratios to the KO annotated genes (Dunn’s post-hoc test, P < 0.05, Additional file 4: Figure S3f). Lower β diversity was identified between E1 and adults, as compared to that of adults vs. E2 or E3 (Kruskal-Wallis test, P < 0.05, Additional file 4: Figure S3g), overall suggesting that the gut microbiota development may continue beyond 6–9 years of age. Altogether, our data points to a more adult-like gut microbiome in children belonging to E1 and much simpler structured gut microbiome in children belonging to E3.

KEGG pathway enrichment analysis revealed distinct differences in microbial functional patterns between enterotypes. Metabolic modules involved in biosynthetic pathways for leucine, lysine, serine, methionine, proline, and histidine production, together with amino acid transport systems for glutamate, branched-chain amino acid (BCAA), and d-methionine, were highly enriched in the microbiome of E3 (Reporter score > 1.96, Additional file 5: Figure S4a). Conversely, the module for cysteine biosynthesis was less represented in E3 (Reporter score < − 1.96, Additional file 5: Figure S4a), in agreement with prior studies reporting the lack of cysteine synthetase in the genomes of Bifidobacterium species [30, 31]. Furthermore, compared to the other two enterotypes, the gut microbiota of children belonging to E3 showed higher enrichment of functions involved in metabolism of simple sugars including glycolysis and the pentose-phosphate pathway, while functions for utilizing complex carbohydrates such as pectin, uronic acids, and glycosaminoglycan degradation were depleted (Additional file 1: Table S10, Additional file 5: Figure S4a). Consistently, Bacteroides and Prevotella have been reported to possess large numbers of genes for fermentation and utilization of complex polysaccharides [32]. In contrast, in vitro culture experiments reported that Bifidobacterium strains grew well on glucose or ribose containing media but exhibited poor or no growth on media containing non-HMO-derived complex carbohydrates such as inulin or exopolysaccharide [33, 34]. Besides these differences, we also observed enterotype-dependent differences in the potential for biosynthesis of several water and lipid-soluble vitamins. For instance, children with a gut microbiota belonging to E1 showed higher potential for biosynthesis of cobalamin (B12) and biotin (B7), whereas children belonging to E2 exhibited higher potential for biosynthesis of menaquinone (vitamin K), pantothenate (B5), and riboflavin (B2) (Additional file 5: Figure S4a). A list of KEGG modules that differed significantly in abundance between enterotypes is provided in Additional file 1: Table S10. A heatmap showed that the relative abundance profiles of eight selected KOs, which are responsible for key steps in amino acid biosynthesis and carbohydrate metabolism, allowed to distinguish children in E3 from the two other groups (Additional file 5: Figure S4b).

Associations between early events and pre-school lifestyle and gut microbiota in school-age children

In order to identify interactions between early environmental factors and the microbiota, we first conducted a PCA to assess the multivariate variation in children of the early events, pre-school dietary, and non-dietary lifestyle factors (Additional file 1: Table S1). We found that breastfeeding duration, educational level of mother at childbirth, and pre-school dietary patterns including intake of protein, fiber, and milk products contributed most to the variability in PC1 (15.05%, Fig. 3a), and total intake of carbohydrates and fat represented the second most important variation among children, as displayed in PC2 (12.74%, Fig. 3a). Interestingly, children in E3 exhibited lower PC1 scores but higher PC2 scores than children in E1 (Fig. 3a, Wilcoxon rank-sum test, P < 0.05). This inter-enterotype difference was governed by specific major contributors of the PC1 scores, including shorter breastfeeding duration and less intake of dietary fiber and plant-based protein in E3 as compared to the two other enterotypes (Kruskal-Wallis test, P < 0.05) (Fig. 3b, Additional file 1: Table S11). Next, we performed permutational multivariate analysis of variance (PERMANOVA) to assess the interactions between early events, pre-school lifestyle, and microbial gene profiles among all individuals and then went on to evaluate their interactions in each enterotype. Based on the entire population of children, early events including breastfeeding duration, the pre-school lifestyle including intake of plant-based protein and dietary fiber were significantly correlated with the microbial composition at 6–9 years of age (adjusted P < 0.05, Bray-Curtis distance, Fig. 3c). In addition, we identified several different correlation patterns within each enterotype (Fig. 3c). For instance, plant-based protein intake was significantly correlated with the gut microbiota in E2 (adjusted P < 0.05), but not in E1 and E3 (adjusted P > 0.05).

Fig. 3 Multiple early events and pre-school lifestyle associated with the school-age gut microbiota. a PCA showing the multivariate variation of children and the major contributions of different factors to PC1 and PC2. A total of 18 factors including early events and pre-school lifestyle (Additional file 1: Table S1) were subjected to PCA, and those factors with component scores for PC1 or PC2 ≥ 0.2 were shown as major contributors. Box plots showing the overall distribution of PC1 and PC2 scores within each enterotype (#P<0.05; Wilcoxon rank-sum test). b Significant major contributors in PCA between enterotypes (#P < 0.05, Wilcoxon rank-sum test; *P < 0.05, Dunn’s post-hoc test). The details of statistical results for all factors are shown in Additional file 1: Table S11. c PERMANOVA of the influence of single-factor early events and pre-school lifestyle on the gut microbial gene profile in the entire cohort and within each enterotype (#P < 0.05; * adjusted P < 0.05) Full size image

Spearman’s rank correlation analyses further revealed significant correlations between early events, pre-school lifestyle, and the relative abundance of certain microbial species (Additional file 6: Figure S5a, adjusted P < 0.05). Interestingly, bacteria enriched in E3 such as B. adolescentis, B. angulatum, and B. breve were negatively correlated with breastfeeding duration (Additional file 6: Figure S5a). The relative abundances of several Streptococcus species, reported to be part of the oral cavity microbiota of children [35,36,37], were positively related to maternal BMI (Additional file 6: Figure S5a, adjusted P < 0.05). The relative abundances of B. angulatum, B. dentium, and Streptococcus mitis were negatively correlated with plant-based protein intake (adjusted P < 0.05). Moreover, several species showed correlations with both early events and pre-school dietary patterns before age 5. For instance, Streptococcus parasanguinis and Streptococcus gordonii that associated with maternal BMI were also positively correlated with animal-based protein intake (adjusted P < 0.05). Similarly, the relative abundances of Collinsella intestinalis showed a positive correlation with maternal BMI and a negative correlation with dietary fiber intake (adjusted P < 0.05). Collectively, our findings indicate that both early events and pre-school dietary lifestyle contribute to shaping of the gut microbiota in school-age Dutch children, with different factors influencing each enterotype.

Enterotype-dependent associations between school-age metabolic phenotypes and pre-school lifestyle

An increasing number of studies has linked the gut microbiota with host metabolic phenotypes, including glucose and insulin homeostasis, and amino acid metabolism [38,39,40]. We examined the association between the gut microbiota and metabolic phenotypes including glucose, insulin, and amino acids levels measured in blood samples collected 3.5 h after the last meal on the same day as the fecal samples were collected (Additional file 6: Figure S5b). We observed that the relative abundances of Bacteroides spp., including Bacteroides xylanisolvens, B. dorei, B. vulgatus, and B. eggerthii, enriched in E1, were negatively correlated with plasma branched-chain amino acid (BCAA) levels. It has been demonstrated in a single mouse study that certain Bacteroides strains may contribute to BCAA degradation and hence reduce the circulating levels of BCAA in the host [41]. In addition, we also observed positive correlations between Dorea longicatena, Coprococcus comes, and plasma glutamate levels, which is in line with a recent study reporting an enrichment of these species and plasma glutamate levels in young obese Chinese adults [42]. The abundance of S. gordonii positively correlated with plasma total triglyceride and glucose.

Next, we examined if the measured levels of blood metabolic parameters, such as insulin and glucose, in school-age children differed between the three enterotypes. No differences were observed between individuals within the three enterotypes (Additional file 1: Table S11, Kruskal-Wallis test, P > 0.05).

Given the multi-relationships between gut microbiota and the above-mentioned early events and pre-school lifestyle, we next asked whether these factors would impact metabolic responses in the early school-age children. We first conducted Spearman’s rank correlation analyses between early factors and school-age metabolic phenotypes across the entire children cohort. We found that insulin levels exhibited negative associations with the pre-school lifestyle related to intake of plant-based protein and dietary fiber (Spearman’s correlation, adjusted P < 0.05, Additional file 7: Figure S6a). However, by stratifying according to enterotypes, we observed different correlation patterns among children within the different enterotypes. For instance, only children within E1 and E2 exhibited negative associations between insulin levels and intake of plant-based protein and dietary fiber (adjusted P < 0.05, Additional file 7: Figure S6b, c). In addition, the negative correlations between increased pre-school dietary fiber intake and reduced total serum triglyceride (TG) levels at school-age were solely present in E2 (Additional file 7: Figure S6c). Despite no such relationships were seen for E3 children, we found plasma free fatty acid (FFA) levels positively correlated with animal-based protein (Additional file 7: Figure S6d), and additionally, a few Streptococcus spp. to be positively correlated with the plasma FFAs levels only in this group (Additional file 8: Figure S7). Most of these enterotype-dependent associations persisted after adjusting for multiple covariates including sex, age, and BMI z-score at stool collection and early events by linear regression models (Additional file 1: Table S12).

Further, the children in E3 who did not exhibit a negative correlation between pre-school dietary lifestyle of plant-based protein and dietary fiber intake and blood insulin levels at school-age exhibited a lower potential for complex carbohydrate metabolism (Additional file 5: Figure S4, Additional file 7: Figure S6). We assumed that the different responses of insulin levels to dietary fiber intake might be due to variations in metabolites reported to influence systemic insulin levels such as butyrate [13] generated by the gut microbiota in response to intake and fermentation of complex carbohydrates. Hence, we compared the abundances of genes involved in the conversion of crotonyl-CoA to butyrate, which is one of the final steps for butyrate production. Interestingly, the abundance of the bcd gene encoding butyryl-CoA dehydrogenase (K00248) was significantly enriched in both E1 and E2 as compared to E3. Moreover, the abundances of ptb (phosphate butyryltransferase, K00634), buk (butyrate kinase, K00929), and but (butyryl CoA: acetate CoA transferase, K01034) were all significantly enriched in E1 as compared to E2 and E3, suggesting a higher potential for butyrate production in E1 (Fig. 4a, b). A taxa distribution analysis revealed that Faecalibacterium prausnitzii, Eubacterium halii, Roseburia inulinivorans, and Odoribacter splanchnicus largely accounted for the prevalence of the bcd genes when focusing on differences of butyrate production potentials between enterotypes (Fig. 4d, Additional file 1: Table S13). We further compared the relative abundance of genes required for succinate production between the three enterotypes, since succinate has been reported as a microbial product produced in response to dietary fiber intake that may contribute to improved plasma glucose and body weight regulation [43]. As shown, E2 exhibited significantly higher abundances of genes encoding the succinate dehydrogenase complex (K00239, K00240, and K00241) than E1, and the dominant contributor in E2 was Prevotella copri (Fig. 4c, e). Despite the fact that propionate is commonly converted from the succinate pathway in the intestine [44], no differences were observed between the three enterotypes with respect to the relative abundance of pct genes (K01026) encoding the propionate CoA-transferase, which is responsible for the last step in propionate production (Additional file 1: Table S14, Kruskal-Wallis test, P > 0.5).

Fig. 4 Enterotype-associated differences in potential for butyrate, succinate, and propionate production. a Pathway showing the genes involved in final biosynthetic steps from crotonyl-CoA to butyrate, including bcd (butyryl-CoA dehydrogenase, K00248), ptb genes (phosphate butyryltransferase, K00634), buk genes (butyrate kinase, K00929), and but (butyryl CoA:acetate CoA transferase, K01034). b The relative abundance of genes involved in butyrate production within each enterotype. c The relative abundance of genes involved in succinate production (succinate dehydrogenase complex (K00239, K00240 and K00241) within each enterotype. d Mean relative abundance of bcd genes (K00248) listed according to annotated bacterial species within each enterotypes. e Mean relative abundance of sdhA genes (K00239) listed according to annotated bacterial species within each enterotypes. Dunn’s post-hoc test; **P < 0.01, ***P < 0.001 Full size image

Collectively, our results illustrate that early events and pre-school dietary lifestyle impact on the development of the gut microbiota, which in turn seems to affect host metabolic responses.