A longitudinal study of the gut microbiota during the first three years of life: Links with problem behavior and executive functions at preschool age

Abstract

Early life is a sensitive period when microbiota-gut-brain interactions may have important impact on development. This study investigated the associations of the gut microbiota in the first three years of life (two, six, and 12 weeks, and one and three years) with problem behavior and executive functions in N = 64 three-year-old children. Higher relative abundance of Streptococcus at the age of two weeks, as well as its trajectory over time (including ages two, six and 12 weeks, and one and three years), was related to worse executive functions. Higher relative abundance of [Ruminococcus] torques group at the age of three years, as well as its trajectory from one to three years, was associated with less internalizing behavior. Besides, several robust age-specific associations were identified: higher Bifidobacterium relative abundance (age three years) was associated with more internalizing and externalizing issues; higher Blautia relative abundance (age three years) was linked to less internalizing behavior; and increased relative abundance of an unidentified Enterobacteriaceae genus (age two weeks) was related to more externalizing behavior. Our findings provide important longitudinal evidence that early-life gut microbiota may be linked to behavioral and cognitive development in low-risk children.

Introduction

The human gut harbors a great number of microorganisms, of which bacteria are an essential part. These microorganisms are collectively termed the ‘gut microbiota’ (Thursby & Juge, 2017). Not only has the gut microbiota been involved in many health outcomes, such as obesity, type 2 diabetes, and irritable bowel syndrome (Vos et al., 2022), but it has also been linked to mental health (Cryan et al., 2019). Accumulating evidence from animal and adult human studies has uncovered several key bidirectional communication pathways between the gut microbiota and brain functioning, named the microbiota-gut-brain axis (MGBA) (Cryan et al., 2019). Remarkably, the MGBA is not only functional in adults but starts playing an equally or even more important role at early ages with regard to child behavior and cognition (Cryan et al., 2019). Both the gut microbiota and the brain develop at a rapid pace during early life, however, only few studies investigated associations between the gut microbiota and behavior in such sensitive periods. Therefore, this study aimed to investigate the relations of the gut microbiota in the first three years of life with child problem behavior and executive functions at the age of three years.

The bidirectional interactions of the MGBA occur through intricately innervated and highly adaptable neuronal pathways, and extremely delicate and difficult-to-measure molecular communication systems (Cryan et al., 2019; de Weerth, 2017). For instance, short-chain fatty acids (SCFAs), mainly produced through dietary fiber fermentation by the gut microbiota, likely affect the brain via the vagus nerve, immunity, and the endocrine system (Dalile et al., 2019). Furthermore, specific microbial taxa can generate γ-aminobutyric acid (GABA), which is the main inhibitory neurotransmitter of the central nervous system and regulates many physiological functions (Mazzoli & Pessione, 2016; Silva et al., 2020). The symporter that mediates the uptake of microbiota-derived GABA is present through the gastrointestinal tract, suggesting that luminal GABA is able to cross the gut barrier and influence extra-gut targets. Although remaining controversial, recent studies suggest the permeability of the blood-brain barrier to GABA, implying its direct impact on the central nervous system (Mazzoli & Pessione, 2016). Besides, GABA receptors are widely expressed in enteric neurons and immune cells, indicating the role of GABA in regulating the gut-to-brain signaling and neuroinflammation (Auteri et al., 2015; Hyland & Cryan, 2010). Such pathways along the MGBA may partially explain how the gut microbiota impacts mental health.

The colonization of the gut by microorganisms mostly commences soon after birth and continues in the following years. The general consensus is that the gut microbiota develops into an adult-like configuration around the age of three years (Derrien et al., 2019), while some studies suggested that this step toward maturation may take longer than previously thought (Ou et al., 2022, 2023). Gut microbial disturbances during the early dynamic and sensitive colonization period can result in subsequent health problems, such as developing allergies and obesity (Zhuang et al., 2019). This is explained by the early-life programing theory, which refers to long-lasting changes and disruptions as a consequence of environmental exposures at a young age (Tarry-Adkins & Ozanne, 2011). In early life, the brain experiences numerous quick developments in neuronal proliferation, migration, differentiation, synaptogenesis, myelination, and apoptosis (Rice & Barone, 2000), largely impacting brain functioning, cognition, and behavior (Erus et al., 2015). Simultaneously, the microbiota is becoming established in the gut of infants and young children (de Weerth, 2017; Wang et al., 2018). Thus, alterations of the gut microbiota in early life may exert considerable effects on the development of the brain. Indeed, there is compelling evidence from animal studies supporting such a hypothesis (Clarke et al., 2014; Leclercq et al., 2017; O’Mahony et al., 2014; Stilling et al., 2015). This marks early life as a sensitive time window to obtain and maintain microbiota composition that will promote normal physical and mental development. However, we know little about early-life gut microbiota in association with child behavior and cognition. Specifically, how the gut microbiota and brain functioning, in particular host behavior, are interconnected in low-risk community infants and children (i.e., generally healthy and neurotypically developing) is underexplored. Knowledge of these associations, particularly when uncovered by comprehensive longitudinal studies, can provide insight into the typical early development of the gut microbiota in relation to child behavior and cognition.

First studies have found evidence for associations between the gut microbiota and child behavior and cognition. Regarding behavior, Loughman et al., reported that increased relative abundances of taxa belonging to the genus Prevotella at one year of age were associated with less internalizing behavior at age two (i.e., problem behavior affecting internal psychological conditions, characterized by withdrawal, anxiety, and emotional problems (T. M. Achenbach, 1966)) (Loughman et al., 2020). In our previous study, we found that the rise of Prevotella 9 in middle childhood was related to more externalizing behavior at age ten (i.e., problem behavior exhibited in the external environment, including features like impulsivity, aggression, and hyperactivity (T. M. Achenbach, 1966)) (Ou et al., 2022). Besides, Laue et al., observed a negative relation between Streptococcus peroris and internalizing behavior in girls before school age, and a positive association between Lachnospiraceae species and externalizing behavior in both genders (Laue et al., 2021). Furthermore, Lachnospiraceae species and Veillonella were linked to more internalizing behavior in preschoolers; interestingly, Veillonella was positively related to externalizing behavior as well (Van De Wouw et al., 2022). Additionally, increased alpha diversity was observed in preschool children with less internalizing (Laue et al., 2021; Van De Wouw et al., 2022).

Four other studies have found an underlying link between infant gut microbiota and child cognition (Aatsinki et al., 2020; Carlson et al., 2017; Rothenberg et al., 2021; Streit et al., 2021; Tamana et al., 2021). Cognition is fundamental for the development of executive functions, including higher-level cognitive processes like inhibitory control (Diamond, 2013). Specifically, a cross-sectional study found more Enterobacteriaceae species in relation to worse cognition at age 45 months (Streit et al., 2021). Longitudinal research reflected that high relative abundances of Bacteroides at age one year were related to better cognition at age two (Carlson et al., 2017; Tamana et al., 2021). Furthermore, Faecalibacterium at one year of age was associated with reduced cognitive functions at age two (Carlson et al., 2017). Additionally, a lower relative abundance of Bifidobacterium and a higher relative abundance of Clostridium at two-and-a-half months were linked to increased attention at eight months (Aatsinki et al., 2020). Moreover, Rothenberg et al., found that children with better cognition showed enriched Faecalibacterium, Sutterella, and Clostridium cluster XIVa at age three years. Finally, high alpha diversity at age one year was reported in two-year-old children with worse cognition (Carlson et al., 2017).

To conclude, a number of associations have been observed between the gut microbiota and problem behavior and cognition in early life, but findings are variable and inconsistent across studies, mainly due to different methodologies used regarding microbiota analyses, genomics, epidemiology, and statistics. Furthermore, most of the previous studies have assessed problem behavior and cognition by using only one questionnaire of a single reporter. In the current longitudinal study in a community sample of children, we investigated the gut microbiota in relation to problem behavior (i.e., internalizing and externalizing behavior) and executive functions (i.e., advanced cognitive abilities, including inhibitory control (Diamond, 2013)) using questionnaires of multiple reporters and behavioral tasks. We had the following two hypotheses: (1) relative abundances and alpha diversity (i.e., Chao1, Shannon, and phylogenetic diversity) of the gut microbiota at age three years are associated with reported problem behavior and executive functions at the same age; (2) relative abundances and alpha diversity of the gut microbiota at early ages (i.e., two, six, and 12 weeks, and one year) are associated with reported problem behavior and executive functions at age three.

We investigated these hypotheses in three ways: (1) as the gut microbiota is highly dynamic in early life, its composition at different ages may be differently associated with problem behavior and executive functions later in life. For this reason, we analyzed the overall gut microbiota composition in relation to preschool-aged cognitive measures in an age-specific manner; (2) for the same reason, relations regarding a single taxon and an alpha diversity index were analyzed in an age-specific manner; (3) based on the age-specific analyses, we explored the trajectories of taxa and alpha diversity parameters in association with mental outcomes over the whole study period. Figure 1 shows the workflow of our analyses. Considering that most published findings were at the genus level, we performed our analyses at the same taxonomic level. However, given the paucity of studies on these relations at such early ages, we did not hypothesize specific associations between microbial taxa and mental outcomes.

Figure 1. Workflow of the analyses.

Materials and methods

Participants

The current study is part of the longitudinal Dutch study named BINGO where early factors associated with child development were investigated. Participants were healthy children and their parents living in the Netherlands. Detailed in- and exclusion criteria are described in a previous publication (Hechler et al., 2018). At baseline, 88 pregnant women were recruited in the Arnhem-Nijmegen region for the BINGO study; 96% were born in the Netherlands. Postnatal exclusion criteria included: complications during pregnancy, gestational age at birth < 37 weeks, birth weight < 2500 g, 5-min Apgar score < 7, and congenital malformations. Seventy-seven mothers were followed up after postnatal exclusion. At three years of child age, 76 families were approached (one drop-out occurred during the previous measurement rounds). Among them, two families could not be contacted, six families did not participate due to time constraints, and one family dropped out due to personal reasons. Parental demographics did not differ significantly between participating and nonparticipating families. This resulted in a final sample of 67 families. Of them, 64 families participated in home visits when their children reached age three, and the other three families were unable to join home visits but filled out questionnaires in this assessment round. Both parents participated in 54 families (81%, 54/67), and only mothers participated in 13 families (19%, 13/67).

Ethics

The BINGO study was independently reviewed by the Ethics Committee of Social Sciences of Radboud University, and no formal objection to this research was made [ECSW2014–1003–189 and amendment: ECSW–2018–034]. The current study was preregistered on the Open Science Framework: https://osf.io/vwgef with amendment: https://osf.io/nyeb4.

Data collection procedure

Collection of child stool samples was done at the ages of two, six, and 12 weeks, and one and three years. Stool samples were stored in the participant’s freezer (−20°C) until they were collected with a portable freezer. The stool samples were stored at -80°C at Radboud University prior to being processed at Laboratory of Microbiology at Wageningen University & Research.

Home visits took place when the child turned three years old. Prior to the home visit, mothers and their partners independently filled in digital questionnaires about their child’s problem behavior and executive functions. During the home visit, the child performed inhibitory control tasks. Tasks were video recorded and afterward rated by two trained observers. Observers were trained by use of a coding manual specific to each task.

Measures

Gut microbiota composition

Stool samples were collected with a polystyrene 10 mL stool container. Total DNA was extracted from 0.01 to 0.15 g of stool sample with 300 µL of Stool Transport and Recovery Buffer by double bead-beating steps as previously described (Gu et al., 2018). The variable V4 region of prokaryotic 16S ribosomal RNA (rRNA) genes was then amplified by PCR in duplicate reactions, by using primers 515F-n (5’-GTGCCAGCMGCCGCGGTAA) and 806R-n (5’-GGACTACHVGGGTWTCTAAT) (Gu et al., 2018). The 16S rRNA gene sequencing was completed on the Illumina HiSeq platform by Eurofins Genomics Germany GmbH.

Behavioral measures

Parental questionnaires

Mothers and their partners filled in all questionnaires mentioned below. However, because fewer partners completed the questionnaires, we used partner reports for sensitivity analyses to validate the maternal reports by calculating Kendall correlations between both. The non-parametric Kendall method was chosen due to its better performance in handling non-normally distributed data and tied values (Kendall, 1945). Maternal reports were used as the final measure of reported problem behavior and executive functions.

To assess child problem behavior, the Child Behavior Checklist for ages of one and a half to five years (CBCL, 103 items) (Thomas M. Achenbach & Ruffle, 2000) and the Strengths and Difficulties questionnaire (SDQ, 25 items) (Goodman, 1997) were used. The CBCL and the SDQ include internalizing and externalizing subscales, consisting of items scored on a three-point Likert scale. The SDQ can detect problem behavior as accurately as the CBCL does (Goodman & Scott, 1999). Raw scores for both questionnaires were used as outcome measure in order to compare and possibly aggregate the measures. However, given that the Kendall correlations on the same subscales of the CBCL and the SDQ were lower than 0.5 (Table S1), we included both instruments separately in the analyses. In both instruments, higher scores on subscales indicate more problem behavior.

To evaluate child executive functions, the Behavior Rating Inventory of Executive Function – preschool version (BRIEF-P, 63 items) questionnaire for preschoolers (Sherman & Brooks, 2010) and the Ratings of Everyday Executive functions (REEF, 77 items) (Nilsen et al., 2017) were used. The BRIEF-P and the REEF are scored on three- and four-point scales, respectively. A higher score on the BRIEF-P indicates worse executive functions, while a higher score on the REEF indicates better executive functions. The BRIEF-P is a commonly used questionnaire that measures general executive functions and does not differentiate between different situations. The REEF rates executive functions in different situations (e.g., executive functions around friends, during grocery shopping, or in the community) and determines an average score. Raw scores for both questionnaires were used as outcome measure in order to compare and possibly aggregate the measures. However, Kendall correlations between the BRIEF-P and the REEF were lower than 0.5 (Table S1), hence both instruments were included in the analyses.

Parental questionnaires were considered acceptable and reliable based on their ω_total (ranging between 0.65 and 0.94) or Cronbach’s α values (ranging between 0.83 and 0.96) (Table S2) (Revelle & Condon, 2019).

Inhibitory control tasks

Six different behavioral tasks with good reliability (i.e., Flanker, Whisper, Gift Wrap, Gift Delay, Snack Delay, and Bear Dragon) were performed to measure inhibitory control as previously stated in detail (Willemsen et al., 2021). Observer reliability was determined by the Intraclass Correlation Coefficient (ICC) relying on absolute agreement. The ICC’s for the inhibitory control tasks ranged from r = 0.84 to r = 0.96 (p < 0.001). Snack Delay and Bear Dragon were excluded from the analyses due to insufficient variation and low number of children that passed the practice trials, respectively. The other four tasks were included in our study. Higher scores on these tasks indicate better inhibitory control.

Statistical analyses

Pre-processing of sequence data

Sequence data were processed via NG-Tax 2.0 with default settings (Poncheewin et al., 2020; Ramiro-Garcia et al., 2018), with SILVA SSU 16S rRNA gene reference database (version 132) (Quast et al., 2012). The raw amplicon sequence variant (ASV) count data were used to calculate alpha diversity by the ape (Paradis, 2020) and the picante (Kembel, 2020) packages. Then, ASV count data were glommed at the genus level prior to analyses.

Gut microbiota composition and development over the first three years of life

For descriptive purposes, we first delineated gut microbiota composition and development in the first three years after birth (including all samples at the age of two, six, and 12 weeks, and one and three years). We compared differences in alpha diversity indices, including Chao1, Shannon, and phylogenetic diversity, between ages using Wilcoxon rank-sum tests corrected with the False Discovery Rate (FDR) method. Next, we also compared beta diversity between ages by conducting Principal Coordinate Analysis (PCoA) via the vegan package (Oksanen et al., 2020). Considering that PCoA can be applied to different dissimilarity and distance metrics that all differ in specific aspects and corresponding interpretation, we included the Bray-Curtis, weighted Jaccard (formula = 2*Bray-Curtis dissimilarity / (1 + Bray-Curtis dissimilarity)), unweighted UniFrac, weighted UniFrac, and Aitchison (the Euclidean distance based on centered-log-transformed ASV count data) methods, to comprehensively describe the compositional differences. Except for the Aitchison distance, we transformed genus-level count data into relative abundances before calculating other dissimilarity and distance metrices. Significance was determined as a p value lower than 0.05 for non-multiple tests and an FDR-adjusted p value lower than 0.05 for multiple tests.

Additionally, we visualized average and individual relative abundances at the genus level over the study period by using a bar plot and a heatmap, respectively. To identify differentially abundant microbial taxa at the genus level between ages, we conducted the Linear Discriminant Analysis Effect Size (LEfSe) method by using the microbiomeMarker R package (Segata et al., 2011), with a log-transformed Linear Discriminant Analysis (LDA) score higher than two indicating significance.

Confounding effects

In our original preregistration, we considered child age and diet quality as potential confounders (i.e., variables that influence both the independent variables and the outcome). After reconsideration, both variables were removed as potential confounders due to two major reasons (amendment can be found via https://osf.io/nyeb4): (1) low variation in child age (see Figure 2 for notes regarding ages); (2) our previous study using the same cohort found no significant associations of diet quality with behavior and executive functions (Willemsen et al., 2021). Given these considerations, no confounders were accounted for in the models performed in this study. Note that potential covariates of the independent variables only (i.e., the gut microbiota) were not accounted for in downstream analyses (Cinelli et al., 2020), as they would remove variation in the gut microbiota data, which was not the purpose of this study. These potential confounders and covariates as well as their relations to the gut microbiota and behavioral outcomes are displayed in a directed acyclic graph (Figure S1).

Figure 2. Alpha and beta diversity of the gut microbiota in the first three years of life. (a–c) alpha diversity as measured by Chao1, Shannon, and phylogenetic diversity indices. Wilcoxon rank-sum tests were conducted between ages and corrected with the FDR method (ns, not significant; *, <0.01). Age2w_mean±sd = 2.08 ± 0.28. Age6w_mean±sd = 6.23 ± 0.55. Age12w_mean±sd = 12.27 ± 0.42. Age1y_mean±sd = 1.04 ± 0.08. Age3y_mean±sd = 3.18 ± 0.10. (d–h) principal coordinate plots of beta diversity, based on different pairwise dissimilarity (Bray-curtis and weighted Jaccard) and distance (UniFrac and Aitchison) matrices, with points and ellipses colored by ages (Lake blue, two weeks; orange, six weeks; purple, 12 weeks; Pink, one year; grass green, three years).

Data imputation and transformation

Missing values (proportion of missing values is shown in Table S3) in problem behavior, executive functions, and inhibitory control were imputed ten times together, by using the predictive mean match method in the R package mice (Buuren, 2021). The imputation model was conducted separately at each age. For instance, at the age of three, 64 children provided gut microbiota data, and their missing values in aforementioned behavioral measures were imputed jointly in one model. No auxiliary variables (i.e., variables that are not included in analyses, but are correlated with imputed variables) were considered in the imputation.

For both random forest models and the Bayesian linear regression models, genus-level relative abundance data were used. Numeric variables were standardized (i.e., subtracting the mean and dividing by the standard deviation) for the Bayesian models only, as random forest models rely on decision trees for which standardization is considered unnecessary.

Main analyses

To determine whether gut microbiota composition in the first three years of life (i.e., two, six, and 12 weeks, and one and three years) is associated with problem behavior (i.e., internalizing and externalizing behavior) and executive functions (including inhibitory control) at age three, we conducted random forest models and the Bayesian linear regression models (Bürkner, 2017; Kuhn et al., 2020). Random forest is first of all well suited to analyze microbiome data as it is appropriate for high-dimensional data, invariant to scaling of inputs, computationally efficient, and able to uncover nonlinear relationships (Belk et al., 2018; Louppe, 2014; Namkung, 2020). The first random forest model was applied to assess the contribution of the total gut microbiota composition on our behavioral outcomes. This was done for the purpose of exploring the gut microbial community as a whole to account for the complex interplay between taxa. The following random forest model was applied as a preselection tool, to select possibly important taxa from high-dimensional data, before passing them on to the Bayesian linear regression model. The Bayesian model was first used to determine age-specific relations (i.e., directions and strengths) of a selected taxon and an alpha diversity index with each outcome measure. By looking at the different time points separately, these analyses can help identify periods of development that are sensitive to certain microbial compositions. Although not preregistered, after reconsideration, we decided to perform an additional analysis to optimize the use of our longitudinal data. Based on the age-specific results, we implemented a multilevel Bayesian model to determine whether trajectories of change in the gut microbiota were associated with the outcome measures at age three. Figure 1 shows the workflow of our analyses.

Age-specific analyses - Determining the contribution of the overall gut microbiota to each behavioral measure through random forest models

Data were imputed ten times: data were randomly split into a training dataset (including 80% participants) and a testing dataset (including 20% participants), leading to ten training datasets and ten corresponding testing datasets. The procedure of data splitting was applied to children who provided gut microbiota information at each age separately. To prevent data leakage, the missing values of behavioral measures were imputed in training and testing datasets separately (ten times) as described earlier. Next, we included genus-level relative abundances of overall gut microbiota as independent variables and one behavioral measure as an outcome. This step was performed on each individual behavioral measure separately. To train the model, a ten-repeated ten-fold cross-validation was conducted on each complete training dataset including imputed values via the caret package (Kuhn et al., 2020). Afterward, we used the trained model to obtain predicted behavioral outcomes of each corresponding complete testing dataset including imputed values. Similarity between predicted and actual behavioral outcomes of the complete testing dataset was measured by the Pearson correlation with its p value obtained from a permutation test (N = 1000). Considering that data splitting and imputation resulted in multiple datasets, we used the median value of the Pearson correlation coefficient from multiple cases to represent the final similarity. The p value corresponding to this median was included. P values were adjusted with FDR methods, with corrected values under 0.05 indicating significance.

Age-specific analyses - Preselecting potentially important gut microbiota contributing to each behavioral measure through random forest models

To identify microbial taxa that contribute to each behavioral outcome, we measured the change in the generalized cross-validation (GCV) value in the random forest model. Larger GCV changes indicate more contribution of the independent variable to the model, in other words, this analysis shows which taxa are potentially more important (Kuhn et al., 2020). Unlike the first random forest model, we did not split the data but used the whole dataset here, because we prioritized the structure of the model and a large sample size can provide more valid information. Missing values of behavioral measures in the whole dataset were imputed as described in the section on data imputation and transformation. Then, relative abundances of all taxa were treated as independent variables with one behavioral measure as an outcome. This procedure was performed on each behavioral measure separately. Next, we carried out a ten-repeated ten-fold cross-validation on each complete dataset containing imputed data and calculated average GCV values of multiple datasets acquired from data imputation. Based on the size of average GCV values, we selected the largest 20 taxa as the top 20 in importance. These 20 taxa were then passed to the Bayesian linear regression models to confirm their actual associations with behavioral measures.

Age-specific analyses – Associating the gut microbiota with behavioral measures by using Bayesian linear regression models

We implemented Bayesian linear regression models to estimate the relations of relative abundances of the selected top 20 microbial taxa with a prevalence value higher than 10% and alpha diversity with the child behavioral measures. Compared to standard linear regression models, the Bayesian linear regression models compute the probability of different effects rather than simply reporting single estimates of the “true effect” (Bürkner, 2018). We performed the Bayesian models by using the brms R package built based on the programing language Stan (Bürkner, 2017). The brm function within the brms package was used with the Gaussian distribution (mean = 0, std = 1) as the prior distribution for all beta coefficients and the Student’s t-distribution for error distribution (due to better performance in handling extreme values) (Lange et al., 1989). A list containing multiple complete datasets including imputed data was directly passed to the brm function, which in turn generated a single estimate. Other arguments of the brm function were set as follows: chains = 4, iter = 2000, and warmup = 1000. Under these settings, chains converged properly with Rhat values lower than 1.01. Regarding the outcomes of the Bayesian models, the less the posterior distribution overlaps with zero, the more likely a relation is positive or negative. In the current study, we defined a relation as positive or negative with confidence when its 95% credible interval (CI) excluded the value zero.

Trajectory analyses – Relating the developmental trajectories of the gut microbiota to behavioral measures through multilevel Bayesian linear regression models

To make maximum use of our longitudinal data, we conducted multilevel models to investigate relations between the developmental trajectories of the gut microbiota and behavioral measures. The multilevel models were performed on microbial taxa and alpha diversity with confident age-specific relations to behavioral measures (i.e., as determined by the Bayesian linear regressions described above). In the multilevel models, microbial and behavioral information as well as the actual age were level 1 variables, and the child was the level 2 variable. Note that missing values in behavioral measures and actual age were not imputed and that in these analyses we used the same distributions and arguments as described earlier. Before performing a testing model, we first checked the intraclass correlation (ICC) of an intercept-only model. When the 95% CI of an ICC excluded the value zero in the intercept-only model, multilevel strategies were used. A trajectory relation was considered with confidence when there was no overlap between its 95% CI and zero.

With respect to taxa, when their prevalence was higher than 10% at five-time points (i.e., two, six, and 12 weeks, and one and three years), multilevel models were performed on the pooled data of all ages. When only the first three time points met the 10% criteria, multilevel models were carried out by pooling samples at these three ages together. When only the prevalences at the last two ages were higher than 10%, multilevel approaches were done in the pooled aged-one-and-three years samples. Rhat values were used to check chain convergence.

Results

Population characteristics and descriptives

Demographic data and descriptives of study variables are shown in Table 1. Roughly 50% of the children were girls. Mothers were mostly highly educated (86.2%). Scores on the questionnaires measuring child problem behavior and executive functions did not differ significantly between mothers and their partners, and they were significantly positively intercorrelated (Table S1).

Table 1. Descriptives of study subjects

In the assessment round at age three, 64 children, 64 mothers, and 51 partners participated in the study. In total, 66, 70, 73, 72, and 64 fecal samples were collected at ages two, six, and 12 weeks, and one, and three years, respectively. Completion percentages were based on number of participating individuals (i.e., completion rates for mothers are based on 64 participating mothers, and completion rates for partners are based on 51 participating partners). P = partner; M = mother; CBCL = the child behavioral checklist; SDQ = the strengths and difficulties questionnaire; BRIEF-P = behavior rating inventory of executive functions – preschool; REEF = ratings of everyday executive functioning. Differences were compared between mother and partner reports by Wilcoxon tests. None of them were significant before or after FDR adjustments.

Gut microbiota composition and development over the first three years of life

We analyzed microbial composition of 345 fecal samples taken at five-time points. A total of 42,056,591 high-quality reads were obtained after being processed with NG-Tax 2.0. Within these reads, 220 microbial taxa were identified at the genus level mainly belonging to the phyla Firmicutes, Actinobacteria, Bacteroidetes, Proteobacteria, and Verrucomicrobia.

For descriptive purposes, we compared alpha and beta diversity between ages (Fig. 2; beta-diversity comparisons between the first three ages and between the last two ages are displayed in Figures S2 and S3) and delineated a general developmental trajectory of the gut microbiota over time (Fig. 3a). Diversity comparisons reflected profound compositional differences between infancy and preschool age. These differences were visualized by the heatmap showing individual relative abundance data (Figure S4). LEfSe identified a total of 106 differentially abundant microbial taxa between ages (log-transformed LDA scores higher than two; Table S4). Due to the large number of significant taxa, only the taxa with log-transformed LDA scores higher than four are highlighted and displayed in Figure 3b, such as an unidentified genus within Enterobacteriaceae, Lactobacillus, Bifidobacterium, Faecalibacterium, and Blautia.

Figure 3. Characteristics of the gut microbiota in the first three years of life. ( a ) Average relative abundances of the gut microbiota at the genus level over time. Others represent genera with relative abundances lower than 1%. ( b ) Differentially abundant genus-level taxa between ages, identified by linear discriminant analysis effect size (LEfSe) with log-transformed linear discriminant analysis (LDA) scores higher than four.

Age-specific analyses

Determining the contribution of the overall gut microbiota to each behavioral measure through random forest models

To explore whether the overall microbial composition in the first three years (i.e., at ages two, six, and 12 weeks, and one and three years) contributes to problem behavior and executive functions at age three, we compared the similarity between the actual and the predicted behavioral results. As shown in Table S5, 92% (46/50) of the models showed insignificant absolute correlation coefficients (i.e., lower than 0.3), indicating a low likelihood that the gut microbiota can contribute to behavioral outcomes. Regarding the 8% (4/50) models with correlation coefficients higher than 0.3, the similarity remained insignificant between the actual and predicted data, implying the same low likelihood. The random forest models showed that the overall gut microbiota did not contribute to problem behavior and executive functions in the present study.

Preselecting potentially important gut microbiota contributing to each behavioral measure through random forest models

As planned in our preregistration, we preselected the microbial taxa that may contribute to the behavioral outcomes the most (i.e., top important taxa) based on GCV values, by performing separate random forest models at each age. The top 20 important taxa at the genus level are depicted in Figures 4 and 5, with the following observations:

1. (1) Bacteroides and Clostridium sensu stricto 1 were the most frequent contributors to CBCL internalizing behavior;

2. (2) Bacteroides and Bifidobacterium were the most frequent contributors to CBCL externalizing behavior, SDQ internalizing and externalizing behavior, and BRIEF-P executive functions;

3. (3) Bacteroides and Blautia were the most frequent contributors to REEF executive functions;

4. (4) Additionally, Bacteroides and Bifidobacterium were the most frequent contributors to the behavioral measures of inhibitory control (i.e. Flanker, Whisper, Gift Wrap, and Gift Delay) (Fig. 6).

Figure 4. Heatmap showing the top 20 important microbial taxa over time and their associations to problem behavior at age three as reported by the mother. The top 20 important genus-level taxa within each age (i.e., 2w, two weeks; 6w, six weeks; 12w, 12 weeks; 1y, one year; 3y, three years) per behavioral measure are shown on the right side of the figure. Behavioral measures include: CBCL-I = internalizing behavior measured by the CBCL; SDQ-I = internalizing behavior measured by the SDQ; CBCL-E, externalizing behavior measured by the CBCL; SDQ-E, externalizing behavior measured by the SDQ. The orange scale indicates the importance of the taxa, with darker color referring to increased importance. The importance was determined by the generalized cross-validation value, with a larger value change indicating more contribution of a taxon to the model, i.e., which taxon is more important. As not all taxa appeared in the top 20 list at each time point, these absent taxa are colored in gray. Numbers on the left side of the figure show how many times a taxon appeared to be in the top 20 list of a behavioral measure over time. The frequently appearing taxa are bolded and colored in orange (five times), yellow (four times), or green (three times).

Figure 5. Heatmap showing the top 20 important microbial taxa over time and their associations to executive functions at age three as reported by the mother. The top 20 important genus-level taxa within each age (i.e., 2w, two weeks; 6w, six weeks; 12w, 12 weeks; 1y, one year; 3y, three years) per cognitive measure are shown on the right side of the figure. The measures include: BRIEF-P = executive functions measured by the BRIEF-P; REEF =, executive functions measured by the REEF. The orange scale indicates the importance of the taxa, with darker color referring to increased importance. The importance was determined by the generalized cross-validation value, with a larger value change indicating more contribution of a taxon to the model, i.e., which taxon is more important. As not all taxa appeared in the top 20 list at each time point, these absent taxa are colored in gray. Numbers on the left side of the figure show how many times a taxon appeared to be in the top 20 list of a measure over time. The frequently appearing taxa are bolded and colored in orange (five times), yellow (four times), or green (three times).

Figure 6. Heatmap showing the top 20 important microbial taxa over time and their associations to observed inhibitory control behavior at age three. The top 20 important genus-level taxa within each age (i.e., 2w, two weeks; 6w, six weeks; 12w, 12 weeks; 1y, one year; 3y, three years) per inhibitory control task are shown on the right side of the figure. The tasks include: FL = flanker; WHSP = whisper; GW = gift wrap; GD = gift delay. The orange scale indicates the importance of the taxa, with darker color referring to increased importance. The importance was determined by the generalized cross-validation value, with a larger value change indicating more contribution of a taxon to the model, i.e., which taxon is more important. As not all taxa appeared in the top 20 list at each time point, these absent taxa are colored in gray. Numbers on the left side of the figure show how many times a taxon appeared to be in the top 20 list of a task over time. The frequently appearing taxa are bolded and colored in orange (five times), yellow (four times), or green (three times).

Associating the gut microbiota with behavioral measures by using Bayesian linear regression models

To confirm whether the aforementioned top 20 important microbial taxa were associated with problem behavior and executive functions, we performed the Bayesian linear regression model on each genus-level taxa (relative abundance) and behavior pair. Table 2 shows the strongest observed associations of these pairs (i.e. estimates higher than 0.2 or lower than −0.2).

Table 2. Associations of the gut microbiota in the first three years of life with behavioral measures at age three

Associations of estimates with 95% credible intervals (CIs) excluding 0 are presented. Behavioral measures include: CBCL = the child behavioral checklist; SDQ = the strengths and difficulties questionnaire; BRIEF − P = behavior rating inventory of executive functions − preschool; REEF = ratings of everyday executive functioning.

Remarkably, there were several highly present taxa (i.e. prevalent in more than 80% of the samples, relative abundance higher than 10%) in relation to the outcome measures: Bifidobacterium at age three years was associated with more internalizing and externalizing behavior (est. = 0.27 for both), Blautia at three years was linked to less internalizing behavior (est. = −0.25), and an unidentified taxa within the Enterobacteriaceae family was related to more externalizing behavior (est. = 0.25).

Next, we checked for consensus between the questionnaires assessing the same construct. For internalizing behavior, there was no consensus between the associations found for the CBCL and the SDQ. For externalizing behavior, more Parabacteroides at two weeks were associated with less externalizing behavior in both the CBCL (est. = −0.30) and the SDQ (est. = −0.28). An opposite finding was found for Butyricicoccus at one year in relation to more externalizing behavior by the CBCL (est. = 0.23), while at three years, it was associated with less externalizing behavior by the SDQ (est. = −0.35).

Within the CBCL results, Barnesiella at age three years was associated with more internalizing (est. = 0.31) and externalizing behavior (est. = 0.33). Within the SDQ results, Bifidobacterium at age three years was associated with more internalizing and externalizing behavior (est. = 0.27 for both).

Regarding executive functions, Ruminococcus 2 at one year and [Ruminococcus] torques group at age three years, were associated with better executive functions as measured by the BRIEF-P (est. = −0.30, note that higher scores on the BRIEF-P indicate worse executive functions) and worse executive functions measured by the REEF (est. = −0.24), respectively. Lastly, Halomonas at six weeks was associated with worse executive functions as measured by the BRIEF-P (est. = 0.24) and the REEF (est. = −0.24).

Different associations were found for the Flanker and the Gift Delay tasks. For the Flanker, relations were identified at the age of six weeks, and one and three years, while for the Gift Wrap, associations were observed at age one year only. This may be due to a highly dynamic gut microbiota ecosystem in early life, of which composition at different ages may be variously linked to executive functions.

There were some overlapping associations between the questionnaires on problem behavior and executive functions. Parabacteroides at two weeks were associated with better executive functions (REEF, est. = 0.25), and less externalizing behavior (CBCL and SDQ, est. = −0.30 and −0.28, respectively). Another consistent result was Streptococcus at two weeks in relation to worse executive functions (BRIEF-P, est. = 0.40) and more externalizing behavior (CBCL, est. = 0.26).

We also measured behavioral relations to alpha diversity, including Chao 1, Shannon and phylogenetic diversity by using the Bayesian linear regression models (strongest results are displayed in Table 3). Interestingly, relations were only observed for alpha diversity at age two weeks. Higher Chao 1 values were associated with less internalizing behavior (CBCL) (est. = −0.28). Furthermore, Chao1 values were in positive relation to better executive function performance (REEF and Gift Wrap, est. = 0.31 and 0.43, respectively). Lastly, higher phylogenetic diversity at age two weeks was also linked to better inhibitory control during the Gift Wrap task (est. = 0.32).

Table 3. Associations of alpha diversity in the first three years of life with problem behavior, executive functioning, and inhibitory controls at age three

Associations of estimates with 95% credible intervals (CIs) excluding 0 are presented. CBCL = the child behavioral checklist; REEF = ratings of everyday executive functioning.

Trajectory analyses

Relating the developmental trajectories of the gut microbiota to behavioral measures through multilevel Bayesian linear regression models

Based on the results of age-specific Bayesian models and the 10% prevalence rule applied to microbial taxa (Table S6), we identified 16 pairs (including 12 pairs of taxa and behavioral measures, and four pairs of alpha diversity and behavioral measures) available at all five ages (i.e., two, six, and 12 weeks, and one and three years), three at the first three ages, and 12 at the last two ages (Table S7). Higher relative abundances of Streptococcus over the first three years of life were weakly related to worse executive functions reported by the BRIEF-P (est. = 0.05; higher scores on the BRIEF-P indicating worse performance), conforming to earlier age-specific findings. We also found that the trajectory of [Ruminococcus] torques group from age one to three was negatively related to internalizing behavior (SDQ, est. = −0.22), implying that higher relative abundances were associated with fewer internalizing difficulties during this period. No enduring associations were observed with confidence regarding alpha diversity.

Discussion

In this longitudinal study, we investigated associations of the gut microbiota during early life with problem behavior and executive functions, including inhibitory control, at three years of age. Several associations with behavior and cognition were found for relative abundances of microbial taxa and alpha diversity throughout the first three years of life. Table S8 provides an overview of the different microbiota taxa and microbial diversity index at different ages relative to the developmental findings. In addition, Table S8 shows the existing literature relative to our findings. Based on this table, we discuss the most prominent findings below.

We found evidence that increased relative abundance of Streptococcus, specifically at the age of two weeks and over the first three years after birth, was associated with worse executive functions at the age of three years. This result indicates that Streptococcus might affect cognitive development throughout early life. Relations between early-life relative abundances of Streptococcus and behavior and cognition in typically developing children have not been observed in previous literature. For comparative purposes, we, therefore, examined studies on microbiota composition in children diagnosed with neurodevelopmental disorders as they mostly have comorbid behavioral and cognitive issues (Schoemaker et al., 2014). According to a systematic review, children with ASD (autism spectrum disorder) frequently show an overgrowth of Streptococcus (Bundgaard-Nielsen et al., 2020). Although gut microbiota dysbiosis in ASD was seemingly partially attributed to an altered dietary pattern (Li et al., 2022; Welberg, 2022), diet was not correlated to the mental outcomes of our community samples. In addition to Streptococcus, both age-specific and trajectory relations were discerned for the [Ruminococcus] torques group: higher relative abundances at the age of three years and throughout the period from one to three years of age were associated with fewer child internalizing problems. Previous research showed excessive absolute abundances of fecal [Ruminococcus] torques group in children with ASD (Wang et al., 2013) and this taxon was strikingly increased in patients with inflammatory bowel disease (Png et al., 2010).

With respect to microbial taxa that only showed age-specific associations, we observed higher relative abundances of Bifidobacterium at age three years related to more internalizing and externalizing behavior at the same age. Interestingly, a systematic review showed decreased Bifidobacterium in ASD children compared to neurotypically developing controls (Xu et al., 2019). Besides, supplementing ASD children with a prebiotic galacto-oligosaccharide increased Bifidobacterium populations in the gut and alleviated autistic symptoms (Grimaldi et al., 2017). However, opposite roles of Bifidobacterium have been described in major depressive disorder (MDD) (Cheung et al., 2019; Knudsen et al., 2021). Such inconsistency also takes place within ADHD studies. Two studies reported that Bifidobacterium longum mitigated ADHD (Finegold et al., 2010; Pärtty et al., 2015), while another study found overgrowing Bifidobacterium species in ADHD subjects (Aarts et al., 2017). In addition, in the current study, a higher relative abundance of Blautia at the age of three years was related to fewer internalizing difficulties (as well as to better executive functions at the same age). Blautia is suggested to play an important role in nutrient absorption and digestion (Eren et al., 2015), and in child gut microbiota development towards a normal adult-like configuration (Hsiao et al., 2014). In line with our findings, depleted Blautia was seen in ASD populations aged from two to 18 years old by several studies as concluded in a systematic review (Liu et al., 2019). However, elevated levels of Blautia were reported in relation to MDD in adults (Cheung et al., 2019) and ADHD in three-year-old children (Laue et al., 2021), indicating that different mechanisms may be involved depending on the psychopathology and chronological age. Lastly, we observed a positive relation between one unidentified Enterobacteriaceae genus at the age of two weeks and externalizing problems at the age of three years. Of interest, more Enterobacteriaceae species were cross-sectionally related to decreased cognitive functioning at the age of 45 months (Streit et al., 2021).

Another of our findings was that higher alpha diversity at age two weeks was linked to fewer internalizing problems and better executive functions at age three years. In accordance with our internalizing behavior result, Laue et al., observed that higher alpha diversity in the first two months of life was related to less internalizing behavior in three-year-old boys (Laue et al., 2021). Furthermore, van de Wouw et al., found lower alpha diversity in three-to-five-year-old children with clinically relevant CBCL cutoff scores for internalizing behavior (Van De Wouw et al., 2022). Besides, Eckermann et al., observed higher Shannon alpha diversity at the age of one, three, and four months in relation to better cognitive ability as measured by Digit Span forwards test at the age of ten years, although not for Shannon alpha diversity at the age of six and ten years (Eckermann et al., 2022). Despite the generally weak or absent relations found between child gut microbiota and executive functions, the abovementioned studies may indicate that higher alpha diversity in the first years of life is related to improved subsequent mental outcomes at later ages. On the contrary, Carlson et al., found higher alpha diversity at one year of age related to worse cognition at two years of age in typically developing toddlers (Carlson et al., 2017). Additionally, a recent paper from van de Wouw et al. (2023) found evidence for a weak to modest cross-sectional relation between higher Shannon diversity and worse verbal comprehension in three- to four-year-old children (van de Wouw et al., 2023). Together, these findings illustrate the fact that we are yet to understand the potential impact of alpha diversity levels at different developmental stages. In general, alpha diversity levels of newborns start increasing immediately after birth due to colonization of microorganisms. With time, breastfed-infants tend to form a Bifidobacterium-predominated configuration which is often less diverse than formula-fed infants. After the introduction of solid food, child alpha diversity starts to increase, gradually reaching a steady state resembling gut microbiota composition of adults. Given these apparent normative fluctuations in levels of alpha diversity in the first months and years of life, having a comparatively high (or low) alpha diversity will potentially impact a child’s development differently depending on the child’s specific age.

Our study contributes to the growing body of literature on the gut microbiota, problem behavior, and executive functions. A strength of our study is the longitudinal design, which covered the period from birth to age three years and allowed for the assessment of multiple developmental stages of the gut microbiota. Another advantage was that questionnaires were filled in by both mother and partner. Partner reports were used for sensitivity analyses and because they showed positive correlations with maternal reports, they enhanced the reliability of our study measures. Furthermore, problem behavior and executive functions were assessed with multiple questionnaires (i.e., CBCL and SDQ for problem behavior, and BRIEF-P and REEF for executive functions), allowing us to determine conformity and consistency between various questionnaires. Finally, we used standardized behavioral tasks as a tool to objectively determine child executive functions.

A limitation of our study is the possible overreliance on the compositional features of the gut microbiota using relative abundances instead of absolute abundances. This approach may increase the chance of spurious associations as relative abundances are dependent on each other. Besides, 16S rRNA gene sequence data are limited at species-level resolution and profiling precise gene functions (Durazzi et al., 2021). Hence, although multiple associations were identified, it is worth noting that these relations do not indicate causality. The results await follow-up studies, preferably preclinical experimental studies, to determine if individual microbes (e.g., aforementioned Streptococcus and Bifidobacterium) and microbial communities as a whole influence behavior and cognition, e.g., by generating neurotransmitters (e.g., GABA and serotonin) and their precursors (e.g., tryptophan and phenylalanine) in the gut (Altaib et al., 2021; Barandouzi et al., 2022; Biederman & Spencer, 1999; Gizer et al., 2009; Kandel et al., 2000; Kanehisa et al., 2022; Staller & Faraone, 2007). Further research including quantitative PCR, whole-genome shotgun metagenomic sequencing, targeted fecal metabolomics, and experimental studies in animal models, would improve the understanding of current correlational results and provide insight into microbial functions and even causality. An analytic limitation of our study is LEfSe, which was used to identify differentially abundant taxa in this study and has recently been pointed out to have higher sensitivity to false positive rates compared to other microbial composition analyses, such as ALDEx2, ANCOM − II, and DESeq2 (Nearing et al., 2022). Due to such methodological limitations, LEfSe-based significant findings should be carefully validated in future studies. Another limitation of our study is the relatively small sample size and mostly highly educated study population, limiting the generalizability of the findings. The restricted sample size may also hamper deep inference with respect to taxa with low prevalence rates to some degree. Our findings on microbial relations to the mental outcomes need to be confirmed in a larger, more representative cohort. Also, the participants were highly educated, which may hamper translating of our findings to individuals with a lower socioeconomic status. Lastly, recent developments in longitudinal analytical approaches (Kodikara et al., 2022), such as zero-inflated beta regression models, block bootstrap methods, and SplinectomeR, will better facilitate the identification of differentially abundant microbial taxa between groups (e.g., below and above clinical cutoffs) over time.

To conclude, our results provide tentative evidence supporting the idea that in a child’s first years of life, the gut microbiota might play a vital role in the development of the brain, in line with the early-life programing theory (Tarry-Adkins & Ozanne, 2011). Potential mechanisms are likely related to microbiota-derived metabolites (Ahmed et al., 2022). As the nature of this study was exploratory and the body of similar research needs to grow to a large extent, it is still premature to translate our correlational findings into clinical implications. Replications in other longitudinal studies on healthy community children are necessary to confirm our findings. Ideally, to avoid inconsistent results caused by different methods used, replication studies should apply the same methodology regarding microbiology, genomics, epidemiology, and statistics (Ou et al., 2024). This will shed more light on key microbial taxa and latent pathways of associations between early gut microbiota and child behavior and cognition.